[R] Solve in maximum likelihood estimation

Harry Ho hassanysabbah at hotmail.com
Sat Feb 17 20:18:58 CET 2007


Hi,

I got the following problem.

I am doing a maximum likelihood estimation for a Kalman Filter.

For this purpose, I have to invert an error matrix Ffast of dimension
"no. parameters X no.parameters". The usualy optim methods often find only 
local minima, so I decided to make the optimization using the SANN 
algorithm, which is very slow already.

However, this becomes a real problem because the "reciprocal condition 
number" of Ffast becomes extremely small (up to 1e-1000 for the amount of 
data I have) for non-sensible paramter combinations.
Thus, I need to extend the tolerance of the solve algorithm to this level 
(tol=1e-1000), which means that calculations can become incredibly slow for 
many parameters.$

Is there a way to address this issue? I could imagine replacing the matrix 
by a dummy matrix each time I get the singularity error, hoping that the 
program recognizes the implausibility. Is there a function that allows doing 
so?

Another way might be, as the temperature decreases, to decrease the 
tolerance as well, as estimates become more reasonable (would have to verify 
that).
How could that work?

I hope this characterization is enough. I still provided the code used, just 
in case somebody bothers to take a closer look at it performance wise.

It certainly has a lot of bad programming style though.


Nevertheless, thanks a lot for all replies.

Remarks: n is very small as to keep the code working in an acceptable time 
frame.
Up to the +-line, the data are generated. After that, the Kalman Filter 
estimates them.
I have restricted this estimation to 4 parameters, the true values are given 
for the other ones.
The estimation results in this form are pretty bad, but with more 
observations, maturities and above all computation they get quite good.









set.seed(13234)

n           <- 20
matur       <- c(3,5,12)


no.state            <- 2
no.obs              <- length(matur)

Xav                 <- matrix(  nrow=no.state,ncol=n)
FF <- FF.kal        <- matrix(  nrow=no.state,ncol=no.state)

GG <- GG.kal        <- matrix(  nrow=no.obs,ncol=no.state)

V                   <- matrix(0,nrow=no.state)
W                   <- matrix(0,nrow=no.obs)

SigV <- SigV.kal    <- matrix(0,ncol=no.state,nrow=no.state)
SigW <- SigW.kal    <- matrix(0,nrow=no.obs,ncol=no.obs)

x                   <- rep(0,n)
A <- B  <- G.B      <- numeric(matur[no.obs])
Z                   <- matrix(ncol=n,nrow=no.obs)


kappa       <- 0.97
theta       <- 0.07
sigma       <- 0.005
lambda      <- -0.162
beta        <- sigma*lambda
rho         <- 0.5*sigma^2*lambda^2

Q           <- sigma^2
R           <- sigma^2*0.125



rho         <- 0.5*lambda^2*sigma^2


Xav[1,]               <- rep(1,n)

FF[1,] <- FF.kal[1,]<- c(1,rep(0,no.state-1))
FF[2,]              <- c(theta-theta*kappa,kappa)

diag(SigW)          <- R
SigV[2,2]           <- Q

for (i in 1:matur[length(matur)]){
B[i] <- (1-kappa^(i-1))/(1-kappa)
G.B[i] <- rho+B[i]*theta*(1-kappa)-0.5*beta^2-B[i]*beta*sigma-B[i]^2*sigma^2
}
A <- cumsum(G.B)

GG[1:no.obs,]    <- c(A[matur],B[matur])/matur


X.0         <- c(1,theta)


for (i in 1:n)
                                                {
                                                    V[,1] <- 
c(0,rnorm(no.state-1,0,sqrt(diag(SigV)[2:no.state])))
                                                    W[,1] <- 
rnorm(no.obs,0,sqrt(diag(SigW)))

                                                    if (i==1) Xav[,i]   <- 
FF%*%X.0         +       V
                                                    if (i>1)  Xav[,i]   <- 
FF%*%Xav[,i-1]     +       V
                                                              Z[,i]     <- 
GG%*%Xav[,i]       +       W
                                                }

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Omega       <- matrix(0,nrow=2,ncol=2)

