Dear R-users,
I'm running an MLE procedure where some ODEs are solved for each iteration
in the maximization process. I use mle2 for the Maximum likelihood and
deSolve for the ODEs.
The problem is that somewhere along the way the ODE solver crashes and I get
the following error message:
DLSODA- Warning..Internal T (=R1) and H (=R2)
are
such that in the machine, T + H = T on the next
step
(H = step size). Solver will continue
anyway.
In above, R1 = 0.2630651367072D+01 R2 = 0.2055665829864D-15
....
Error in optim(par = par, fn = U, method = "Nelder-Mead", control =
list(maxit = 100), :
function cannot be evaluated at initial parameters
...
Error in BBsolve(par = par, fn = Bo, s = s, outmat = outmat, method = c(1,
:
object 'ans.best' not found
In addition: Warning messages:
1: In lsoda(y, times, func, parms, ...) :
an excessive amount of work (> maxsteps ) was done, but integration was
not successful - increase maxsteps
2: In lsoda(y, times, func, parms, ...) :
Returning early. Results are accurate, as far as they go
For each iteration of the MLE my code gives me the parameters. If I use the
parameters from the last iteration before the above error message appears I
would expect that it should crash immediately, but this isn't the case. Any
suggestions why this is the case and how to avoid the ODE solver to crash
would be much appreciated.
Thank you for your time.
Kristian
the code is below.
#Loading the needed packages
library(deSolve)
library(stats4)
library(BB)
library(bbmle)
c2 <- c(3.02, 2.88, 2.91, 2.83, 2.85, 2.88, 2.91, 2.90, 2.94, 3.09, 3.17,
3.14, 3.37, 3.40, 3.50, 3.58, 3.55, 3.70, 3.90, 3.77)
c10 <- c(4.39, 4.22, 4.27, 4.21, 4.25, 4.34, 4.47, 4.40, 4.46, 4.64, 4.73,
4.60, 4.87, 4.96, 5.09, 5.10, 5.08, 5.26, 5.54, 5.37)
T = 20
#definining -log-likelihood function
minusloglik <- function(K_vv, K_rv, K_rr, theta_v, theta_r, Sigma_rv,
Sigma_rr, lambda_v, lambda_r){
# #solving ODEs as functions of "parameters"
parameters <- c(K_vv,
K_rv,
K_rr,
theta_v,
theta_r,
Sigma_rv,
Sigma_rr,
lambda_v,
lambda_r)
state <- c( b_1 = 0,
b_2 = 0,
a = 0)
#declaring ODEs
Kristian <- function(t, state, parameters){
with(as.list(c(state, parameters)),
{
db_1 = -((K_vv+lambda_v)*b_1+(K_rv+Sigma_rv*lambda_v
+Sigma_rr*lambda_r)*b_2+0.5*(b_1)^2+Sigma_rv*b_1*b_2+0.5*
((Sigma_rv)^2+(Sigma_rr)^2)*(b_2)^2 )
db_2 = -K_rr*b_2+1
da = K_vv*theta_v*b_1+(K_rv*theta_v+K_rr*theta_r)*b_2
list(c(db_1, db_2, da))
})
}
# time making a sequence from t to T evaluated at each delta seq(t, T,
by = delta)
times <- seq(0, 10, by = 0.5)
outmat <- ode(y = state, times = times, func = Kristian, parms =
parameters)
#solving to equations in to unknowns.
Bo <- function(x, s, outmat)
{
f <- rep(NA, length(x))
z <- - outmat[,2:4] %*% c(x[1],x[2],1)
f[1] <- (1-exp(z[5]))/sum(exp(z[2:5])) - s[1]
f[2] <- (1-exp(z[21]))/sum(exp(z[2:21])) - s[2]
f
}
#extracting statevariables
v<- numeric(T)
r<- numeric(T)
for(i in 1:T)
{
s <- c(c2[i],c10[i] )
par <- c(50, 5)
BBsolveans <- BBsolve(par=par, fn=Bo, s=s, outmat=outmat, method
=c(1,2,3), control =list(maxit =1000), quiet =TRUE)
v[i] <- BBsolveans$par[1]
r[i] <- BBsolveans$par[2]
}
r <- r[!v < 0]
v <- v[!v < 0]
#calculating transition density as a function of parameters
transdens <- function(v1, v2, r1, r2, t1,t2)#where v1 is the volatily at
time i and v2 is at time i+1
{
delta <- 7/365
alpha_r <- 0
c <- 2*K_vv*(1-exp(-K_vv*delta))^(-1)
q <- 2*K_vv*theta_v-1
m <-
-0.5*(v2*(K_rv*v2-2*K_rv*theta_v+K_rr*theta_r)-v1*(K_rv*v1-2*K_rv*theta_v+K_rr*theta_r))+exp(-K_rr)*r1
var <-
abs(0.5*((v2*(2*alpha_r+(Sigma_rr+Sigma_rv)^2*v2)-(v1*(2*alpha_r+(Sigma_rr+Sigma_rv)^2*v1)))))
f <- rep(NA, 1)
f[1]<-pchisq(2*c*v2, 2*q+2, ncp=2*c*v1*exp(-K_vv*delta), lower.tail =
TRUE, log.p = FALSE)*pnorm(r2,mean = m, sd= sqrt(var))
f
}
p<- numeric(length(v)-1)
for(i in 1:(length(v)-1)){
p[i] <- transdens(v1=v[i], v2=v[i+1], r1=r[i], r2=r[i+1], t1=i,
t2=i+1)
}
# Calculating Jacobian determinant.
gfunc <- function(x)
{
z <- - outmat[,2:4] %*% c(x[1],x[2],1)
x[1] <- x1
x[2] <- x2
c((1-exp(z[5]))/sum(exp(z[2:5])), (1-exp(z[21]))/sum(exp(z[2:21])))
}
J<- numeric(length(v)-1)
for(i in 1:(length(v)-1)){
x1 <- v[i]
x2 <- r[i]
x <- cbind(x1, x2)
Jac <- jacobian(func=gfunc, x=x)
J[i] <- abs(det(Jac))
}
J <- J[!p == 0]
p <- p[!p == 0]
#removing NAs
J <- J[!is.na(J)]
p <- p[!is.na(p)]
print(parameters)
f<-rep(NA, 1)
f[1] <- -1/(T-1)*sum(log(p)-log(J)) # the 20 should be equal to the length
of the sample.
f
}
# process crashes with DLSODA error for these parameter values:
fit <- mle2(minusloglik, start = list(K_vv =0.523023048, K_rv=-0.421004051,
K_rr = 0.082203320 , theta_v = 0.373992355, theta_r=0.290026071,
Sigma_rv=-0.655919272, Sigma_rr=-0.072001589,lambda_v= 0.006053411,
lambda_r=0.607157538), method = "L-BFGS-B")
#but can be restarted for these only to crash later. Notice that it is only
lambda_r that is different in the two.
#fit <- mle2(minusloglik, start = list(K_vv =0.523023048, K_rv=-0.421004051,
K_rr = 0.082203320 , theta_v = 0.373992355, theta_r=0.290026071,
Sigma_rv=-0.655919272, Sigma_rr= -0.072001589,lambda_v= 0.006053411,
lambda_r=0.606157538), method = "L-BFGS-B")
print(summary(fit))
[[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.
