Hello all,
I am trying to produce the closed-form Ait-Sahalia approximation to the
log-likelihood function of a diffusion process but the function arguments I
am passing to the function result in NaN being returned. I've checked the
'Transform','S' and 'M' functions and the problem doesn't seem to be with
these. Does anyone have any ideas why this isn't working? My code is
reproduced below for info.
Thanks in advance.
Regards,
Steven
## simulate sample of diffusion process to be estimated
# define drift and diffusion coefficients for SDE
d <- expression(x^(-1) - 1 + x - x^2)
s <- expression(x^1.3)
s.x <- expression(1.3*(x^0.3))
# simulate process using the Milstein approximation
X <-
sde.sim(t0=0,T=100,X0=1,N=1000,drift=d,sigma=s,sigma.x=s.x,method="milstein")
## define the Lamperti transform function
Transform <- function(t,x,theta)
(1/(theta[5]*(theta[6]-1)))*(x^(1-theta[6]))
## define the diffusion function
S <- function(t,x,theta) theta[5]*(x^theta[6])
## define the transformed drift function and it's first six derivatives
m0 <- function(t,x,theta)
{
b <- theta[5]*(theta[6]-1)
p1 <- (theta[6]+1)/(theta[6]-1)
p2 <- theta[6]/(theta[6]-1)
p3 <- (2-theta[6])/(1-theta[6])
(-1/theta[5])*(theta[1]*((b*x)^p1) - theta[2]*((b*x)^p2) + theta[3]*b*x -
theta[4]*((b*x)^p3) - 0.5*theta[6]*(theta[5]^2)*((b*x)^(-1)))
}
m1 <- function(t,x,theta)
{
b <- theta[5]*(theta[6]-1)
p1 <- (theta[6]+1)/(theta[6]-1)
p2 <- theta[6]/(theta[6]-1)
p3 <- (2-theta[6])/(1-theta[6])
(-b/theta[5])*(theta[1]*(p1*(b*x)^(p1-1)) - p2*theta[2]*((b*x)^(p2-1)) +
theta[3] - p3*theta[4]*((b*x)^(p3-1)) +
0.5*theta[6]*(theta[5]^2)*((b*x)^(-2)))
}
m2 <- function(t,x,theta)
{
b <- theta[5]*(theta[6]-1)
p1 <- (theta[6]+1)/(theta[6]-1)
p2 <- theta[6]/(theta[6]-1)
p3 <- (2-theta[6])/(1-theta[6])
(-(b^2)/theta[5])*(theta[1]*(p1*(p1-1)*(b*x)^(p1-2)) -
p2*(p2-1)*theta[2]*((b*x)^(p2-2)) - p3*(p3-1)*theta[4]*((b*x)^(p3-2)) -
theta[6]*(theta[5]^2)*((b*x)^(-3)))
}
m3 <- function(t,x,theta)
{
b <- theta[5]*(theta[6]-1)
p1 <- (theta[6]+1)/(theta[6]-1)
p2 <- theta[6]/(theta[6]-1)
p3 <- (2-theta[6])/(1-theta[6])
(-(b^3)/theta[5])*(theta[1]*(p1*(p1-1)*(p1-2)*(b*x)^(p1-3)) -
p2*(p2-1)*(p2-2)*theta[2]*((b*x)^(p2-3)) -
p3*(p3-1)*(p3-2)*theta[4]*((b*x)^(p3-3)) +
3*theta[6]*(theta[5]^2)*((b*x)^(-4)))
}
m4 <- function(t,x,theta)
{
b <- theta[5]*(theta[6]-1)
p1 <- (theta[6]+1)/(theta[6]-1)
p2 <- theta[6]/(theta[6]-1)
p3 <- (2-theta[6])/(1-theta[6])
(-(b^4)/theta[5])*(theta[1]*(p1*(p1-1)*(p1-2)*(p1-3)*(b*x)^(p1-4)) -
p2*(p2-1)*(p2-2)*(p2-3)*theta[2]*((b*x)^(p2-4)) -
p3*(p3-1)*(p3-2)*(p3-3)*theta[4]*((b*x)^(p3-4)) -
12*theta[6]*(theta[5]^2)*((b*x)^(-5)))
}
m5 <- function(t,x,theta)
{
b <- theta[5]*(theta[6]-1)
p1 <- (theta[6]+1)/(theta[6]-1)
p2 <- theta[6]/(theta[6]-1)
p3 <- (2-theta[6])/(1-theta[6])
(-(b^5)/theta[5])*(theta[1]*(p1*(p1-1)*(p1-2)*(p1-3)*(p1-4)*(b*x)^(p1-5)) -
p2*(p2-1)*(p2-2)*(p2-3)*(p2-4)*theta[2]*((b*x)^(p2-5)) -
p3*(p3-1)*(p3-2)*(p3-3)*(p3-4)*theta[4]*((b*x)^(p3-5)) +
60*theta[6]*(theta[5]^2)*((b*x)^(-6)))
}
m6 <- function(t,x,theta)
{
b <- theta[5]*(theta[6]-1)
p1 <- (theta[6]+1)/(theta[6]-1)
p2 <- theta[6]/(theta[6]-1)
p3 <- (2-theta[6])/(1-theta[6])
(-(b^6)/theta[5])*(theta[1]*(p1*(p1-1)*(p1-2)*(p1-3)*(p1-4)*(p1-5)*(b*x)^(p1-6))
- p2*(p2-1)*(p2-2)*(p2-3)*(p2-4)*(p2-5)*theta[2]*((b*x)^(p2-6)) -
p3*(p3-1)*(p3-2)*(p3-3)*(p3-4)*(p3-5)*theta[4]*((b*x)^(p3-6)) -
360*theta[6]*(theta[5]^2)*((b*x)^(-7)))
}
# now create list consisting of m0..m6
M <- list(m0,m1,m2,m3,m4,m5,m6)
## use HPloglik function to calculate log-likelihood function at the true
parameter values
HPloglik(X,c(1,1,1,1,1,1.3),Moments,Transform,S)
--
S
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