F <- array(dim=c(no.obs,no.obs,n))
Fminus <- matrix(ncol=no.obs,nrow=no.obs)

Ffast <- matrix(ncol=no.obs,nrow=no.obs)

P <- Pminus <- array(dim=c(no.state,no.state,n))
Pfast <- Pfastmin <- matrix(NA,ncol=no.state,nrow=no.state)

K <- array(dim=c(no.state,no.obs,n))
Kfast <- matrix(nrow=no.state,ncol=no.obs)
X <- Xhatminus <- array(dim=c(no.state,1,n))
Xfast <- Xfastmin <- matrix(ncol=1,nrow=no.state)


v <- array(dim=c(no.obs,1,n))
vfast <- matrix(nrow=no.obs,ncol=1)
Y <- matrix(ncol=1,nrow=no.obs)

log.lik <- numeric(n)

Rprof()

kalman <- function(a)
                                                    {

                                                    #(a <- startwerte)

                                                    theta.kal <- a[1]
                                                    kappa.kal <- a[2]
                                                    sigma.kal <- a[3]

                                                    x.kal     <- theta
                                                    Omega.kal <- sqrt(Q)
                                                    beta.kal  <- beta

                                                    Omega[2,2]      <- 
Omega.kal
                                                    Xhat            <- 
c(1,x.kal)
                                                    FF.kal[2,]      <- 
c(theta.kal-theta.kal*kappa.kal,kappa.kal)
                                                    SigV.kal[2,2]   <- 
sigma.kal^2
                                                    diag(SigW.kal)  <- 
a[4]^2

                                                    for (i in 
1:matur[length(matur)]){
                                                    B[i] <- 
(1-kappa.kal^(i-1))/(1-kappa.kal)
                                                    G.B[i] <- 
rho+B[i]*theta.kal*(1-kappa.kal)-0.5*beta.kal^2-B[i]*beta.kal*sigma.kal-B[i]^2*sigma.kal^2
                                                    }
                                                    A <- cumsum(G.B)

                                                    GG.kal[1:no.obs,]    <- 
c(A[matur],B[matur])/matur

                                                    for(i in 1:n){
                                                                  if 
(i==1){Xhat            <- c(1,x.kal); Omega[2,2] <- Q}else{Xhat  <- 
c(1,Xfast[2,1]);Omega[2,2] <- Pfast[2,2]}
                                                                  (Y[,1] <- 
Z[,i])
                                                                  (Xfastmin 
<- FF.kal%*%Xhat)
                                                                  (Pfastmin 
<- FF.kal%*%Omega%*%t(FF.kal)+SigV.kal )
                                                                  (Ffast <- 
GG.kal%*%Pfastmin%*%t(GG.kal)+SigW.kal)

                                                                  Fminus <- 
solve(Ffast,tol=1e-40)



                                                                  (Kfast <- 
Pfastmin%*%t(GG.kal)%*%Fminus)
                                                                  (vfast <- 
Y-GG.kal%*%Xfastmin)
                                                                  (Xfast <- 
Xfastmin+Kfast%*%vfast)
                                                                  (Pfast <- 
abs(Pfastmin-Kfast%*%GG.kal%*%Pfastmin))
                                                                  log.lik[i] 
<- -1/2*log(2*pi)-0.5*log(det(Ffast))-0.5*t(vfast)%*%Fminus%*%vfast
                                                                  }
                                                    sum(-log.lik)

                                                    }

kalman(c(theta,kappa,sqrt(Q),sqrt(R)))

startwerte <- c(theta,kappa,sqrt(Q),sqrt(R))
optim(startwerte*0.9,kalman,method="SANN")
c(theta,kappa,sqrt(Q),sqrt(R))
#nlmestimate$estimate
Rprof(NULL)
gamma2 <- summaryRprof()

_________________________________________________________________
Haben Spinnen Ohren? Finden Sie es heraus – mit dem MSN Suche Superquiz via



More information about the R-help mailing list