\documentclass[notitlepage]{article} \usepackage{../../Math545} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{bbm} \usepackage{bm} \usepackage{verbatim} \usepackage{amsfonts,amsmath,amssymb} \newtheorem{alg}{Algorithm} \newtheorem{ex}{Example} \newtheorem{note}{Note} \newenvironment{proof}{\par\noindent{\normalfont{\bf{Proof }}}} {\hfill \small$\blacksquare$\\[2mm]} \setlength{\parindent}{0in} \def\E{\mathbb{E}} \def\Var{\text{Var}} \def\Cov{\text{Cov}} \def\Corr{\text{Corr}} \def\bW{\mathbf{W}} \newcommand{\bs}[1]{\bm{#1}} \newcommand{\Rho}{\mathrm{P}} \newcommand{\cif}{\text{ if }} \newcommand{\Ra}{\Longrightarrow} \newcommand{\ra}{\longrightarrow} \def\bphi{\boldsymbol \phi} \lstloadlanguages{R} \definecolor{keywordcolor}{rgb}{0,0.6,0.6} \definecolor{delimcolor}{rgb}{0.461,0.039,0.102} \definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432} \lstdefinestyle{Rsettings}{ basicstyle=\ttfamily, breaklines=true, showstringspaces=false, keywords={if, else, function, theFunction, tmp}, % Write as many keywords otherkeywords={}, commentstyle=\itshape\color{Rcommentcolor}, keywordstyle=\color{keywordcolor}, moredelim=[s][\color{delimcolor}]{"}{"}, } \lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}} \def\ctitle{MATH 545: Introduction to Time Series Analysis} \def\ttitle{Winter 2025} \def\stitle{State Space and Dynamic Linear Models} \begin{document} <>= 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 A \textit{state space model} couples two possibly vector-valued time series processes, $\{\bX_t\}$ and $\{\bY_t\}$ to extend the types of models that can be achieved via single processes. A general linear state space model (or \textit{dynamic linear model}) is specified by two equations \[ \begin{array}{rlcll} \textrm{OBSERVATION EQUATION} \qquad\qquad& \bY_t &=& \bF_t \bX_t + \bd_t + \bV_t & \qquad K \times 1\\[6pt] \textrm{STATE EQUATION} \qquad\qquad & \bX_{t+1} &=& \bG_t \bX_t + \be_t + \bW_t & \qquad L \times 1 \end{array} \] where \begin{itemize} \item $\bF_t$ is a $K \times L$ deterministic matrix; \smallskip \item $\bd_t$ is a $K \times 1$ deterministic vector; \smallskip \item $\bG_t$ is a $L \times L$ deterministic matrix; \smallskip \item $\be_t$ is a $L \times 1$ deterministic vector. \end{itemize} and the stochastic terms are $\{\bV_t\}$ and $\{\bW_t\}$ which are independent vector white noise processes. Thus $\{\bX_t\}$ is an $L \times 1$ time series process that can be thought to represent unobserved (or latent) stochastic variation, and $\{\bY_t\}$ is a $K \times 1$ process that represents the observed process. To initialize, we also specify the initial state $\bX_1$. \bigskip \textbf{EXAMPLE:} Suppose $K=L=1$, and \[ \begin{array}{rlcll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& X_t + V_t \\[6pt] \textrm{STATE EQUATION} \qquad\qquad & X_{t+1} &=& X_t + 0.01 + W_t \end{array} \] where $V_t \sim \text{Normal}(0,\sigma_V^2)$ and $W_t \sim \text{Normal}(0,\sigma_W^2)$. We choose $\sigma_V=1, \sigma_W = 0.1$ and $X_1 = 0$. <>= set.seed(3834) n<-1000;X<-Y<-rep(0,n) F<-G<-1;dvec<-0;evec<-0.01 sigW<-0.1;sigV<-1 W<-rnorm(n)*sigW V<-rnorm(n)*sigV X[1]<-0;Y[1]<-F*X[1]+dvec+V[1] for(i in 2:n){ X[i]<-G*X[i-1]+evec+W[i] Y[i]<-F*X[i]+dvec+V[i] } par(mar=c(3,3,0,0));plot(X,type='l',ylim=range(Y));points(Y,pch=19,cex=0.3,col='red');lines(X) legend(0,max(Y),c('X','Y'),col=c('black','red'),pch=c(-1,19),lty=c(1,-1)) @ It is possible to form a \textbf{deterministic} latent structure by setting $\sigma_W = 0$. <>= F<-G<-1;dvec<-0;evec<-0.01; sigV<-1;sigW<-0.0 W<-rnorm(n)*sigW;V<-rnorm(n)*sigV X[1]<-0; Y[1]<-F*X[1]+evec+V[1] for(i in 2:n){ X[i]<-G*X[i-1]+evec+W[i] Y[i]<-F*X[i]+dvec+V[i] } par(mar=c(2,3,1,0));plot(X,type='l',ylim=range(Y));points(Y,pch=19,cex=0.3,col='red');lines(X) legend(0,max(Y),c('X','Y'),col=c('black','red'),pch=c(-1,19),lty=c(1,-1)) @ \bigskip \textbf{EXAMPLE:} Suppose $K=1$ and $L=2$, with $\bX_t = (X_{1,t},X_{2(t-1)})^\top$, and \[ \begin{array}{rlcll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& [1 \; 0] \bX_t + V_t \\[6pt] \textrm{STATE EQUATION} \qquad\qquad & \bX_{t+1} &=& \begin{bmatrix} \phi_1 & \phi_2 \\ 1 & 0 \end{bmatrix}\bX_t + \bW_t \end{array} \] where $V_t \sim \text{Normal}(0,\sigma_W^2)$ and $\bW_t = (W_{1t},W_{2t})^\top$, $W_{1t} \sim \text{Normal}(0,\sigma_W^2)$, and $W_{2t} \equiv 0$ for all $t$, with $\E[V_t W_{1s}] = 0$ for all $t,s$. Note that this model specifies an $\AR(2)$ model for $\{X_{1,t}\}$, the first component of $\{\bX_t\}$, as \[ X_{1,(t+1)} = \phi_1 X_{1,t} + \phi_2 X_{1,t-1} + W_{1t}. \] with $Y_t = X_{1,t} + V_t$. \medskip In the following simulation, we choose $\phi_1 = 1.7$, $\phi_2 = -0.72$, $\sigma_W=1, \sigma_V = 0.5$ and $\bX_1 = (0,0)^\top$. The process for $\bX_{1,t}$ is stationary, but highly autocorrelated as the $\AR(2)$ polynomial is \[ \Phi(z) = 1-1.7z + 0.72 z^2 = (1-0.8z)(1-0.9z) \] so both reciprocal roots are close to one. In this case, $\{Y_t\}$ has zero mean, and $\Var[Y_t] = \gamma_X(0) + \sigma_V^2$. Then \begin{align*} \E[Y_t Y_{t+h}] & = \E[(X_{1,t}+V_t)(X_{1,t+h}+V_{t+h}) ] \\[6pt] & = \E[X_{1,t}X_{1,t+h}] + \E[ X_{1,t}V_{t+h} ] + \E[V_t X_{1,t+h}] + \E[V_t V_{t+h}] \\[6pt] & = \gamma_X(h) \end{align*} as $\E[V_t] = 0$ for all $t$, and $\E[ X_{1,t}V_{s} ] = \E[ X_{1,t}] \E[V_{s} ] = 0$ for all $t$ and $s$. Therefore the resulting process $\{Y_t\}$ is \textbf{stationary}, with autocovariance function \[ \gamma_Y(h) \equiv \gamma_X(h) \qquad |h| \geq 1. \] <>= set.seed(3834) n<-1000;X<-matrix(0,ncol=n,nrow=2) Y<-rep(0,n) F<-c(1,0) phi1<-1.7;phi2<--0.72;G<-matrix(c(phi1,1,phi2,0),2,2) sigV<-1;sigW<-0.5 W1<-rnorm(n)*sigW V<-rnorm(n)*sigV Y[1]<-(F%*%X[,1])[1,1]+V[1] for(i in 2:n){ X[,i]<-G%*%X[,i-1]+c(W1[i],0) Y[i]<-(F%*%X[,i])[1,1]+V[i] } par(mar=c(2,3,0,0));plot(X[1,],type='l',ylim=range(-max(abs(Y)),max(abs(Y)))) points(Y,pch=19,cex=0.3,col='red');lines(X) legend(800,max(abs(Y)),c('X1','Y'),col=c('black','red'),pch=c(-1,19),lty=c(1,-1)) @ <>= set.seed(3834) n<-100000 X<-matrix(0,ncol=n,nrow=2) Y<-rep(0,n) F<-c(1,0) phi1<-1.7;phi2<--0.72 G<-matrix(c(phi1,1,phi2,0),2,2) sigW<-1;sigV<-0.5 W1<-rnorm(n)*sigW V<-rnorm(n)*sigV Y[1]<-(F%*%X[,1])[1,1]+V[1] for(i in 2:n){ X[,i]<-G%*%X[,i-1]+c(W1[i],0) Y[i]<-(F%*%X[,i])[1,1]+V[i] } psi1<-ARMAtoMA(ar=c(phi1,phi2),lag.max=50) gammaX0<-sigV^2*(1+sum(psi1^2)) true.acv.Y<-gammaX0*ARMAacf(ar=c(phi1,phi2),lag.max=50) true.acv.Y[1]<-true.acv.Y[1]+sigW^2; true.acf.Y<-true.acv.Y/true.acv.Y[1] acf.Y<-acf(Y,lag.max=50,plot=FALSE)$acf rbind(acf.Y[1:7],true.acf.Y[1:7]) @ \pagebreak \textbf{EXAMPLE: LOCAL LEVEL MODEL} $K=L=1$, and \[ \begin{array}{rlcll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& M_t + V_t \\[6pt] \textrm{STATE EQUATION} \qquad\qquad & M_{t+1} &=& M_t + W_t \end{array} \] where $V_t \sim \text{Normal}(0,\sigma_V^2)$ and $W_t \sim \text{Normal}(0,\sigma_W^2)$. In the following simulation, we choose $\sigma_V=0.1, \sigma_W = 0.02$ and $M_1 = 0$. <>= set.seed(3834) n<-1000;M<-Y<-rep(0,n) G<-F<-1 sigW<-0.02;sigV<-0.1 W<-rnorm(n)*sigW V<-rnorm(n)*sigV M[1]<-0;Y[1]<-F*X[1]+V[1] for(i in 2:n){ M[i]<-G*M[i-1]+W[i] Y[i]<-F*M[i]+V[i] } par(mar=c(2,3,0,0));plot(M,type='l',ylim=range(Y));points(Y,pch=19,cex=0.3,col='red');lines(X) legend(800,2,c('X','Y'),col=c('black','red'),pch=c(-1,19),lty=c(1,-1)) @ \bigskip Here, the process $\{Y_t\}$ is \textbf{non-stationary}: we clearly have that \[ \Var[Y_t] = \Var[M_t] + \Var[W_t] = \Var[M_t] + \sigma_W^2 \] but $\{M_t\}$ is a random walk, so it has a variance that is dependent on $t$ \[ \Var[M_t] = t \sigma_V^2. \] \pagebreak \textbf{EXAMPLE: LOCAL LINEAR MODEL} $K=1$, $L=2$ and \[ \begin{array}{rlcll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& M_t + V_t \\[6pt] \textrm{STATE EQUATION} \qquad\qquad & M_{t+1} &=& M_t + B_t + W_{1t}\\[3pt] & B_{t+1} &=& B_t + W_{2t} \end{array} \] where we define the state $\bX_t = (M_t,B_t)^\top$, and $\bW_t = (W_{1t},W_{2t})^\top$, so that \[ \bG_t = \bG_L = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} \qquad \bF_t = [1 \; 0]. \] We have $V_t \sim \text{Normal}(0,\sigma_V^2)$ and $W_{1t} \sim \text{Normal}(0,\sigma_{W_1}^2), W_{2t} \sim \text{Normal}(0,\sigma_{W_2}^2)$ with $W_{1t}$ and $W_{2t}$ independent. In the following simulation, we choose $\sigma_V=5, \sigma_{W_1} = 0.001$ and $\sigma_{W_2} = 0.005$, with $M_1 = B_1 = 0$. <>= set.seed(3834) n<-1000;X<-matrix(0,ncol=n,nrow=2);Y<-rep(0,n) F<-c(1,0) G<-matrix(c(1,0,1,1),2,2) sigV<-5;sigW1<-0.001;sigW2<-0.005 V<-rnorm(n)*sigV;W1<-rnorm(n)*sigW1;W2<-rnorm(n)*sigW2 X[,1]<-0; Y[1]<-(F%*%X[,1])[1,1]+V[1] for(i in 2:n){ X[,i]<-G%*%X[,i-1]+c(W1[i],W2[i]) Y[i]<-(F%*%X[,i])[1,1]+V[i] } mt<-apply(X,2,sum) par(mar=c(2,3,1,0));plot(mt,type='l',ylim=range(Y));points(Y,pch=19,cex=0.3,col='red') lines(mt);legend(0,max(Y),c('X','Y'),col=c('black','red'),pch=c(-1,19),lty=c(1,-1)) @ \bigskip This model is suitable for carrying out non-linear regression; the smoothness of the underlying function as modelled by the states is determined by the magnitude of $\sigma_{W_1}$ and $\sigma_{W_2}$. \pagebreak \textbf{EXAMPLE: SEASONAL MODEL} In a deterministic purely seasonal model, with seasonality $r$ lags, we have $Y_t = s_t + V_t$, say, where $s_t = s_{t-r}$ for all $t$, such that $\sum\limits_{t=1}^r s_t = 0$. That is, \[ s_{t+1} = -s_t - s_{t-1} - \cdots - s_{t-r+2}. \] To construct a stochastic seasonal model, we may write \[ S_{t+1} = -S_t - S_{t-1} - \cdots - S_{t-r+2} + W_{1t} \] where $\{W_{1t}\} \sim \WN(0,\sigma_S^2)$, say. We write this as a state space model as follows: we have $K=1$, $L=(r-1)$ and \[ \begin{array}{rlcll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& \bF_t \bX_t + V_t \\[6pt] \textrm{STATE EQUATION} \qquad\qquad & \bX_{t+1} &=& \bG_t \bX_t + \bW_t \end{array} \] where we define the state $\bX_t = (Y_t,Y_{t-1},\ldots,Y_{t-r+2})^\top$, and $\bW_t = (S_t,0,0,\ldots,0)^\top$, so that \[ \bG_t \equiv \bG_S = \begin{bmatrix} -1 & -1 & -1 & -1 & \cdots & -1 &-1\\[3pt] 1 & 0 & 0 & 0 & \cdots & 0 &0\\[3pt] 0 & 1 & 0 & 0 & \cdots & 0 &0\\[3pt] 0 & 0 & 1 & 0 & \cdots & 0 &0\\[3pt] \vdots & \vdots & \vdots & \ddots & \ddots & \ddots &\vdots\\[3pt] 0 & 0 & 0 & 0 & \cdots & 1 &0 \end{bmatrix} \] and \[ \bF_t = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \end{bmatrix} \] with $V_t \sim \text{Normal}(0,\sigma_V^2)$ and $S_{t} \sim \text{Normal}(0,\sigma_{S}^2)$. In the following simulation, with $r=12$, we choose $\sigma_V=0.05, \sigma_{S} = 0.2$, and with $\bX_t = (-0.1,0.1,0.3,0.7,0.4,0.2,-0.1,0.0,-0.1,0.0,-1.4)^\top$. <>= set.seed(3834);n<-200;r<-12 X<-matrix(0,ncol=n,nrow=r-1);Y<-rep(0,n) F<-c(1,rep(0,r-2)) G<-matrix(0,r-1,r-1) G[1,]<-rep(-1,r-1);for(j in 2:(r-1)){G[j,j-1]<-1} sigV<-0.05;sigS<-0.2 V<-rnorm(n)*sigV;S<-rnorm(n)*sigS X[,1]<-c(-0.1,0.1,0.3,0.7,0.4,0.2,-0.1,0.0,-0.1,0.0,-1.4) Y[1]<-(F%*%X[,1])[1,1]+V[1] for(i in 2:n){ Wi<-c(S[i],rep(0,r-2)) X[,i]<-G%*%X[,i-1]+Wi Y[i]<-(F%*%X[,i])[1,1]+V[i] } par(mar=c(2,3,0,0));plot(X[1,],type='l',ylim=range(Y));points(Y,pch=19,cex=0.3,col='red') legend(0,max(Y),c('X','Y'),col=c('black','red'),pch=c(-1,19),lty=c(1,-1)) @ \bigskip \textbf{EXAMPLE: LOCAL LINEAR + SEASONAL MODEL} We can combine the two previous models to construct a state space model as follows: we have $K=1$, $L=(r+1)$ and \[ \begin{array}{rlcll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& \bF_t \bX_t + V_t \\[6pt] \textrm{STATE EQUATION} \qquad\qquad & \bX_{t+1} &=& \bG_t \bX_t + \bW_t \end{array} \] where we define the state $\bX_t = (M_t, B_t, S_t,S_{t-1},\ldots,S_{t-r+2})^\top$, and $\bW_t = (W_{1t}, W_{2t}, W_{3t},0,\ldots,0)^\top$, so that \[ \bG_t \equiv \bG = \begin{bmatrix} \bG_L & \mathbf{0}\\[3pt] \mathbf{0} & \bG_S \end{bmatrix} \] and \[ \bF_t = \begin{bmatrix} 1 & 0 & 1 & 0 & \cdots & 0 \end{bmatrix} \] with $V_t \sim \text{Normal}(0,\sigma_V^2)$ and $W_{3t} \sim \text{Normal}(0,\sigma_{S}^2)$. In the following simulation, we choose $\sigma_V=1, \sigma_{W_1} = 0.1$, $\sigma_{W_2} = 0.05$, and $\sigma_S=0.5$, with $M_1 = B_1 = 0$, and with $\bX_1 = (-0.1,0.1,0.3,0.7,0.4,0.2,-0.1,0.0,-0.1,0.0,-1.4)^\top$. <>= set.seed(3834) n<-200;r<-12 X<-matrix(0,ncol=n,nrow=r+1);Y<-mt<-rep(0,n) F<-c(1,0,1,rep(0,r-2)) G<-matrix(0,r+1,r+1) G[1:2,1:2]<-matrix(c(1,0,1,1),2,2) G[3,]<-c(0,0,rep(-1,r-1));for(j in 4:(r+1)){G[j,j-1]<-1} sigV<-1 sigW1<-0.1 sigW2<-0.05 sigS<-0.5 V<-rnorm(n)*sigV;W1<-rnorm(n)*sigW1;W2<-rnorm(n)*sigW2;S<-rnorm(n)*sigS X[,1]<-c(0,0,-0.1,0.1,0.3,0.7,0.4,0.2,-0.1,0.0,-0.1,0.0,-1.4) mt[1]<-(F%*%X[,1])[1,1] Y[1]<-mt[1]+V[1] for(i in 2:n){ Wi<-c(W1[i],W2[i],S[i],rep(0,r-2)) X[,i]<-G%*%X[,i-1]+Wi mt[i]<-(F%*%X[,i])[1,1] Y[i]<-mt[i]+V[i] } par(mar=c(2,3,0,0)) plot(mt,type='l',ylim=range(Y));points(Y,pch=19,cex=0.3,col='red') lines(X[1,],col='blue') legend(0,max(Y),c(expression(X[t]),expression(M[t]),expression(Y[t])), col=c('black','blue','red'), pch=c(-1,-1,19),cex=c(1,1,1),lty=c(1,1,-1)) @ \end{document}