| Time Series Analysis and Examples |
Example 10.3: Kalman Filtering: SSM Estimation With the EM Algorithm
The following example estimates the normal SSM of the mink-muskrat data by using the EM algorithm. The mink-muskrat series are detrended. Refer to Harvey (1989) for details of this data set. Since this EM algorithm uses filtering and smoothing, you can learn how to use the KALCVF and KALCVS calls to analyze the data. Consider the bivariate SSM:
is a 
identity matrix, the
observation noise has a time-invariant covariance matrix
, and the covariance matrix of the transition
equation is also assumed to be time invariant.
The initial state 
has mean 
and covariance 
.
For estimation, the matrix is fixed as
![[ 0.1 & 0.0 \ 0.0 & 0.1 ]](images/timeseriesexpls_timeseriesexplseq441.gif){kind=small}
is updated by the
smoothing procedure such that 
.
Note that this estimation requires an extra smoothing step
since the usual smoothing procedure does not produce 
.The EM algorithm maximizes the expected log-likelihood function given the current parameter estimates. In practice, the log-likelihood function of the normal SSM is evaluated while the parameters are updated by using the M-step of the EM maximization
![{f}^{i+1} & = & {s}_t(1)[{s}_{t-1}(0)]^{-1} \ {v}^{i+1} & = & \frac{1}t ( {s}_t... ...y}_t - {h}{z}_{t| t})^' + {h}{p}_{t| t} {h}^' ] \ \mu^{i+1} & = & {z}_{0| t}](images/timeseriesexpls_timeseriesexplseq444.gif)
represents the current iteration number, and
recursively such that
and the
initial value 
is derived by using the formula
![{p}_{t,t-1| t} = [ {i}- {p}_{t| t-1} {h}^' ({h}{p}_{t| t-1} {h}^' + {{r}}) {h} ] {f}{p}_{t-1| t-1}](images/timeseriesexpls_timeseriesexplseq450.gif){kind=small}
in the IML module LIK.
The log-likelihood function is written
.The iteration history is shown in Output 10.3.1.
As shown in Output 10.3.2, the eigenvalues
of 
are within the unit circle, which
indicates that the SSM is stationary.
However, the muskrat series (Y1) is
reported to be difference stationary.
The estimated parameters are almost
identical to those of the VAR(1) estimates.
Refer to Harvey (1989, p. 469).
Finally, multistep forecasts of 
are computed by using the KALCVF call.
Here is the code:
call kalcvf(pred,vpred,filt,vfilt,y,15,a,f,b,h,var,z0,vz0);
The predicted values of the state vector 
and
their standard errors are shown in Output 10.3.3.
Here is the code:
title 'SSM Estimation using EM Algorithm';
data one;
input y1 y2 @@;
datalines;
. . . data lines omitted . . .
;
proc iml;
start lik(y,pred,vpred,h,rt);
n = nrow(y);
nz = ncol(h);
et = y - pred*h`;
sum1 = 0;
sum2 = 0;
do i = 1 to n;
vpred_i = vpred[(i-1)*nz+1:i*nz,];
et_i = et[i,];
ft = h*vpred_i*h` + rt;
sum1 = sum1 + log(det(ft));
sum2 = sum2 + et_i*inv(ft)*et_i`;
end;
return(sum1+sum2);
finish;
start main;
use one;
read all into y var {y1 y2};
/*-- mean adjust series --*/
t = nrow(y);
ny = ncol(y);
nz = ny;
f = i(nz);
h = i(ny);
/*-- observation noise variance is diagonal --*/
rt = 1e-5#i(ny);
/*-- transition noise variance --*/
vt = .1#i(nz);
a = j(nz,1,0);
b = j(ny,1,0);
myu = j(nz,1,0);
sigma = .1#i(nz);
converge = 0;
logl0 = 0.0;
do iter = 1 to 100 while( converge = 0 );
/*--- construct big cov matrix --*/
var = ( vt || j(nz,ny,0) ) //
( j(ny,nz,0) || rt );
/*-- initial values are changed --*/
z0 = myu` * f`;
vz0 = f * sigma * f` + vt;
/*-- filtering to get one-step prediction and filtered value --*/
call kalcvf(pred,vpred,filt,vfilt,y,0,a,f,b,h,var,z0,vz0);
/*-- smoothing using one-step prediction values --*/
call kalcvs(sm,vsm,y,a,f,b,h,var,pred,vpred);
/*-- compute likelihood values --*/
logl = lik(y,pred,vpred,h,rt);
/*-- store old parameters and function values --*/
myu0 = myu;
f0 = f;
vt0 = vt;
rt0 = rt;
diflog = logl - logl0;
logl0 = logl;
itermat = itermat // ( iter || logl0 || shape(f0,1) || myu0` );
/*-- obtain P*(t) to get P_T_0 and Z_T_0 --*/
/*-- these values are not usually needed --*/
/*-- See Harvey (1989 p154) or Shumway (1988, p177) --*/
jt1 = sigma * f` * inv(vpred[1:nz,]);
p_t_0 = sigma + jt1*(vsm[1:nz,] - vpred[1:nz,])*jt1`;
z_t_0 = myu + jt1*(sm[1,]` - pred[1,]`);
p_t1_t = vpred[(t-1)*nz+1:t*nz,];
p_t1_t1 = vfilt[(t-2)*nz+1:(t-1)*nz,];
kt = p_t1_t*h`*inv(h*p_t1_t*h`+rt);
/*-- obtain P_T_TT1. See Shumway (1988, p180) --*/
p_t_ii1 = (i(nz)-kt*h)*f*p_t1_t1;
st0 = vsm[(t-1)*nz+1:t*nz,] + sm[t,]`*sm[t,];
st1 = p_t_ii1 + sm[t,]`*sm[t-1,];
st00 = p_t_0 + z_t_0 * z_t_0`;
cov = (y[t,]` - h*sm[t,]`) * (y[t,]` - h*sm[t,]`)` +
h*vsm[(t-1)*nz+1:t*nz,]*h`;
do i = t to 2 by -1;
p_i1_i1 = vfilt[(i-2)*nz+1:(i-1)*nz,];
p_i1_i = vpred[(i-1)*nz+1:i*nz,];
jt1 = p_i1_i1 * f` * inv(p_i1_i);
p_i1_i = vpred[(i-2)*nz+1:(i-1)*nz,];
if ( i > 2 ) then
p_i2_i2 = vfilt[(i-3)*nz+1:(i-2)*nz,];
else
p_i2_i2 = sigma;
jt2 = p_i2_i2 * f` * inv(p_i1_i);
p_t_i1i2 = p_i1_i1*jt2` + jt1*(p_t_ii1 - f*p_i1_i1)*jt2`;
p_t_ii1 = p_t_i1i2;
temp = vsm[(i-2)*nz+1:(i-1)*nz,];
sm1 = sm[i-1,]`;
st0 = st0 + ( temp + sm1 * sm1` );
if ( i > 2 ) then
st1 = st1 + ( p_t_ii1 + sm1 * sm[i-2,]);
else st1 = st1 + ( p_t_ii1 + sm1 * z_t_0`);
st00 = st00 + ( temp + sm1 * sm1` );
cov = cov + ( h * temp * h` +
(y[i-1,]` - h * sm1)*(y[i-1,]` - h * sm1)` );
end;
/*-- M-step: update the parameters --*/
myu = z_t_0;
f = st1 * inv(st00);
vt = (st0 - st1 * inv(st00) * st1`)/t;
rt = cov / t;
/*-- check convergence --*/
if ( max(abs((myu - myu0)/(myu0+1e-6))) < 1e-2 &
max(abs((f - f0)/(f0+1e-6))) < 1e-2 &
max(abs((vt - vt0)/(vt0+1e-6))) < 1e-2 &
max(abs((rt - rt0)/(rt0+1e-6))) < 1e-2 &
abs((diflog)/(logl0+1e-6)) < 1e-3 ) then
converge = 1;
end;
reset noname;
colnm = {'Iter' '-2*log L' 'F11' 'F12' 'F21' 'F22'
'MYU11' 'MYU22'};
print itermat[colname=colnm format=8.4];
eval = teigval(f0);
colnm = {'Real' 'Imag' 'MOD'};
eval = eval || sqrt((eval#eval)[,+]);
print eval[colname=colnm];
var = ( vt || j(nz,ny,0) ) //
( j(ny,nz,0) || rt );
/*-- initial values are changed --*/
z0 = myu` * f`;
vz0 = f * sigma * f` + vt;
free itermat;
/*-- multistep prediction --*/
call kalcvf(pred,vpred,filt,vfilt,y,15,a,f,b,h,var,z0,vz0);
do i = 1 to 15;
itermat = itermat // ( i || pred[t+i,] ||
sqrt(vecdiag(vpred[(t+i-1)*nz+1:(t+i)*nz,]))` );
end;
colnm = {'n-Step' 'Z1_T_n' 'Z2_T_n' 'SE_Z1' 'SE_Z2'};
print itermat[colname=colnm];
finish;
run;
Output 10.3.1: Iteration History
|
|
| Iter | -2*log L | F11 | F12 | F21 | F22 | MYU11 | MYU22 |
| 1.0000 | -154.010 | 1.0000 | 0.0000 | 0.0000 | 1.0000 | 0.0000 | 0.0000 |
| 2.0000 | -237.962 | 0.7952 | -0.6473 | 0.3263 | 0.5143 | 0.0530 | 0.0840 |
| 3.0000 | -238.083 | 0.7967 | -0.6514 | 0.3259 | 0.5142 | 0.1372 | 0.0977 |
| 4.0000 | -238.126 | 0.7966 | -0.6517 | 0.3259 | 0.5139 | 0.1853 | 0.1159 |
| 5.0000 | -238.143 | 0.7964 | -0.6519 | 0.3257 | 0.5138 | 0.2143 | 0.1304 |
| 6.0000 | -238.151 | 0.7963 | -0.6520 | 0.3255 | 0.5136 | 0.2324 | 0.1405 |
| 7.0000 | -238.153 | 0.7962 | -0.6520 | 0.3254 | 0.5135 | 0.2438 | 0.1473 |
| 8.0000 | -238.155 | 0.7962 | -0.6521 | 0.3253 | 0.5135 | 0.2511 | 0.1518 |
| 9.0000 | -238.155 | 0.7962 | -0.6521 | 0.3253 | 0.5134 | 0.2558 | 0.1546 |
| 10.0000 | -238.155 | 0.7961 | -0.6521 | 0.3253 | 0.5134 | 0.2588 | 0.1565 |
Output 10.3.2: Eigenvalues of F Matrix
|
Output 10.3.3: Multistep Prediction
Copyright © 2009 by SAS Institute Inc., Cary, NC, USA. All rights reserved.