[R] DLSODA error

From: Kristian Lind <kristian.langgaard.lind_at_gmail.com>
Date: Thu, 28 Apr 2011 11:17:36 +0200


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_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. Received on Thu 28 Apr 2011 - 09:33:59 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 28 Apr 2011 - 21:40:34 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive