I'm trying to implement a recursive function using integrate, and I
suspect I need a Vectorize somewhere,
but I can't suss it out. Any help would be appreciated. I've tried
traceback() and various debugging ideas to no avail (most likely due to
my inexperience with these tools.)
Here's what I have.
Nk <- function(m, C) {
if (length(m) > 1) {
rho <- C[1, -1]
Rmat <- C[-1, -1]
B <- diag(1/sqrt(1 - rho*rho)) %*%
(-rho %*% t(rho) + Rmat) %*%
diag(1/sqrt(1 - rho*rho))
integrate( function(x) dnorm(x) * Nk((m[-1] - rho*x)/sqrt(1 -
rho*rho), B), -10, m[1] )$value
} else {
pnorm(m[1])
}
}
my example is> x2 <- c(0.0781292, -1.6403152)
> sigma2 <- matrix(c(1, -0.5054781, -0.5054781, 1), nrow=2)
> Nk(x2, sigma2)
Error in integrate(function(x) dnorm(x) * Nk((m[-1] - rho * x)/sqrt(1 -
:
non-finite function value
All the pieces work outside of the function, and the integrand is finite
as far as I can see.
[Yes, this is a recursive function for multivariate cumulative normal.
It seems to match (so far for k=2 without recursion) the existing R
functions from packages mvtnorm and mnormt.
It is from D. Cassimon, et al. Closed-form valuation of American call
options on stocks paying multiple dividends. Finance Research Letters 4
(2007) 33-48.]
Thank you to anyone who can shed some light on this.
David L. Reiner, PhD
Rho Trading Securities, LLC
davidr at rhotrading.com
2008-Feb-20 19:54 UTC
[R] recursive function help SOLVED (sort of)
Well, it turns out to be very simple - just insert a Vectorize between
integrate and function(x).
However, the special cases where C[i,j]==1 make the actual code quite
messy.
It matches pmvnorm {mvtnorm} and pmnorm {mnormt} quite well.
And the recursive method is incredibly slow for higher dimensions.
And still some cases blow up or don't converge.
So, never mind. It looked clever, but I would recommend pmvnorm for
speed and accuracy,
even though it is non-deterministic for higher dimensions.
One little note: for the bivariate, this method (without recursion) is
as accurate as the existing methods
and a bit faster than pmvnorm.
-- David
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of davidr at rhotrading.com
Sent: Tuesday, February 19, 2008 3:08 PM
To: r-help at r-project.org
Subject: [R] recursive function help
I'm trying to implement a recursive function using integrate, and I
suspect I need a Vectorize somewhere,
but I can't suss it out. Any help would be appreciated. I've tried
traceback() and various debugging ideas to no avail (most likely due to
my inexperience with these tools.)
Here's what I have.
Nk <- function(m, C) {
if (length(m) > 1) {
rho <- C[1, -1]
Rmat <- C[-1, -1]
B <- diag(1/sqrt(1 - rho*rho)) %*%
(-rho %*% t(rho) + Rmat) %*%
diag(1/sqrt(1 - rho*rho))
integrate( function(x) dnorm(x) * Nk((m[-1] - rho*x)/sqrt(1 -
rho*rho), B), -10, m[1] )$value
} else {
pnorm(m[1])
}
}
my example is> x2 <- c(0.0781292, -1.6403152)
> sigma2 <- matrix(c(1, -0.5054781, -0.5054781, 1), nrow=2)
> Nk(x2, sigma2)
Error in integrate(function(x) dnorm(x) * Nk((m[-1] - rho * x)/sqrt(1 -
:
non-finite function value
All the pieces work outside of the function, and the integrand is finite
as far as I can see.
[Yes, this is a recursive function for multivariate cumulative normal.
It seems to match (so far for k=2 without recursion) the existing R
functions from packages mvtnorm and mnormt.
It is from D. Cassimon, et al. Closed-form valuation of American call
options on stocks paying multiple dividends. Finance Research Letters 4
(2007) 33-48.]
Thank you to anyone who can shed some light on this.
David L. Reiner, PhD
Rho Trading Securities, LLC
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.