\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\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}}}

\def\ctitle{MATH 545: Introduction to Time Series Analysis}
\def\ttitle{Winter 2025}
\def\stitle{State Space Models And The Kalman Filter}


\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/"))
@

\maketitle

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: for $t=1,\dots ,n $
\begin{subequations}
\label{SSM}
\begin{eqnarray}
\bY_t &=&\bF_t \bX_{t}+\bd_{t}+\myHonet \bZ_t\text{,\quad
} \label{SSM1} \\
\bX_{t+1} &=&\bG_t \bX_{t}+\evec_{t}+\myHtwot \bZ_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 $\bF_t$ 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 $\bG_t $ 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 assumed to be non-stochastic.  Process $\{\bZ_t\}$ is a  an $\left( \left( K+L\right) \times 1\right)$ process with
\[
\E [\bZ_t] =\boldsymbol{0} \qquad \Var[ \bZ_t] = \mathbf{I}_{K+L}.
\]  
We may write
\[
\bV_t = \myHonet \bZ_t \qquad \bW_t = \myHtwot \bZ_t
\]
as the possibly dependent stochastic error terms in the two equations: if $\myHonet \myHtwot^{\top }=\boldsymbol{0}$ for all $t$, then the error terms 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}. Often, $\{\bZ_t\}$, and initial state vector, $\bX_{1}$, are assumed to be Gaussian random variables.

\vspace{0.1 in} 

The formulation is completed by the assumption that the initial state vector $\bX_{1}$ is independent of $\bZ_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 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}-\bF_t \ba_{t|t-1}\text{,}\label{prefil0}\\[3pt]
\bM_{t}&=&\bF_t \bP_{t|t-1}\bF_t^{\top }+\bSigma_{t}  \label{prefil1} \\[3pt]
\ba_{t|t} &=&\ba_{t|t-1}+\bP_{t|t-1}\bF_t^{\top } \bM_{t}^{-1}\br_t\text{,}
\label{prefil2} \\[6pt]
\bP_{t|t} &=&\bP_{t|t-1}-\bP_{t|t-1}\bF_t^{\top} \bM_{t}^{-1}\bF_t \bP_{t|t-1}\text{,}
\label{prefil3} \\[3pt]
\ba_{t+1|t} &=&\bG_t  \ba_{t|t}\text{,}  \label{prefil4} \\[3pt]
\bP_{t+1|t} &=&\bG_t  \bP_{t|t}\bG_t ^{\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; 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} &=&\bG_t  \ba_{t|t-1}+\bK_{t} \br_t \label{predeq1}\\[3pt]
\bK_{t}&=&\bG_t  \bP_{t|t-1}\bF_t^{\top } \bM_{t}^{-1} \label{predeq2}  \\[3pt]
\bP_{t+1|t} &=&\bG_t  \bP_{t|t-1} \bL_{t}^{\top }+\bOmega_{t} \label{predeq3}\\[3pt]
\bL_{t}& = &\bG_t -\bK_{t}\bF_t\text{.}  \label{predeq4}
\end{eqnarray}
\end{subequations}

\end{itemize}

\pagebreak

\textbf{Parameter Estimation:} 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)}] =\bF_t \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 $\bF_t \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}-\bF_t \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

\textbf{Smoothing:} 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.

\begin{itemize}

\item 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

\item  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}\bG_t ^{\top}\bP_{t+1|t}^{-1} \text{.}
\end{eqnarray*}

