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 (i1) 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
