Hi,
Following convention below:
y(t) = Ax(t)+Bu(t)+eps(t) # observation eq
x(t) = Cx(t-1)+Du(t)+eta(t) # state eq
I modified the following routine (which I copied from:
http://www.stat.pitt.edu/stoffer/tsa2/Rcode/Kall.R) to accommodate u(t), an
exogenous input to the system.
for (i in 2:N){
xp[[i]]=C%*%xf[[i-1]]
Pp[[i]]=C%*%Pf[[i-1]]%*%t(C)+Q
siginv=A[[i]]%*%Pp[[i]]%*%t(A[[i]])+R
sig[[i]]=(t(siginv)+siginv)/2 # make sure sig is symmetric
siginv=solve(sig[[i]]) # now siginv is sig[[i]]^{-1}
K=Pp[[i]]%*%t(A[[i]])%*%siginv
innov[[i]]=as.matrix(yobs[i,])-A[[i]]%*%xp[[i]]
xf[[i]]=xp[[i]]+K%*%innov[[i]]
Pf[[i]]=Pp[[i]]-K%*%A[[i]]%*%Pp[[i]]
like= like + log(det(sig[[i]])) + t(innov[[i]])%*%siginv%*%innov[[i]]
}
like=0.5*like
list(xp=xp,Pp=Pp,xf=xf,Pf=Pf,like=like,innov=innov,sig=sig,Kn=K)
}
I tried to fit my problem and observe that I got positive log likelihood mainly
because the log of determinant of my variance matrix is largely negative.
That's not good because they should be positive. Have anyone experience this
kind of instability?
Also, I realize that I have about 800 sample points. The above routine when
being plugged to optim becomes very slow. Could anyone share a faster way to
compute kalman filter?
And my last problem is, optim with my defined feasible space does not converge.
I have about 20 variables that I need to identify using MLE method. Is there any
other way that I can try out? I tried most of the methods available in optim
already. They do not converge at all...... Thank you.
- adschai
On Thu, 15 Nov 2007, adschai at optonline.net wrote:> Hi, > > Following convention below: > y(t) = Ax(t)+Bu(t)+eps(t) # observation eq > x(t) = Cx(t-1)+Du(t)+eta(t) # state eq > > I modified the following routine (which I copied from: http://www.stat.pitt.edu/stoffer/tsa2/Rcode/Kall.R) to accommodate u(t), an exogenous input to the system. > > for (i in 2:N){ > xp[[i]]=C%*%xf[[i-1]] > Pp[[i]]=C%*%Pf[[i-1]]%*%t(C)+Q > siginv=A[[i]]%*%Pp[[i]]%*%t(A[[i]])+R > sig[[i]]=(t(siginv)+siginv)/2 # make sure sig is symmetric > siginv=solve(sig[[i]]) # now siginv is sig[[i]]^{-1} > K=Pp[[i]]%*%t(A[[i]])%*%siginv > innov[[i]]=as.matrix(yobs[i,])-A[[i]]%*%xp[[i]] > xf[[i]]=xp[[i]]+K%*%innov[[i]] > Pf[[i]]=Pp[[i]]-K%*%A[[i]]%*%Pp[[i]] > like= like + log(det(sig[[i]])) + t(innov[[i]])%*%siginv%*%innov[[i]] > } > like=0.5*like > list(xp=xp,Pp=Pp,xf=xf,Pf=Pf,like=like,innov=innov,sig=sig,Kn=K) > } > > I tried to fit my problem and observe that I got positive log likelihood > mainly because the log of determinant of my variance matrix is largely > negative. That's not good because they should be positive. Have anyone > experience this kind of instability?Why are you expecting that? The magnitude of the log-likelihood depends on the scale of your data, and hence so does the sign of its log.> Also, I realize that I have about 800 sample points. The above routine > when being plugged to optim becomes very slow. Could anyone share a > faster way to compute kalman filter?Try> help.search("kalman", agrep=FALSE)kalsmo.car(cts) Compute Components with the Kalman Smoother kalsmo.comp(cts) Estimate Componenents with the Kalman Smoother nbkal(repeated) Negative Binomial Models with Kalman Update extended(sspir) Iterated Extended Kalman Smoothing kfilter(sspir) Kalman filter for Gaussian state space model kfs(sspir) (Iterated extended) Kalman smoother smoother(sspir) Kalman smoother for Gaussian state space model tsmooth(timsac) Kalman Filter KalmanLike(stats) Kalman Filtering> And my last problem is, optim with my defined feasible space does not > converge. I have about 20 variables that I need to identify using MLE > method. Is there any other way that I can try out? I tried most of the > methods available in optim already. They do not converge at all......Most likely because of errors in your code. Optim solves quite large Kalman filter problems for arima(). -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
Kalman filter for general state space models, especially in its naive version, is known for its numerical instability. This is the reason why people developed square root filters, based on Cholesky decomposition of variance matrices. In package dlm the implementation of Kalman filter is based on the singular value decomposition (SVD) of the relevant variance matrices, providing for a more robust algorithm than the standard square root filter. The lack of convergence of optimization algorithms in finding MLEs of unknown parameters is not surprising to me, especially if you have many parameters. When using state space models it is easy to end up with a fairly flat, or multimodal likelihood function. These are two signs of poor identifiability of the model. You can live with it, but be aware that it is there. My suggestion is to start the optimization from several different initial values and compare maximized values of the likelihood. Simulated annealing may be used to better explore the parameter space. HTH, Giovanni> Date: Thu, 15 Nov 2007 04:41:26 +0000 (GMT) > From: adschai at optonline.net > Sender: r-help-bounces at r-project.org > Priority: normal > Precedence: list > > Hi, > > Following convention below: > y(t) = Ax(t)+Bu(t)+eps(t) # observation eq > x(t) = Cx(t-1)+Du(t)+eta(t) # state eq > > I modified the following routine (which I copied from: http://www.stat.pitt.edu/stoffer/tsa2/Rcode/Kall.R) to accommodate u(t), an exogenous input to the system. > > for (i in 2:N){ > xp[[i]]=C%*%xf[[i-1]] > Pp[[i]]=C%*%Pf[[i-1]]%*%t(C)+Q > siginv=A[[i]]%*%Pp[[i]]%*%t(A[[i]])+R > sig[[i]]=(t(siginv)+siginv)/2 # make sure sig is symmetric > siginv=solve(sig[[i]]) # now siginv is sig[[i]]^{-1} > K=Pp[[i]]%*%t(A[[i]])%*%siginv > innov[[i]]=as.matrix(yobs[i,])-A[[i]]%*%xp[[i]] > xf[[i]]=xp[[i]]+K%*%innov[[i]] > Pf[[i]]=Pp[[i]]-K%*%A[[i]]%*%Pp[[i]] > like= like + log(det(sig[[i]])) + t(innov[[i]])%*%siginv%*%innov[[i]] > } > like=0.5*like > list(xp=xp,Pp=Pp,xf=xf,Pf=Pf,like=like,innov=innov,sig=sig,Kn=K) > } > > I tried to fit my problem and observe that I got positive log likelihood mainly because the log of determinant of my variance matrix is largely negative. That's not good because they should be positive. Have anyone experience this kind of instability? > > Also, I realize that I have about 800 sample points. The above routine when being plugged to optim becomes very slow. Could anyone share a faster way to compute kalman filter? > > And my last problem is, optim with my defined feasible space does not converge. I have about 20 variables that I need to identify using MLE method. Is there any other way that I can try out? I tried most of the methods available in optim already. They do not converge at all...... Thank you. > > - adschai > > ______________________________________________ > R-help at r-project.org 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. > >-- Giovanni Petris <GPetris at uark.edu> Department of Mathematical Sciences University of Arkansas - Fayetteville, AR 72701 Ph: (479) 575-6324, 575-8630 (fax) http://definetti.uark.edu/~gpetris/
adschai at optonline.net wrote:> Hi, > > Following convention below: > y(t) = Ax(t)+Bu(t)+eps(t) # observation eq > x(t) = Cx(t-1)+Du(t)+eta(t) # state eq > > I modified the following routine (which I copied from: http://www.stat.pitt.edu/stoffer/tsa2/Rcode/Kall.R) to accommodate u(t), an exogenous input to the system. > > for (i in 2:N){ > xp[[i]]=C%*%xf[[i-1]] > Pp[[i]]=C%*%Pf[[i-1]]%*%t(C)+Q > siginv=A[[i]]%*%Pp[[i]]%*%t(A[[i]])+R > sig[[i]]=(t(siginv)+siginv)/2 # make sure sig is symmetric > siginv=solve(sig[[i]]) # now siginv is sig[[i]]^{-1} > K=Pp[[i]]%*%t(A[[i]])%*%siginv > innov[[i]]=as.matrix(yobs[i,])-A[[i]]%*%xp[[i]] > xf[[i]]=xp[[i]]+K%*%innov[[i]] > Pf[[i]]=Pp[[i]]-K%*%A[[i]]%*%Pp[[i]] > like= like + log(det(sig[[i]])) + t(innov[[i]])%*%siginv%*%innov[[i]] > } > like=0.5*like > list(xp=xp,Pp=Pp,xf=xf,Pf=Pf,like=like,innov=innov,sig=sig,Kn=K) > } > > I tried to fit my problem and observe that I got positive log likelihood mainly because the log of determinant of my variance matrix is largely negative. That's not good because they should be positive. Have anyone experience this kind of instability? >The determinant should not be negative (the line "# make sure sig is symmetric" should look after that) but the log can be negative.> Also, I realize that I have about 800 sample points. The above routine when being plugged to optim becomes very slow. Could anyone share a faster way to compute kalman filter?The KF recursion you show is not time varying, but the code you show looks like it may allow for time varying models. (I'm just guessing, I'm not familiar with the code.) If you do not need the time varying aspect then this is probably a slow way to implement. Package dse1 in the dse bundle implements the recursion the way you show, with an exogenous input u(t), so you don't even have to modify the code. It is implemented in R for demonstration, but also in fortran for speed. See the dse Users Guide for more details, and also ?SS and ?l.SS. One caveat is that the way the likelihood is cumulated over i in your code above allows for time varying sig, which in theory can be important for a diffuse prior. In dse the likelihood is calculated using the residual to get a sample estimate of sig, which is faster. I have not found cases where this makes much difference (but would be interested to know of one).> And my last problem is, optim with my defined feasible space does not converge. I have about 20 variables that I need to identify using MLE method. Is there any other way that I can try out? I tried most of the methods available in optim already. They do not converge at all...... Thank you.This used to be a well known problem for multivariate state space models, but seems to be largely forgotten. It does not seriously affect univariate models. For a review of the old literature see section 5 of Gilbert 1993, available at <http://www.bank-banque-canada.ca/pgilbert/research.php#MiscResearch>. The bottom line is that you are in trouble trying to use hill climbing methods and state-space models unless you have a very good starting point, or are luck. This is a feature of the parameter space, not a reflection on the optimization routine. One solution is to start by estimating a VAR or ARMA model and then convert it to a state space model, which was one of the original purposes of dse. (See ?bft for example.) Paul Gilbert> - adschai > > ______________________________________________ > R-help at r-project.org 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.=================================================================================== La version fran?aise suit le texte anglais. ------------------------------------------------------------------------------------ This email may contain privileged and/or confidential information, and the Bank of Canada does not waive any related rights. Any distribution, use, or copying of this email or the information it contains by other than the intended recipient is unauthorized. If you received this email in error please delete it immediately from your system and notify the sender promptly by email that you have done so. ------------------------------------------------------------------------------------ Le pr?sent courriel peut contenir de l'information privil?gi?e ou confidentielle. La Banque du Canada ne renonce pas aux droits qui s'y rapportent. Toute diffusion, utilisation ou copie de ce courriel ou des renseignements qu'il contient par une personne autre que le ou les destinataires d?sign?s est interdite. Si vous recevez ce courriel par erreur, veuillez le supprimer imm?diatement et envoyer sans d?lai ? l'exp?diteur un message ?lectronique pour l'aviser que vous avez ?limin? de votre ordinateur toute copie du courriel re?u.