\end{itemize}

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 + v_{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} + W_{1,t-1}
\]
where $\{B_t\}$ is another stochastic process with
\[
B_{t} = B_{t-1} + W_{2,t-1}
\]
where $\{W_{1t}\}$ and $\{W_{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} + W_{3t}
\]
where $\{W_{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 + V_{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 & = & \bF \bX_t + V_{t} \\[6pt]
\bX_{t+1} & = & \bG \bX_t + \bW_{t}
\end{eqnarray*}
where $Y_t$ is scalar ($K=1$).  Here $\bF = \left[1 \: 0 \: 1 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0\right]^\transpose$ and
\[
\bG = \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 $\{V_t\} \sim \WN(0,\sigma_V^2)$ and $\{\bW_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_V^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
F.t<-array(0,c(n,K,L));G.t<-array(0,c(n,L,L))
Sigma.t<-array(0,c(n,K,K));Omega.t<-array(0,c(n,L,L))

F1<-c(1,0,1,rep(0,10))
G1<-matrix(0,L,L)
G1[1,1:2]<-1;G1[2,2]<-1
G1[3,3:13]<--1;for(i in 3:12){G1[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){

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

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

 		att[i,]<-apt[i,] + Ppt[i,,] %*% tFmat %*% Mtinv %*% rt[i,]
		Ptt[i,,]<-Ppt[i,,]-Ppt[i,,] %*% tFmat %*% Mtinv %*% Fmat %*% Ppt[i,,]
		apt[i+1,]<-Gm %*% att[i,]
		Ppt[i+1,,]<-Gm %*% Ptt[i,,] %*% t(Gm) + 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.v<-exp(AirPass.ml$par[4])
c(sig.v,sig.1,sig.2,sig.3)

Y.t<-array(0,c(n,K))
for(i in 1:n){
	F.t[i,,]<-F1
	G.t[i,,]<-G1
	Sigma.t[i,,]<-sig.v
	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.5, 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))
G.t<-array(0,c(n,L,L))
K.t<-array(0,c(n,L,K))
L.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.0
for(i in 1:n){

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

	G.t[i,,]<-G1

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

	###  Second recursions
	new.r.t[i,]<-Y.t[i,] - Fmat %*% new.a.pt[i,]
	K.t[i,,]<-G.t[i,,] %*% new.P.pt[i,,] %*% tFmat %*% solve(M.t[i,,])
	new.a.pt[i+1,]<- G.t[i,,] %*% new.a.pt[i,] + 
                          matrix(K.t[i,,],nrow=L,ncol=K) %*% new.r.t[i,]
	L.t[i,,]<-G.t[i,,] -  matrix(K.t[i,,],nrow=L,ncol=K) %*% Fmat
	new.P.pt[i+1,,]<-G.t[i,,] %*% new.P.pt[i,,] %*% t(L.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):1){
	Pstar.t[i,,]<-P.tt[i,,] %*% t(G.t[i,,]) %*% solve(P.pt[i+1,,]+diag(1e-10,L))
	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,,]
}
#a.tT

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

plot(a.pt[,1],type="l",ylim=range(-50,650),
    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)
@

<<KF6a,comment='+', fig.width=7.25, fig.height=5.5, fig.align='center'>>=
par(mar=c(4,4,2,2))
Ypred<-as.vector(matrix(F1,nrow=1) %*% t(a.pt))
Ysmooth<-as.vector(matrix(F1,nrow=1) %*% t(a.tT))
plot(y,type="p",ylim=range(-50,650),pch=19,cex=0.4,
        ylab="Passengers",xlab="Months",lwd=1)
lines(Ypred,col='red')
lines(Ysmooth,col='blue')
legend(0,600,c("Data","Predicted", "Smoothed"),lty=c(NA,1,1),
        col=c("black","red","blue"),lwd=1,pch=c(19,NA,NA))
@
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_V^2 = \Sexpr{sig.v}
\]
The initial seasonal component values are
\[
\Sexpr{round(xv[5:17],2)}
\]
We note that this model does explain the observed variation in the data perfectly, as the variance of the prediction residuals does not seem to be constant.  Also, the observation noise $\{V_t\}$ has very small variance, $\sigma_V^2 = \Sexpr{sig.v}$, implying that the latent structure is explaining the variance almost perfectly.

\pagebreak
Residuals: There are different ways that we can form residuals, for example, we can inspect the prediction residuals
\[
r_t = y_t - \bF \ba_{t|t-1}
\]
that measure discrepancies in one-step-ahead forecasts, and are used in the likelihood construction, or the filtering residuals
\[
r_t = y_t - \bF \ba_{t|t}
\]
to assess the fit in-sample, or the smoothing residuals
\[
r_t = y_t - \bF \ba_{t|n}
\]
to assess the fit in light of the complete data set.

<<KF6b,comment='+', fig.width=7.25, fig.height=5.0, fig.align='center'>>=
Yresids<-Y.t-Ypred[-n]
par(mar=c(4,4,2,2))
plot(Yresids,type="l",ylim=range(-60,60),
    ylab="Prediction Residuals",xlab="Months",col="red",lwd=1)
lines(Yresids);points(Yresids,pch=19,cex=0.6)
@

<<KF6c,comment='+', fig.width=7.25, fig.height=5.0, fig.align='center',echo=FALSE>>=
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')
@

\pagebreak

\textbf{Example: Wine data}  These data comprise sales of red wine in Australia, from January 1980 to July 1995, measured in thousands of litres.  The data are plotted on the log scale below:

<<wine0,fig.height=5.0,fig.width=7>>=
library(Rssa)
data(AustralianWine)
Wine<-as.numeric(AustralianWine[,4])
n<-length(Wine)
tvec<-c(1:n)
lWine<-log(Wine)
par(mar=c(4,4,2,2))
plot(ts(log(AustralianWine[,4]),start=c(1980,1),freq=12),
        xlab='Year',ylab='log(Wine)',ylim=range(6,8.5))
ind <- 1980 + c(0:186)/12
points(ind, lWine, pch = 19, cex=0.5)
y<-lWine
@

We fit the same model
\[
Y_t = M_t + S_t + V_{t}
\]
with
\[
M_{t} = M_{t-1} + B_{t-1} + W_{1,t-1}
\qquad \qquad
B_{t} = B_{t-1} + W_{2,t-1}
\]
yielding state vector $\bX_t = \left[M_t \: \: B_{t} \: \: S_t \: \: S_{t-1} \: \: \cdots \: \: S_{t-10}\right]^\transpose$ and
\begin{eqnarray*}
Y_t & = & \bF \bX_t + V_{t} \\[6pt]
\bX_{t+1} & = & \bG \bX_t + \bW_{t}
\end{eqnarray*}

<<wine2,comment='+', fig.width=7.25, fig.height=4.1, fig.align='center', echo=FALSE>>=
K<-1;L<-13
F.t<-array(0,c(n,K,L));G.t<-array(0,c(n,L,L))
Sigma.t<-array(0,c(n,K,K));Omega.t<-array(0,c(n,L,L))

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


<<wine3,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center'>>=
#Perform the optimization using optim

xv<-c(0,0,0,0,6,0,rnorm(11,0,0.1))
LogWine.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=-500,maxit=500,abstol=1e-6))
xv<-LogWine.ml$par
LogWine.ml

X.t<-matrix(0,nrow=length(y),ncol=length(LogWine.ml$par)-4)
X.t[1,]<-LogWine.ml$par[5:length(LogWine.ml$par)]
sig.1<-exp(LogWine.ml$par[1])
sig.2<-exp(LogWine.ml$par[2])
sig.3<-exp(LogWine.ml$par[3])
sig.v<-exp(LogWine.ml$par[4])
@
For the Wine 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_V^2 = \Sexpr{sig.v}
\]
with seasonal component estimates $\Sexpr{round(xv[5:17],2)}$.

<<wine5,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center', echo=FALSE>>=
Y.t<-array(0,c(n,K))
for(i in 1:n){
	F.t[i,,]<-F1
	G.t[i,,]<-G1
	Sigma.t[i,,]<-sig.v
	Omega.t[i,,]<-diag(c(sig.1,sig.2,sig.3,rep(0,10)))
	Y.t[i,]<-y[i]
}
#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))
G.t<-array(0,c(n,L,L))
K.t<-array(0,c(n,L,K))
L.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.0
for(i in 1:n){

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

	G.t[i,,]<-G1

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

	###  Second recursions
	new.r.t[i,]<-Y.t[i,] - Fmat %*% new.a.pt[i,]
	K.t[i,,]<-G.t[i,,] %*% new.P.pt[i,,] %*% tFmat %*% solve(M.t[i,,])
	new.a.pt[i+1,]<- G.t[i,,] %*% new.a.pt[i,] + 
                          matrix(K.t[i,,],nrow=L,ncol=K) %*% new.r.t[i,]
	L.t[i,,]<-G.t[i,,] -  matrix(K.t[i,,],nrow=L,ncol=K) %*% Fmat
	new.P.pt[i+1,,]<-G.t[i,,] %*% new.P.pt[i,,] %*% t(L.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):1){
	Pstar.t[i,,]<-P.tt[i,,] %*% t(G.t[i,,]) %*% solve(P.pt[i+1,,]+diag(1e-10,L))
	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,,]
}
#a.tT

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

