[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