Indefinite Hessian in optim() and disabled "start" argument in mle()

238 views Asked by At

I was trying to use mle() to minimize the -log likelihood function, however it always pops out an error message like this:

Error in minuslogl() (from kf.r@7524iUV#2) : argument "psi" is missing, with
no default

My functions is:

KFlogL2 <- function(psi){
.........
return (res)
}

> psi.test
[1] 2.00 0.20 0.20 0.02 0.20 0.20 0.50

> mle(KFlogL2,start = list(psi = psi.test),method = 'L-BFGS-B',lower =
> lb,upper = ub)

> mle(KFlogL2,start = list(psi = c(2,.2,.2,.02,.2,.2,.5)),method =
> 'L-BFGS-B',lower = lb,upper = ub)

I tried both ways above but neither works...

Anyone knows why?

For those who would like to introduce optim(), I've actually tried that, however my Hessian matrix, after taking inverse and extracting diagonals - turns to be negative! The convergence is fine. I have no idea what that means. So I wanted to try with mle() to see if things can change. My Optimization function:

KFopt <- optim(par = psi.test, fn = KFlogL, method = 'L-BFGS-B',
               lower = lb, upper = ub, hessian = TRUE)

My KF function:

KFlogL2 <- function(psi){
    k <- psi[1]
    sigmax <- psi[2]
    lambdax <- psi[3]
    mu <- psi[4]
    sigmae <- psi[5]
    rnmu <- psi[6]
    pxe <- psi[7]
    m <- length(init.state)
    N <- ncol(y.use)
    nobs <- nrow(y.use)
    s <- rep(.01,N)
    dt <- 7/360
    cc <- c(0,mu * dt)
    T <- diag(m)
    T[1,1] <- exp(-k * dt)
    xx <- (1-exp(-2 * k * dt)) * sigmax * sigmax / 2 / k
    xy <- (1-exp(-k * dt)) * pxe * sigmax * sigmae / k
    yx <- (1-exp(-k * dt)) * pxe * sigmax * sigmae / k
    yy <- sigmae * sigmae * dt
    Q <- matrix(c(xx,xy,yx,yy),m,m)
    R <- diag(m)
    # measurement equation
    d <- rep(0,N)
    Z <- matrix(0,N,m)
    for (i in 1:N){
        matur.i <- matur.use[i]
        p1 <- (1-exp(-2 * k * matur.i)) * sigmax * sigmax / 2 / k
        p2 <- sigmae * sigmae * matur.i
        p3 <- 2 * (1-exp(k * matur.i)) * pxe * sigmax * sigmae / k
        d[i] <- rnmu * matur.i - (1-exp(-k * matur.i)) * lambdax / k + 1/2 * (p1 + p2 + p3)
        Z[i,] <- exp(-k * matur.i)
    }
    H <- diag(s)
    # Kalman Filter
    save.ytt1 <- save.vtt <- save.vt <- matrix(0,nobs,N)
    save.att1 <- save.att <- matrix(0,nobs,m)
    save.Ptt1 <- save.Ptt <- matrix(0,nobs,m*m)
    save.Ftt1 <- matrix(0,nobs,N*N)
    save.detF <- save.vFv <- rep(0,nobs)
    Ptt <- init.dist
    att <- init.state
    for (iter in 1:nobs){
        Ptt1 <- T %*% Ptt %*% t(T) + R %*% Q %*% t(R)
        Ftt1 <- Z %*% Ptt1 %*% t(Z) + H
        att1 <- T %*% att + cc
        yt <- y.use[iter,]
        ytt1 <- Z %*% att1 + d
        vt <- yt - ytt1
        temp <- Ptt1 %*% crossprod(Z,solve(Ftt1))
        att <- att1 + temp %*% vt
        Ptt <- Ptt1 - temp %*% Z %*% Ptt1
        ytt <- Z %*% att + d
        vtt <- yt - ytt
        save.vtt[iter,] <- vtt
        save.vt[iter,] <- vt
        save.att[iter,] <- att
        save.Ptt1[iter,] <- as.numeric(Ptt1)
        save.Ptt[iter,] <- as.numeric(Ptt)
        save.detF[iter] <- det(Ftt1)
        save.vFv[iter] <- crossprod(vt,solve(Ftt1,vt))
    }
    logL <- - (N + 1) * nobs/2 * log(2 * pi) - 1/2 * sum(log(save.detF)) - 1/2 * sum(save.vFv)
    res <- -logL
    return (res)
}

I don't know how to attach a file here, so here is just a glimpse:

3.130263167 3.130700134 3.058707073 3.012589391 2.999724295 2.991724252
3.102342009 3.09421922  2.999724295 2.952824773 2.940747965 2.932259851
3.118392286 3.125882958 3.006177531 2.949164638 2.926917958 2.913979772
3.106378794 3.072693315 2.991724252 2.950211758 2.937573359 2.926917958
3.104138147 3.111735949 3.031099417 2.991724252 2.979602892 2.969388298
3.107273648 3.113959655 3.033509638 2.992226134 2.979602892 2.969388298
3.085572978 3.086486637 3.059176446 3.033509638 3.026746327 3.020912572
3.100092289 3.099641737 3.064325065 3.033028058 3.025776395 3.019936962
3.073618812 3.07176696  3.061520014 3.054001182 3.052112607 3.050220459
3.058707073 2.961658293 3.054001182 3.039270576 3.03543364  3.032546247
3.003204288 3.006672214 3.016024977 3.007660844 3.003700443 3.001217204
2.96217549  2.959068289 3.013080912 3.017493765 3.017493765 3.017493765
3.013572192 3.0194488   3.022860941 3.014554028 3.006672214 3.004692015
3.007166651 3.010620886 3.042138646 3.038791763 3.033028058 3.0301337
2.865623588 2.872434057 2.972463647 2.994731773 2.998727783 3.002707887
2.85474458  2.858766418 2.983659692 2.989714201 2.987700102 2.985681938
2.882003508 2.944438979 2.999226163 2.999226163 2.99723115  2.995232149
2.931726944 2.938103161 3.012097628 3.006672214 3.009635179 3.00914196
2.904712875 2.905807566 2.983659692 2.984165637 2.983659692 2.983659692
2.971439581 2.975019232 3.007660844 3.001217204 2.996731774 2.993730271
2.862772146 2.872434057 2.993229143 2.998727783 2.994731773 2.991724252
2.890371758 2.889260029 2.992727765 3.005187432 3.003700443 3.003700443
2.797890905 2.814210397 2.930126516 2.959586827 2.97603964  2.986691529
2.855895328 2.862772146 2.947591898 2.965788397 2.973997781 2.980110893
2.744060639 2.750470917 2.913437031 2.949164638 2.964241606 2.975019232
2.840247371 2.841414913 2.936512914 2.962692419 2.97603964  2.986186861
2.817203515 2.821378886 2.915606229 2.935982269 2.947067102 2.957511061
2.836150204 2.829677689 2.938103161 2.961140829 2.973997781 2.986186861
2.903068589 2.925846146 3.032546247 3.039749159 3.043569603 3.045474365
2.987700102 3.023347441 3.060583246 3.060114532 3.060114532 3.060114532
3.023833704 3.029650492 3.085572978 3.079153882 3.078693794 3.078693794
3.387774361 3.343215099 3.273742726 3.23317313  3.201119103 3.168003494
3.284663565 3.274121299 3.221272949 3.185939325 3.16463081  3.147594623
3.348499593 3.346741196 3.298426104 3.244153632 3.212455257 3.195811885
3.32251486  3.327909586 3.24998682  3.191710157 3.162940193 3.144152279
3.377587516 3.371425223 3.291010423 3.224857897 3.191710157 3.173459961
3.410817625 3.426215145 3.308716529 3.229222117 3.193763124 3.176803048
3.510948246 3.508855256 3.390473418 3.300640127 3.259249719 3.240245851
3.646232879 3.625140613 3.476614021 3.354455119 3.292126287 3.258865473
3.538928277 3.524888854 3.454738149 3.328626689 3.249211025 3.206803244
3.706964922 3.698829785 3.541828511 3.400863993 3.316365446 3.272606147
3.672241813 3.660737148 3.507357577 3.375537635 3.296207168 3.25617161
3.399529325 3.379973745 3.237501289 3.141994781 3.101442728 3.075928816
3.553346059 3.542118073 3.380994674 3.270329106 3.213662258 3.176803048
3.478467017 3.485232111 3.325395668 3.223266173 3.165475048 3.127199036
3.505557397 3.50013733  3.348850901 3.247268899 3.195402469 3.162093811
3.384390263 3.364187556 3.256556892 3.185112195 3.143289838 3.118834471
3.498021566 3.492256113 3.348148161 3.259633817 3.210843653 3.181796817
3.369018483 3.422958873 3.295466427 3.211246798 3.171364842 3.148024084
3.277144733 3.273742726 3.169265324 3.113515309 3.097837496 3.08876714
0

There are 0 answers