Hello All,
I am trying to estimate the parameters of a stochastic differential equation
(SDE) using quasi-maximum likelihood methods but I am having trouble with
the 'optim' function that I am using to optimise the log-likelihood
function.
After simulating the SDE I generated samples of the simulated data of
varying size (I want to see what effect adding more observations has on the
accuracy of the estimates) and ran the 'optim' function to optimise the
log-likelihood function. The optimiser seemed to work fine for sample sizes
up to 100,000 i.e. the optim function converged, but when I tried to run
optim using a sample size of 200,000 the optim failed to converge - the
correspoding error message was "NEW_X". Can someone please advise what this
error message means, and what I can do to avoid producing this error? I
have supplied an extract of my code below for inspection (note that I am
fairly new to programming so my code is probably not the most
efficient/elegant...).
Thank you in advance,
Steven
# I use a Milstein scheme to generate a simulation of the SDE I am trying to
estimate
Mil.sim <- function(X0,drift,diffusion,diffusion.x,horizon,no.steps)
{
set.seed(123)
X <- rep(NA,no.steps+1) # initialise vector to store sample
X[1] <- X0
shocks <- rnorm(no.steps) # generate normal shocks for FEM method
step.size <- horizon/no.steps # define step size to use in simulation
for (i in 2:(no.steps+1)) # loop to generate sample
{
X[i] <- X[i-1] + drift(X[i-1])*step.size +
diffusion(X[i-1])*sqrt(step.size)*shocks[i-1] +
0.5*diffusion(X[i-1])*diffusion.x(X[i-1])*step.size*(shocks[i-1]^2 - 1)
}
X <<- ts(X,start=0,end=horizon,frequency=1/step.size) # coerce
simulated process into 'ts' class for ease of plotting
}
# generate simulation of SDE
Mil.sim(1,d,s,s.x,1000,10000000)
# this sets up the sample and associated variables that I will use to carry
out the estimation
observations.1 <- seq(1,by=50,length(X)) # set up sample of X,
varying 'by=..' varies the size (and spacing) of the sample
X.obs.1 <- X[observations.1]
obs.times.1 <- (observations.1-1)/10000
delta.1 <- obs.times.1[2] - obs.times.1[1]
sample.1 <- data.frame(tt = obs.times.1,values = X.obs.1) # create
data frame to be used in optimisation of log likelihood function
n.1 <- length(sample.1$values)-1
# define the log-likelihood function for optimisation
IEM.loglik.fn <- function(par,data)
{
A <- data[2:(n.1+1)]*(1-par[3]*delta.1) -
par[1]*delta.1*(data[2:(n.1+1)]^(-1)) + par[2]*delta.1 +
par[4]*delta.1*(data[2:(n.1+1)]^2) - data[1:n.1]
B <- data[1:n.1]^(-2*par[6])
C <- 1 - par[3]*delta.1 + par[1]*delta.1*(data[2:(n.1+1)]^(-2)) +
2*par[4]*delta.1*data[2:(n.1+1)]
loglik <- -(n.1/2)*log(2*pi*delta.1) - n.1*log(par[5]) -
par[6]*sum(log(data[1:n.1])) - (1/(2*delta.1*(par[5]^2)))*sum(B*(A^2)) +
sum(log(C))
}
# define the gradient of the parameters to be estimated (6 parameters to
estimate)
IEM.gradient <- function(par,data)
{
A <- data[2:(n.1+1)]*(1-par[3]*delta.1) -
par[1]*delta.1*(data[2:(n.1+1)]^(-1)) + par[2]*delta.1 +
par[4]*delta.1*(data[2:(n.1+1)]^2) - data[1:n.1]
B <- data[1:n.1]^(-2*par[6])
C <- 1 - par[3]*delta.1 + par[1]*delta.1*(data[2:(n.1+1)]^(-2)) +
2*par[4]*delta.1*data[2:(n.1+1)]
diff.1 <- (1/(par[5]^2))*sum((data[2:(n.1+1)]^(-1))*A*B) +
delta.1*sum((data[2:(n.1+1)]^(-2))*(C^(-1)))
diff.2 <- -(1/(par[5]^2))*sum(A*B)
diff.3 <- (1/(par[5]^2))*sum(data[2:(n.1+1)]*A*B) - delta.1*sum(C^(-1))
diff.4 <- -(1/(par[5]^2))*sum((data[2:(n.1+1)]^2)*A*B) +
2*delta.1*sum(data[2:(n.1+1)]*(C^(-1)))
diff.5 <- -(n.1/par[5]) + (1/(delta.1*(par[5]^3)))*sum((A^2)*B)
diff.6 <- -sum(log(data[1:n.1])) +
(1/(delta.1*(par[5]^2)))*sum((log(data[1:n.1]))*(A^2)*B)
obj <- c(diff.1,diff.2,diff.3,diff.4,diff.5,diff.6)
return(obj)
}
# use optim function to find maximum likelihood parameter estimates
opt <-
optim(c(3,3,3,3,3,3),method="L-BFGS-B",lower=c(0,0,0,0,1e-8,0),fn=IEM.loglik.fn,gr=IEM.gradient,control=list(fnscale=-1),hessian=TRUE,data=sample.1$values)
--
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.
