\documentclass[11pt]{article}

\usepackage{../Math545}

\usepackage{amsmath,amssymb}

\pagestyle{plain}

\renewcommand{\labelenumi}{\arabic{enumi}.}
\renewcommand{\labelenumii}{\arabic{enumi}.\arabic{enumii}}
\renewcommand{\labelenumiii}{\arabic{enumi}.\arabic{enumii}.\arabic{enumiii}}

\setlength{\topmargin}{-0.75 in} \setlength{\textheight}{9.5 in}

\parindent0in

\newenvironment{packed_enum}{
\begin{enumerate}
  \setlength{\itemsep}{12pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
  \setlength{\topsep}{0pt}
}{\end{enumerate}}

\def\myGt{\mathbf{G}_t}
\def\myFt{\mathbf{F}_t}
\def\myHonet{\mathbf{H}_{1t}}
\def\myHtwot{\mathbf{H}_{2t}}
\def\myZ{\mathbf{Z}}
\def\evec{\mathbf{e}}
\def\bP{\mathbf{P}}
\def\bM{\mathbf{M}}
\def\bK{\mathbf{K}}
\def\bL{\mathbf{L}}
\def\bv{\mathbf{v}}
\def\ba{\mathbf{a}}
\def\bOmega{{\boldsymbol{\Omega}}}
\def\bSigma{{\boldsymbol{\Sigma}}}
\def\beps{{\boldsymbol{\epsilon}}}


\begin{document}

<<setup, cache=FALSE, echo=FALSE>>=
library(knitr)
# global chunk options
opts_chunk$set(cache=TRUE, autodep=TRUE)
options(scipen=999)
options(repos=c(CRAN="https://cloud.r-project.org/"))
@

\begin{center}
\textsclarge{MATH 545: Introduction to Time Series Analysis}

\smallskip
\textsc{State Space Models and The Kalman Filter}
\end{center}

Let $\by_{1:t}=\left\{ \by_{1},\dots ,\by_{t}\right\} $ be the observed data, and $\bX_{t}$ be the corresponding unobserved state vector at time $t$, for $t=1,\ldots,n$.  Inference goals for state space models are to

\begin{itemize}

\item \textbf{filtering (or tracking)}: estimate $\bX_t$ from all observations up to time $t$ via the distribution
\[
p(\bx_t | \by_{1:t})
\]


\item \textbf{prediction (or forecasting)}: estimate $\bX_t$ from observations up to a previous time point $t-s$ via the distribution
\[
p(\bx_t | \by_{1:(t-s)})
\]


\item \textbf{smoothing}: estimate $\bX_t$ from all past and possibly some future observations via the distribution
\[
p(\bx_t | \by_{1:(t+s)})
\]

\end{itemize}
as well as estimating system parameters. The quantities of interest are computed sequentially in the following recursive calculation:

\[
\begin{array}{llcl}
\textbf{Prediction} & p(\bx_t|\by_{1:(t-1)}) & = & \displaystyle{\int} p(\bx_t|\bx_{t-1}) p(\bx_{t-1} | \by_{1:(t-1)}) \; d\bx_{t-1} \\[9pt]
\textbf{Filtering} & p(\bx_t | \by_{1:t}) & = & \dfrac{p(\by_t|\bx_t) p(\bx_t | \by_{1:(t-1)})}{p(\by_t |\by_{1:(t-1)})}
\end{array}
\]


The above equations only provide an analytical solution if all densities in the state and observation equation are Gaussian, and both the state and the
observation equation are linear. The general form of a linear state space model is given by the following two equations
\begin{subequations}
\label{SSM}
\begin{eqnarray}
\bY_t &=&\myGt \bX_{t}+\bd_{t}+\myHonet \myZ_{t}\text{,\quad
}t=1,\dots
,n  \label{SSM1} \\
\bX_{t+1} &=&\myFt \bX_{t}+\evec_{t}+\myHtwot \myZ_{t}\text{.} \label{SSM2}
\end{eqnarray}
\end{subequations}
where $\bY_t$ is $(K \times 1)$, $\bX_t$ is $(L \times 1)$.  In (\ref{SSM})

\begin{itemize}

\item $\myGt$ is $\left( K\times L\right) $,

\item $\bd_{t}$ is $\left( K\times 1\right) $,

\item $\myHonet$ $\left( K\times \left( K+L\right) \right) $,

\item $\myFt $ is $\left( L\times L\right) $

\item $\evec_{t}$ is $\left( L\times 1\right) $

\item $\myHtwot$ is $\left( L\times \left( K + L\right) \right) $.

\end{itemize}
These matrices are referred to as the \textit{system matrices} and are assumed to be non-stochastic.  Process $\{\myZ_{t}\}$ is an $\left( \left( K+L\right) \times 1\right)$ vector of serially uncorrelated (independent), identically distributed disturbances with $\E [\myZ_{t}] =\boldsymbol{0}$ and $\Var[ \myZ_{t}] = \mathbf{I}_{K+L}$, the identity matrix.   We may define $\bW_t = \myHonet \myZ_{t}$ and $\bV_t = \myHtwot \myZ_{t}$ be the possibly dependent stochastic terms in the two equations; if $\myHonet \myHtwot^{\top }=\boldsymbol{0}$ for all $t$, then the measurement and transition equations are uncorrelated.  If the system matrices do not evolve with time, the state space model is called \textit{time-invariant} or \textit{time-homogeneous}. If the disturbances, $\myZ_{t}$, and initial state vector, $\bX_{1}$, are assumed to have a normal distribution, then the model is termed Gaussian.

\vspace{0.1 in} 

The formulation is completed by the assumption that the initial state vector $\bX_{1}$ is independent of $\myZ_{t}$ at all times $t$ and
has unconditional mean and variance $\ba_{1|0}$ and $\bP_{1|0}$ respectively.




\pagebreak

\textbf{The Kalman Filter:}  The Kalman filter is used for prediction, filtering and smoothing.
\begin{itemize}

\item \textit{Filtering} is concerned with computing $ \E [\bX
    _{t}|\by_{1:t}] $,

\smallskip

\item \textit{Prediction} is concerned with computing $\E [\bX_{t+1}|\by_{1:t}]$.

\smallskip

\item \textit{Smoothing} is concerned with estimating $\E[\bX_{t}|\by_{1:n}] $, for all $ t<n$.

\end{itemize}

The Kalman filter consists of a set of simple recursions that compute the latter quantities.  For simplicity, assume that
\[
\myHonet \myHtwot^{\top }= \boldsymbol{0} \qquad \qquad \bd_t = \evec_t = \boldsymbol{0}
\]
and
\[
\myHonet\myHonet^{\top }=\bSigma_{t} \qquad \qquad
\myHtwot\myHtwot^{\top}=\bOmega_{t}.
\]
The Kalman filter recursively computes the quantities

\begin{itemize}

\item $\ba_{t|t}=\E [\bX_{t}|\by_{1:t}] $,

\item $\ba_{t+1|t}=\E [\bX_{t+1}|\by_{1:t}] $,

\item $\bP_{t|t}=\text{MSE} [\bX_{t}|\by_{1:t}] $,

\item $ \bP_{t+1|t}=\text{MSE} [\bX_{t+1}|\by_{1:t}] $.

\end{itemize}
where $\text{MSE}$ is the \textit{mean-square error} or \textit{one-step ahead prediction variance}.  Then, starting with $\ba_{1|0}$ and $\bP_{1|0}$, $\ba_{t|t}$, $\ba_{t+1|t}$ and their MSEs are obtained by running for $t=1,\dots ,n$, the
recursions
\begin{subequations}
\label{prefil}
\begin{eqnarray}
\br_t &=&\by_{t}-\myGt \ba_{t|t-1}\text{,}\label{prefil0}\\[6pt]
\bM_{t}&=&\myGt \bP_{t|t-1}\myGt^{\top }+\bSigma_{t}  \label{prefil1} \\[6pt]
\ba_{t|t} &=&\ba_{t|t-1}+\bP_{t|t-1}\myGt^{\top } \bM_{t}^{-1}\br_t\text{,}
\label{prefil2} \\[6pt]
\bP_{t|t} &=&\bP_{t|t-1}-\bP_{t|t-1}\myGt^{\top} \bM_{t}^{-1}\myGt \bP_{t|t-1}\text{,}
\label{prefil3} \\[6pt]
\ba_{t+1|t} &=&\myFt  \ba_{t|t}\text{,}  \label{prefil4} \\[6pt]
\bP_{t+1|t} &=&\myFt  \bP_{t|t}\myFt ^{\top }+\bOmega_{t}\text{.}
\label{prefil5}
\end{eqnarray}
\end{subequations}


\textbf{Notes:}

\begin{itemize}

\item The quantities $\br_t$ and $\bM_{t}$ in (\ref{prefil1}) denote the one-step-ahead error in forecasting $\by_{t}$ conditional
    on the information set at time $t-1$ and its MSE respectively. This facilitates computation of the likelihood function as in the innovations algorithm.

\smallskip

\item Under Gaussianity, the quantities $\ba_{t|t}$ and $\ba_{t|t-1}$ are optimal (in MSE) estimators of $\bx_{t}$ conditional on the available information at after and before the $t$th data point is observed.

\smallskip

\item In general, $\ba_{t|t}$ and $\ba_{t|t-1}$ are optimal only within the class of \textbf{linear} estimators, that is, they are the \textit{minimum mean square linear estimators} of $\bx_{t}$ conditional on $\by_{1:t}$ and $\by_{1:(t-1)}$.

\smallskip

\item Finally, equations (\ref{prefil2}) and (\ref{prefil4}) can be combined together to yield recursions that only compute the one-step-ahead prediction estimates of $\bx_{t+1}$ given $\by_{1:t}$. Similarly, by taking together (\ref{prefil3}) and (\ref{prefil5}), a single set of recursions for the MSE is obtained, which goes directly from $\bP_{t|t-1}$ to $\bP_{t+1|t}$.

\vspace{0.1 in}

The resulting recursions, along with (\ref {prefil1}), are
\begin{subequations}
\label{predeq}
\begin{eqnarray}
\ba_{t+1|t} &=&\myFt  \ba_{t|t-1}+\bK_{t} \br_t \label{predeq1}\\[6pt]
\bK_{t}&=&\myFt  \bP_{t|t-1}\myGt^{\top } \bM_{t}^{-1} \label{predeq2}  \\[6pt]
\bP_{t+1|t} &=&\myFt  \bP_{t|t-1} \bL_{t}^{\top }+\bOmega_{t} \label{predeq3}\\[6pt]
\bL_{t}& = &\myFt -\bK_{t}\myGt\text{.}  \label{predeq4}
\end{eqnarray}
\end{subequations}

\end{itemize}

\pagebreak

\textbf{Parameter Estimation:} Another application of the Kalman filter is the estimation of any unknown parameters $\theta $ that appear in the system matrices. The likelihood for data $\mathbf{y}=\left( \by_{1},\ldots,\by_{n}\right)$ can be constructed as
\begin{equation*}
p\left( \by_{1},\dots ,\by_{n}\right) =p\left( \by_{n}| \by_{1:(n-1)} \right)
\cdots p\left( \by_{2}|\by_1\right) p\left( \by_{1}\right)
=\prod_{t=1}^{n} p\left( \by_{t}|\by_{1:(t-1)}\right) \text{.}
\end{equation*}
For a Gaussian model, by taking conditional expectations on both sides of the
observation equation, with $\bd_{t}\equiv 0$, then
\[
\E [\bY_{t}|\by_{1:(t-1)}] =\myGt \ba_{t|t-1} \qquad \qquad
\Var[ \bY_{t}|\by_{1:(t-1)}] =\bM_{t}.
\]
Crucially, the one-step-ahead prediction density $p\left(\by_{t}|\by_{1:(t-1)}\right) $ is the density of a multivariate normal random
variable with mean $\myGt \ba_{t|t-1}$ and covariance matrix $\bM_{t}$.  Thus, the log-likelihood function is given by
\begin{eqnarray*}
\ell = -\frac{n}{2}\log \left( 2\pi \right)
-\frac{1}{2}\sum_{t=1}^{n}\log \left( \det \bM_{t}\right) -
\frac{1}{2}\sum_{t=1}^{n} \br_t^{\top } \bM_{t}^{-1} \br_t\text{.}
\end{eqnarray*}
where
\[
\br_t = \left( \by_{t}-\myGt \ba_{t|t-1}\right)
\]
Typically, numerical procedures are used in order to maximize the
log-likelihood to obtain the ML estimates of the parameters $\theta$ which are consistent and asymptotically normal.

\bigskip

A \textit{smoothing} algorithm can be applied to a state space model; given a fixed set of data, estimates of the state vector are computed at each time $t$ in the sample period taking into account the full information set available.

\medskip
This algorithm computes $ \ba_{t|n}=\E [\bX_{t}|y_{1:n}] $ along with its MSE, $\bP_{t|n}$ , for all $t=1,\ldots ,n-1$ via a set of backward recursions and use the quantities $\ba_{t+1|t}$, $\ba_{t|t}$ and the MSE matrices $\bP_{t|t}$, $\bP_{t+1|t}$, which are obtained from (\ref{prefil}).

\medskip

In particular, to obtain $\ba_{t|n}$ and $\bP_{t|n}$, we start with $\ba_{n|n}$ and $\bP_{n|n}$ and we run backwards for $t=n-1,\dots ,0$
\begin{eqnarray*}
\ba_{t|n} &=&\ba_{t|t}+\bP_{t}^{\ast }\left( \ba_{t+1|n}-\ba_{t+1|t}\right) \\[6pt]
\bP_{t|n} &=&\bP_{t|t}+\bP_{t}^{\ast }\left( \bP_{t+1|n}-\bP_{t+1|t}\right)
\bP_{t}^{\ast }\\[6pt]
\bP_{t}^{\ast }& =& \bP_{t}\myFt ^{\top}\bP_{t+1|t} \text{.}
\end{eqnarray*}


If the state space model is not Gaussian, the likelihood can still be constructed in the same way using the minimum mean square linear estimators of the state vector.  However, the estimators $\widehat{ \theta }$ that maximize the likelihood are the \textit{quasi-maximum likelihood}
(QML) estimators of the parameters, which are also consistent and asymptotically normal.

\pagebreak

\textbf{EXAMPLE: AIR PASSENGERS DATA}   The data below represent monthly counts (in thousands) of international airline passengers, 1949 to 1960 ($n=144$).  There is evidence of a trend, and a twelve month periodicity. The usual model for these data $\{y_t\}$ is
\[
y_t = m_t + s_t + w_{t}
\]
where $m_t$ is the trend, and $s_t$ is the periodic function with period 12 (so that $s_{t+12k} = s_t$ for each $t$, and $k=1,2,\ldots$ such that $s_1 + s_2 + \cdots + s_{12} = 0$; note that the periodicity constraints imply that for $t=12,13,\ldots$.
\[
s_{t+1} = -s_t - s_{t-1} - \cdots - s_{t-10}
\]
Both components are presumed deterministic.  A generalization of the model allows for a \textit{stochastic level}, \textit{stochastic trend}, and \textit{stochastic seasonal} component. 


\begin{itemize}

\item A random level with random trend model has $\{M_t\}$ as a stochastic process
\[
M_{t} = M_{t-1} + B_{t-1} + V_{1(t-1)}
\]
where $\{B_t\}$ is another stochastic process with
\[
B_{t} = B_{t-1} + V_{2(t-1)}
\]
where $\{V_{1t}\}$ and $\{V_{2t}\}$ are independent white noise processes with variances $\sigma_1^2$ and $\sigma_2^2$.  

\item For the stochastic periodic component consider the series $\{S_t\}$ given by
\[
S_{t+1} = -S_t - S_{t-1} - \cdots - S_{t-10} + V_{3t}
\]
where $\{V_{3t}\}$ is another independent white noise process with variance $\sigma_3^2$.
\end{itemize}

\bigskip
Using the stochastic version of the model above
\[
Y_t = M_t + S_t + W_{t}
\]
and defining the $L=13$-dimensional state vector $\bX_t = \left[M_t \: \: B_{t} \: \: S_t \: \: S_{t-1} \: \: \cdots \: \: S_{t-10}\right]^\transpose$, a proposed state-space model for the air passenger data takes the form
\begin{eqnarray*}
Y_t & = & \bG \bX_t + W_{t} \\[6pt]
\bX_{t+1} & = & \bF \bX_t + \bV_{t}
\end{eqnarray*}
where $Y_t$ is scalar ($K=1$).  Here $\bG = \left[1 \: 0 \: 1 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0\right]^\transpose$ and
\[
\bF = \left[\begin{array}{rrrrrrrrrrrrr}
1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 &-1 &-1 &-1 &-1 &-1 &-1 &-1 &-1 &-1 &-1 &-1\\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0
\end{array}
\right]
\]
where $\{W_t\} \sim \WN(0,\sigma_W^2)$ and $\{\bV_t\} \sim \WN(0,\Omega_t)$ with
$\Omega_t = \text{diag}(\sigma_1^2,\sigma_2^2,\sigma_3^2,0,0,0,0,0,0,0,0,0,0)$.

\bigskip

Using the Kalman Filter, this model can be fitted using maximum (quasi) likelihood; the seventeen parameters of the model to be estimated are
\[
(\sigma_1^2,\sigma_2^2,\sigma_3^2,\sigma_W^2,M_1,B_1,S_1,S_2,S_3,\ldots,S_{11})
\]
After maximizing the likelihood, the states can be estimated.

<<KF1,comment='+', fig.width=7.25, fig.height=4.0, fig.align='center'>>=
library(MASS)
par(mar=c(4,4,1,0))
y<-as.numeric(AirPassengers);n<-length(y)
plot(y,pch=19,cex=0.6)
lines(y)
@

<<KF2,comment='+', fig.width=7.25, fig.height=4.1, fig.align='center'>>=
#K is the dimension of Y
#L is the dimension of X

K<-1;L<-13
G.t<-array(0,c(n,K,L));F.t<-array(0,c(n,L,L))
Sigma.t<-array(0,c(n,K,K));Omega.t<-array(0,c(n,L,L))

G1<-c(1,0,1,rep(0,10))
F1<-matrix(0,L,L)
F1[1,1:2]<-1;F1[2,2]<-1
F1[3,3:13]<--1;for(i in 3:12){F1[i+1,i]<-1}
@

<<KF3,comment='+', fig.width=7.25, fig.height=4.1, fig.align='center'>>=

##Compute Estimates

ml.filter.new<-function(x,Fm,Gm,nv,Lval,Kval,yv){

	Yt<-yv
	Sigmat<-exp(x[4])
	Omegat<-diag(c(exp(x[1:3]),rep(0,10)))

	att<-array(0,c(nv,Lval))
	apt<-array(0,c(nv+1,Lval))
	Ptt<-array(0,c(nv,Lval,Lval))
	Ppt<-array(0,c(nv+1,Lval,Lval))
	rt<-array(0,c(nv,Kval))
	Mt<-array(0,c(nv,Kval,Kval))

	apt[1,]<-x[5:length(x)]
	Ppt[1,,]<-Omegat

	#Filter Forward

	loglike<-0
	for(i in 1:nv){

		Gmat<-matrix(Gm,nrow=Kval,ncol=Lval)
		rt[i,]<-Yt[i] - Gmat %*% apt[i,]
		Mt[i,,]<-Gmat %*% Ppt[i,,] %*% t(Gmat) + Sigmat

		if(Sigmat == Inf || Mt[i,,] == 0 || is.na(Mt[i,,])){
			return(Inf)
		}
		Mtinv<-solve(Mt[i,,])

 		att[i,]<-apt[i,] + Ppt[i,,] %*% t(Gmat) %*% Mtinv %*% rt[i,]
		Ptt[i,,]<-Ppt[i,,]-Ppt[i,,] %*% t(Gmat) %*% Mtinv %*% Gmat %*% Ppt[i,,]
		apt[i+1,]<-Fm %*% att[i,]
		Ppt[i+1,,]<-Fm %*% Ptt[i,,] %*% t(Fm) + Omegat

		loglike<-loglike-0.5*log(det(matrix(Mt[i,,],Kval,Kval)))-
                         0.5*t(rt[i,]) %*% Mtinv %*% rt[i,]
	}

	return(as.numeric(loglike))

}
@

<<KF4,comment='+', fig.width=7.25, fig.height=4.1, fig.align='center'>>=
#Perform the optimization using optim
xv<-c(log(170),-10,log(11),log(0.01),rep(0,13))
AirPass.ml<-optim(xv,ml.filter.new,Fm=F1,Gm=G1,nv=n,Kval=K,Lval=L,yv=y,
            method="BFGS",control=list(trace=10,fnscale=-1,maxit=500,abstol=1e-6))
xv<-AirPass.ml$par
@

<<KF5,comment='+', fig.width=7.25, fig.height=4.1, fig.align='center'>>=
AirPass.ml

X.t<-matrix(0,nrow=length(AirPassengers),ncol=length(AirPass.ml$par)-4)
X.t[1,]<-AirPass.ml$par[5:length(AirPass.ml$par)]
sig.1<-exp(AirPass.ml$par[1])
sig.2<-exp(AirPass.ml$par[2])
sig.3<-exp(AirPass.ml$par[3])
sig.w<-exp(AirPass.ml$par[4])
c(sig.w,sig.1,sig.2,sig.3)

Y.t<-array(0,c(n,K))
for(i in 1:n){
	G.t[i,,]<-G1
	F.t[i,,]<-F1
	Sigma.t[i,,]<-sig.w
	Omega.t[i,,]<-diag(c(sig.1,sig.2,sig.3,rep(0,10)))
	Y.t[i,]<-y[i]
}
@

<<KF6,comment='+', fig.width=7.25, fig.height=4.1, fig.align='center'>>=
#Produce filtering and smoothing estimates

a.tt<-array(0,c(n,L))
a.pt<-array(0,c(n+1,L))
new.a.pt<-array(0,c(n+1,L))
P.tt<-array(0,c(n,L,L))
P.pt<-array(0,c(n+1,L,L))
new.P.pt<-array(0,c(n+1,L,L))
r.t<-array(0,c(n,K))
new.r.t<-array(0,c(n,K))
M.t<-array(0,c(n,K,K))
F.t<-array(0,c(n,L,L))
L.t<-array(0,c(n,L,K))
K.t<-array(0,c(n,L,L))

a.pt[1,]<-X.t[1,]
P.pt[1,,]<-diag(c(sig.1,sig.2,sig.3,rep(0,10)))
new.a.pt[1,]<-X.t[1,]
new.P.pt[1,,]<-diag(c(sig.1,sig.2,sig.3,rep(0,10)))

#Filter Forwards
log.like<-0
for(i in 1:n){

	Gmat<-matrix(G.t[i,,],nrow=K,ncol=L)
	r.t[i,]<-Y.t[i,] - Gmat %*% a.pt[i,]
	M.t[i,,]<-Gmat %*% P.pt[i,,] %*% t(Gmat) + Sigma.t[i,,]
	Mt.inv<-solve(M.t[i,,])

	F.t[i,,]<-F1

 	a.tt[i,]<-a.pt[i,] + P.pt[i,,] %*% t(Gmat) %*% Mt.inv %*% r.t[i,]
	P.tt[i,,]<-P.pt[i,,]-P.pt[i,,] %*% t(Gmat) %*% Mt.inv %*% Gmat %*% P.pt[i,,]
	a.pt[i+1,]<- F.t[i,,] %*% a.tt[i,]
	P.pt[i+1,,]<-F.t[i,,] %*% P.tt[i,,] %*% t(F.t[i,,]) + Omega.t[i,,]

	###
	new.r.t[i,]<-Y.t[i,] - Gmat %*% new.a.pt[i,]
	L.t[i,,]<-F.t[i,,] %*% new.P.pt[i,,] %*% t(Gmat) %*% solve(M.t[i,,])
	new.a.pt[i+1,]<- F.t[i,,] %*% new.a.pt[i,] + 
                          matrix(L.t[i,,],nrow=L,ncol=K) %*% new.r.t[i,]
	K.t[i,,]<-F.t[i,,] -  matrix(L.t[i,,],nrow=L,ncol=K) %*% Gmat
	new.P.pt[i+1,,]<-F.t[i,,] %*% new.P.pt[i,,] %*% t(K.t[i,,]) + Omega.t[i,,]

	log.like<-log.like-0.5*log(det(matrix(M.t[i,,],K,K)))-
                       0.5*t(r.t[i,]) %*% Mt.inv %*% r.t[i,]
}

#Smooth Backwards

Pstar.t<-array(0,c(n,L,L))
a.tT<-array(0,c(n,L))
P.tT<-array(0,c(n,L,L))
Pstar.t[n,,]<-P.tt[n,,]
a.tT[n,]<-a.tt[n,]
P.tT[n,,]<-P.tt[n,,]

for(i in (n-1):(L-3)){
	Pstar.t[i,,]<-P.tt[i,,] %*% t(F.t[i,,]) %*% P.pt[i+1,,]
	a.tT[i,]<-a.tt[i,]+Pstar.t[i,,] %*% (a.tT[i+1,] - a.pt[i+1,])
	P.tT[i,,]<-P.tt[i,,] + Pstar.t[i,,] %*% (P.tT[i+1,,]-P.pt[i+1,,]) %*% Pstar.t[i,,]
}

par(mar=c(4,4,1,0))

plot(a.pt[,1],type="l",ylim=range(-50,600),ylab="Passengers",xlab="Months",col="red",lwd=1)
lines(a.pt[,2],lty=2,col="blue",lwd=1)
lines(a.pt[,3],lty=3,col="black",lwd=1)
legend(0,600,c("Local Level","Long Term Slope","Seasonal"),
       lty=c(1,2,3),col=c("red","blue","black"),lwd=1)

plot(a.pt[,3],type="l",ylim=range(-100,100),ylab="Passengers",xlab="Months",col="black",lwd=1)
legend(0,100,c("Seasonal component"),col=c("black"),lwd=1)


Yfit<-as.vector(matrix(G1,nrow=1) %*% t(a.pt))
plot(Yfit,type="l",ylim=range(-50,600),ylab="Passengers",xlab="Months",col="red",lwd=1)
lines(y);points(y,pch=19,cex=0.6)
legend(0,600,c("Data","One-step Forecasts"),lty=c(1,1),col=c("black","red"),lwd=1,pch=c(19,NA))

Yresids<-Y.t-Yfit[-n]
plot(Yresids,type="l",ylim=range(-60,60),ylab="Residuals",xlab="Months",col="red",lwd=1)
lines(Yresids);points(Yresids,pch=19,cex=0.6)

acf(Yresids,lag.max=30,ylim=range(-1,1),main='ACF of Residuals')

pacf(Yresids,lag.max=30,ylim=range(-1,1),main='PACF of Residuals')

@
For the Air Passenger data the estimates are
\[
\sigma_1^2 = \Sexpr{sig.1} \qquad \sigma_2^2 = \Sexpr{sig.2} \qquad \sigma_3^2 = \Sexpr{sig.3} \qquad \sigma_W^2 = \Sexpr{sig.w}
\]
The initial seasonal component values are
\[
\Sexpr{round(xv[5:17],2)}
\]
\end{document}
plot(a.pt[,1],type="l",ylim=range(-50,600),ylab="Passengers",xlab="Months",lwd=1)
lines(a.pt[,2],lty=2,lwd=1)
lines(a.pt[,3],lty=4,lwd=1)
legend(0,600,c("Local Level","Long Term Slope","Seasonal"),lty=c(1,2,4),lwd=1)

Yfit<-as.vector(matrix(G1,nrow=1) %*% t(a.pt))  #Fitted values
plot(Yfit,type="l",ylim=range(-50,600),ylab="Passengers",xlab="Months",lwd=1,lty=4)
lines(y)
points(y,pch=19,cex=0.6)
legend(0,600,c("Data","One-step Forecasts"),lty=c(1,4),lwd=1,pch=c(19,NA))