plot(a.pt[,1],type="l",ylim=range(-2,9),
    ylab="Sales",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,5,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(-5,5),
    ylab="Sales",xlab="Months",col="black",lwd=1)
legend(0,5,c("Seasonal component"),col=c("black"),lwd=1)
@

<<wine6a,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center', echo=FALSE>>=
par(mar=c(4,4,2,2))
Ypred<-as.vector(matrix(F1,nrow=1) %*% t(a.pt))
Ysmooth<-as.vector(matrix(F1,nrow=1) %*% t(a.tT))
plot(y,type="p",ylim=range(4.75,9),pch=19,cex=0.4,
        ylab="Sales",xlab="Months",lwd=1)
lines(Ypred,col='red')
lines(Ysmooth,col='blue')
legend(0,9,c("Data","Predicted", "Smoothed"),lty=c(NA,1,1),
        col=c("black","red","blue"),lwd=1,pch=c(19,NA,NA))
@

<<wine6b,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center',echo=FALSE>>=
Yresids<-Y.t-Ypred[-n]
par(mar=c(4,4,2,2))
plot(Yresids,type="l",ylim=range(-2,2),
    ylab="Prediction Residuals",xlab="Months",col="red",lwd=1)
lines(Yresids);points(Yresids,pch=19,cex=0.6)
@

