That would likely be of very little help. RATS is written in C, so there is often very little resemblance to the matrix operations that would be used in (say) Matlab - calculations get rearranged for efficiency, specialized matrix routines are employed, etc.
RATS uses Koopman's algorithm, adapted to the dating scheme in our transition equation rather than his. This requires only a different placement on the t+0 post-data to t+1 pre-data calculation. (He does that at the end of our iteration; we do it at the beginning of the next). We're quite sure it's doing this correctly.
This was a procedure for doing the exact KF which we had before this was put directly into RATS. Note that the part of this beginning with if rank<>0 {... is excessively complicated. Because there's only one observable, the rank has to be either zero or one.
Code: Select all
proc KFExact start end
type integer start end
*
* Kalman filter conditional on initial observations
* Formulas are from "Exact Initial Kalman Filtering and Smoothing
* for Nonstationary Time Series Models" by Koopman, JASA, December 1997
*
* The basic state space model is as in the RATS instruction DLM, but with
* some of the components required to be fixed matrices rather than FRML's. Also,
* KFExact is restricted to a single measurement equation.
*
* X(t)=AX(t-1)+w(t)
*
* y(t)=c(t)'X(t)+v(t)
*
* Syntax:
*
* @KFExact( options ) start end
*
* Parameters:
* start end Range to estimate
*
* Options:
* A = A matrix (fixed nxn RECTANGULAR)
* C = FRML[RECT] giving the c matrices (should evaluate to a 1xn matrix)
* Y = FRML giving the y values
* SV = real number giving the (fixed) variance of v(t)
* SW = SYMMETRIC array giving the (fixed) variance of w(t)
* G = rxn matrix which transforms X to stationarity
*
* @KFExact sets %%LOGLIKE to the log likelihood value, and %RSS and %SEESQ to
* the sum of squared residuals and %RSS divided by the number of observations.
*
* Revision schedule:
* June, 2002, Written by Tom Doan, Estima
*
option RECT A
option FRML[RECT] c
option REAL sv 0.0
option SYMM sw
option FRML y
option RECT g ;* Transformation which gives stationarity
*
local vect nx hx
local symm nv hv nv1 hv1
local real norm resid sigmasq varterm
local integer n r rank time nobs
local rect p vc
local symm stv
local symm gw
local symm basehv
local symm d e tt ll nn etxx c22
local rect c12 revise
local vect cx
local integer trank
local integer i j
*
compute n = %rows(A)
if %defined(g)
compute r = %rows(G)
else
compute r = 0
*
dim nx(n) hx(n)
dim nv(n,n) hv(n,n)
dim d(n,n) e(n,n)
*
compute nx=%const(0.0)
*
* If r<>0, we have a stationary part. Figure out the finite and infinite components of
* the initial variance. The mean is assumed zero in either case.
*
if r>0 {
compute p =(g*a*tr(g))*inv(g*tr(g))
compute gw =g*sw*tr(g)
compute stv=%psdinit(p,gw)
compute nv=%mqform(stv,inv(g*tr(g))*g)
compute hv=(%identity(n)-tr(g)*inv(g*tr(g))*g)
compute basehv=hv
}
else {
compute nv=%const(0.0)
compute hv=%identity(n)
}
*
* This is the initial rank of the infinite part. This gets decremented as information
* from early observations resolves the uncertainty.
*
compute rank=n-r
*
* Initialize variables needed for computing the log likelihood function.
*
compute sigmasq=0.0,varterm=0.0,nobs=0
*
do time=start,end
*
compute cx =c(time)
*
* The finite part works the same as the usual KF.
*
compute nv1=a*nv*tr(a)+sw
compute nx =a*nx
*
* If rank is positive, we need to attack the infinite component. Once rank drops to
* zero (generally when we've seen the same number of data points as the initial rank,
* though it could be longer), we can just proceed as if it were a standard finite KF.
*
if rank<>0 {
compute hv1 =a*hv*tr(a)
*
* Compute a matrix T for which T HV T' is diagonal with 1's and 0's on the diagonal.
* %psddiag uses the initial matrix basehv to determine the natural scale of the
* diagonal values so machine zeros can be detected.
*
compute tt=%psddiag(hv1,basehv)
compute ll=tt*hv1*tr(tt)
compute nn=tt*nv1*tr(tt)
*
* Since ll is 0's and 1's, its rank is easy to determine.
*
compute trank=n
do i=1,n
if ll(i,i)<0.5 {
compute trank=i
break
}
end do i
*
dim c22(n-trank,n-trank) c12(trank,n-trank)
ewise c22(i,j)=nn(i+trank,j+trank)
ewise c12(i,j)=nn(i ,j+trank)
if (trank <> n) {
compute c22=inv(c22)
compute c12=c12*c22
compute etxx=%innerxx(c12)
}
*
* Finite component
*
ewise d(i,j)=%if(i<=trank.or.j<=trank,0.0,c22(i-trank,j-trank))
*
* h**-1 component
*
ewise e(i,j)=%if(i<=trank.and.j<=trank,(i==j),%if(i>trank.and.j>trank,etxx(i-trank,j-trank),-c12(i-trank,j)))
*
compute revise=nv*tr(a)*tr(tt)*d*t+hv*tr(a)*tr(tt)*e*t
*
compute norm=%qform(hv1,cx)
*
* If norm is (effectively) zero, the observation doesn't interact with the infinite
* components of the variance, so we drop through to the purely finite analysis.
*
if abs(norm)<1.e-12
compute hv=hv1
else {
if %valid(y(time)) {
compute nx=nx +(1.0/norm)*hv1*cx*(y(time)-%dot(cx,nx))
compute nv=nv1-(1.0/norm)*(hv1*cx*tr(cx)*nv1+nv1*cx*tr(cx)*hv1)+(1.0/norm**2)*(hv1*cx)*tr(cx)*hv1*(sv+%qform(nv1,cx))
compute hv=hv1-(1.0/norm)*(hv1*cx*tr(cx)*hv1)
compute rank =rank-1
compute varterm=varterm+log(norm)
compute nobs =nobs+1
}
else {
compute nv=nv1
compute hv=hv1
}
next
}
}
*
* Finite only. Just standard KF.
*
if %valid(y(time)) {
compute norm =sv+%qform(nv1,cx)
compute resid =y(time)-%dot(cx,nx)
compute sigmasq=sigmasq+resid**2/norm
compute varterm=varterm+log(norm)
compute nobs =nobs+1
compute vc =nv1*cx
compute nx =nx+(1.0/norm)*vc*resid
compute nv =nv1-(1.0/norm)*vc*tr(vc)
}
else
compute nv=nv1
end do time
*
* Figure out the log likelihood, etc.
*
compute %%loglike=-.5*nobs*log(sigmasq/nobs)-.5*varterm
compute %seesq =sigmasq/nobs
compute %rss =sigmasq
*
end