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 ______________________________________________ R-help@stat.math.ethz.ch 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.