<<wine6c,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center',echo=FALSE>>=
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')
@

\pagebreak
For comparison, we now fit the simpler local-level only model
\[
Y_t = M_t + S_t + V_{t}
\]
with
\[
M_{t} = M_{t-1} + W_{1,t-1}
\]
yielding state vector $\bX_t = \left[M_t \: \: S_t \: \: S_{t-1} \: \: \cdots \: \: S_{t-10}\right]^\transpose$ and
\begin{eqnarray*}
Y_t & = & \bF \bX_t + V_{t} \\[6pt]
\bX_{t+1} & = & \bG \bX_t + \bW_{t}
\end{eqnarray*}
with $\bF = \left[1 \: 1 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0 \: 0\right]^\transpose$ and
\[
\bG = \left[\begin{array}{rrrrrrrrrrrrr}
1 &  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & -1 &-1 &-1 &-1 &-1 &-1 &-1 &-1 &-1 &-1 &-1\\
0 &  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 & 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 & 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 & 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 & 1 & 0 & 0\\
0 &  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0
\end{array}
\right]
\]
where $\{V_t\} \sim \WN(0,\sigma_V^2)$ and $\{\bW_t\} \sim \WN(0,\Omega_t)$ with
$\Omega_t = \text{diag}(\sigma_1^2,\sigma_2^2,0,0,0,0,0,0,0,0,0,0)$.

<<wine2-2,comment='+', fig.width=7.25, fig.height=4.1, fig.align='center', echo=FALSE>>=
K<-1;L<-12
F.t<-array(0,c(n,K,L));G.t<-array(0,c(n,L,L))
Sigma.t<-array(0,c(n,K,K));Omega.t<-array(0,c(n,L,L))

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


<<wine2-3,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center',echo=FALSE>>=
#Perform the optimization using optim
ml.filter.new<-function(x,Fm,Gm,nv,Lval,Kval,yv){

	Yt<-yv
	Sigmat<-exp(x[3])
	Omegat<-diag(c(exp(x[1:2]),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[4:length(x)]
	Ppt[1,,]<-Omegat

	#Filter Forward

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

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

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

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

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

	return(as.numeric(loglike))

}
@


<<wine2-4,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center'>>=
##Compute Estimates
d<-12
ncyc<-round(n/12)+1
wvec<-rep(0,12)
for(j in 1:12){
    ind<-c(0:ncyc)*12+j
    ind<-ind[ind<=n]
    wvec[j]<-mean(y[ind])
}
svec<-wvec-wvec[1]
xv<-c(0,0,0,6,rnorm(11,0,0.2))
LogWine.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=-500,maxit=500,abstol=1e-6))
xv<-LogWine.ml$par
LogWine.ml

X.t<-matrix(0,nrow=length(y),ncol=length(LogWine.ml$par)-3)
X.t[1,]<-LogWine.ml$par[4:length(LogWine.ml$par)]
sig.1<-exp(LogWine.ml$par[1])
sig.2<-exp(LogWine.ml$par[2])
sig.v<-exp(LogWine.ml$par[3])
@
For the new model the estimates are
\[
\sigma_1^2 = \Sexpr{sig.1} \qquad \sigma_2^2 = \Sexpr{sig.2} \qquad \sigma_V^2 = \Sexpr{sig.v}
\]
with seasonal component estimates $\Sexpr{round(xv[4:15],2)}$.

<<wine2-5,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center', echo=FALSE>>=
Y.t<-array(0,c(n,K))
for(i in 1:n){
	F.t[i,,]<-F1
	G.t[i,,]<-G1
	Sigma.t[i,,]<-sig.v
	Omega.t[i,,]<-diag(c(sig.1,sig.2,rep(0,10)))
	Y.t[i,]<-y[i]
}
#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))
G.t<-array(0,c(n,L,L))
K.t<-array(0,c(n,L,K))
L.t<-array(0,c(n,L,L))

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

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

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

	G.t[i,,]<-G1

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

	###  Second recursions
	new.r.t[i,]<-Y.t[i,] - Fmat %*% new.a.pt[i,]
	K.t[i,,]<-G.t[i,,] %*% new.P.pt[i,,] %*% tFmat %*% solve(M.t[i,,])
	new.a.pt[i+1,]<- G.t[i,,] %*% new.a.pt[i,] + 
                          matrix(K.t[i,,],nrow=L,ncol=K) %*% new.r.t[i,]
	L.t[i,,]<-G.t[i,,] -  matrix(K.t[i,,],nrow=L,ncol=K) %*% Fmat
	new.P.pt[i+1,,]<-G.t[i,,] %*% new.P.pt[i,,] %*% t(L.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):1){
	Pstar.t[i,,]<-P.tt[i,,] %*% t(G.t[i,,]) %*% solve(P.pt[i+1,,]+diag(1e-10,L))
	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,,]
}
#a.tT

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

plot(a.pt[,1],type="l",ylim=range(-2,9),
    ylab="Sales",xlab="Months",col="red",lwd=1)
lines(a.pt[,2],lty=3,col="black",lwd=1)
legend(0,5,c("Local Level","Seasonal"),
       lty=c(1,2,3),col=c("red","black"),lwd=1)

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

<<wine2-6a,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center', echo=FALSE>>=
par(mar=c(4,4,2,2))
Ypred<-as.vector(matrix(F1,nrow=1) %*% t(a.pt))
Ysmooth<-as.vector(matrix(F1,nrow=1) %*% t(a.tT))
plot(y,type="p",ylim=range(4.75,9),pch=19,cex=0.4,
        ylab="Sales",xlab="Months",lwd=1)
lines(Ypred,col='red')
lines(Ysmooth,col='blue')
legend(0,9,c("Data","Predicted", "Smoothed"),lty=c(NA,1,1),
        col=c("black","red","blue"),lwd=1,pch=c(19,NA,NA))
@

<<wine2-6b,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center',echo=FALSE>>=
Yresids<-Y.t-Ypred[-n]
par(mar=c(4,4,2,2))
plot(Yresids,type="l",ylim=range(-2,2),
    ylab="Prediction Residuals",xlab="Months",col="red",lwd=1)
lines(Yresids);points(Yresids,pch=19,cex=0.6)
@

<<wine2-6c,comment='+', fig.width=7.25, fig.height=4.75, fig.align='center',echo=FALSE>>=
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')
@
\end{document}