\documentclass[11pt]{article} \usepackage{fancyhdr} \usepackage{lastpage} \usepackage{amsfonts,amssymb,amsmath} \usepackage{undertilde} \usepackage{../MATH545} \usepackage{enumitem} \usepackage{booktabs} \usepackage{bm} \usepackage{hyperref} \parindent 0in \pagestyle{empty} \def\ExNumber{3} \lfoot{MATH 545 EXERCISES \ExNumber~Solutions} \cfoot{ } \rfoot{{\small{\em Page \thepage \ of \pageref{LastPage}}}} \renewcommand{\headrulewidth}{0.0pt} \renewcommand{\footrulewidth}{0.0pt} \def\ctitle{MATH 545: Introduction to Time Series Analysis} \def\ttitle{Winter 2025} \def\stitle{Exercises \ExNumber: Solutions} \begin{document} \pagestyle{fancy} \maketitle <>= library(knitr) library(kableExtra) suppressMessages(library(dplyr)) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) options(repos=c(CRAN="https://cloud.r-project.org/")) @ <>= set.seed(11) myrnd <- function(x,d) {formatC(x, format="f", digits=d)} @ \begin{enumerate}[label=\arabic*.,itemsep=0.2in] \item \textit{The data set contained in the file \texttt{LakeErie.txt} is a collection of 600 monthly measurements of the level of Lake Erie, from January 1921 to December 1970.} <>= erie.url<-'https://www.math.mcgill.ca/dstephens/545/Data/LakeErie.txt' erie<-scan(erie.url) par(mar=c(4,4,2,2)) plot(erie,type='l',xlab='t',ylab='Level (ft)') points(erie,pch=19,cex=0.5) @ \medskip \begin{enumerate}[label=(\alph*)] \item \textit{Carry out an analysis of these data in \texttt{R} by examining ACF and PACF plots of the original and and any appropriately differenced series.} \item \textit{Attempt to characterize the stochastic (random) component of the data using a SARIMA model. } \item \textit{Using the code provided in \texttt{knitr} sheets, formulate and fit a seasonal state space model to these data.} \end{enumerate} \pagebreak \textbf{Solutions:} \begin{enumerate}[label=(\alph*),itemsep=0.2in] \item We can assess whether there is any global trend by using \texttt{lm} on the original data <>= x<-c(1:600) y<-erie summary(lm(y~x))$coef @ It seems that there is a global increasing trend as the slope coefficient is significantly different from zero. This linear trend can be removed by first differencing. <>= erie.diff<-diff(erie) par(mar=c(4,4,2,2)) plot(erie.diff,type='l',xlab='t',ylab='Differenced data') points(erie.diff,pch=19,cex=0.5) @ This series is more plausibly stationary, so we examine the ACF <>= par(mar=c(4,4,3,2)) acf(erie.diff,lag.max=50,main='ACF of Erie data') @ This suggests an analysis of the lag-12 differenced data may be more appropriate. <>= erie.diff<-diff(erie.diff, 12) par(mar=c(4,4,3,2)) plot(erie.diff,type='l',xlab='t',ylab='Seasonally differenced') points(erie.diff,pch=19,cex=0.5) acf(erie.diff,lag.max=50,main='ACF of lag 1/lag 12 differenced Erie data') @ We can conclude that if $\{X_t\}$ denotes the original data, then \[ Y_t = (1-B^{12}) (1-B) X_t \] is a plausibly stationary series, albeit one that still exhibits some dependence at lag 1 and lag 12. \item \textit{Attempt to characterize the stochastic (random) component of the data using a SARIMA model. } Given the final ACF plot from (a), it seems that the $\SARIMA(1,1,1) \times (1,1,1)_{12}$ model, or some simpler related model may provide a good fit. We try fitting several SARIMA models, and compare using the AIC and BIC model selection criteria. <>= fit.1111 <- arima(erie,c(1,1,1), seasonal=list(order=c(1,1,1),period=12)) fit.1101 <- arima(erie,c(1,1,1), seasonal=list(order=c(0,1,1),period=12)) fit.1001 <- arima(erie,c(1,1,0), seasonal=list(order=c(0,1,1),period=12)) fit.1010 <- arima(erie,c(1,1,0), seasonal=list(order=c(1,1,0),period=12)) fit.0101 <- arima(erie,c(0,1,1), seasonal=list(order=c(0,1,1),period=12)) AIC(fit.1111,fit.1101,fit.1001,fit.1010,fit.0101) AIC(fit.1111,fit.1101,fit.1001,fit.1010,fit.0101,k=log(length(erie))) @ The results are as follows: \begin{table}[ht] \centering \begin{tabular}{c|c|c} Model & AIC & BIC \\ \hline $\SARIMA(1,1,1) \times (1,1,1)_{12}$ & \Sexpr{AIC(fit.1111)} & \Sexpr{AIC(fit.1111,k=log(length(erie)))} \\ $\SARIMA(1,1,1) \times (0,1,1)_{12}$ & \Sexpr{AIC(fit.1101)} & \Sexpr{AIC(fit.1101,k=log(length(erie)))} \\ $\SARIMA(1,1,0) \times (0,1,1)_{12}$ & \Sexpr{AIC(fit.1001)} & \Sexpr{AIC(fit.1001,k=log(length(erie)))} \\ $\SARIMA(1,1,0) \times (1,1,0)_{12}$ & \Sexpr{AIC(fit.1010)} & \Sexpr{AIC(fit.1010,k=log(length(erie)))} \\ $\SARIMA(0,1,1) \times (0,1,1)_{12}$ & \Sexpr{AIC(fit.0101)} & \Sexpr{AIC(fit.0101,k=log(length(erie)))} \\ \hline \end{tabular} \end{table} and hence the selected model seems to be $\SARIMA(1,1,0) \times (0,1,1)_{12}$ as this has the lowest AIC and BIC. <>= coef(fit.1001) @ If we take the residuals from this model, we see that the white noise requirement is approximately met <>= acf(residuals(fit.1001),lag.max=30,main='ACF of residuals') @ \item \textit{Using the code provided in \texttt{knitr} sheets, formulate and fit a seasonal state space model to these data.} It seems here that a local-level model plus a stochastic seasonal component with period 12 should be adequate. That is, \[ 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)$. <>= y<-erie n<-length(y) 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,L-2)) G1<-matrix(0,L,L) G1[1,1]<-1 G1[2,2:12]<--1;for(i in 2:11){G1[i+1,i]<-1} @ <>= #Perform the optimization 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,Lval-2))) 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)) } @ <>= ##Compute Estimates xv<-c(0,0,0,y[1],rnorm(11,10,0.2)) Erie.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<-Erie.ml$par Erie.ml X.t<-matrix(0,nrow=length(y),ncol=length(Erie.ml$par)-3) X.t[1,]<-Erie.ml$par[4:length(Erie.ml$par)] sig.1<-exp(Erie.ml$par[1]) sig.2<-exp(Erie.ml$par[2]) sig.v<-exp(Erie.ml$par[3]) @ For this 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)}. \] <>= 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,L-2))) 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,L-2))) new.a.pt[1,]<-X.t[1,] new.P.pt[1,,]<-diag(c(sig.1,sig.2,rep(0,L-2))) #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,20), ylab="Level",xlab="Months",col="red",lwd=1) lines(a.pt[,2],lty=3,col="black",lwd=1) legend(400,10,c("Local Level","Seasonal"), lty=c(1,2,3),col=c("red","black"),lwd=1) plot(a.pt[,2],type="l",ylim=range(-3,3), ylab="Level",xlab="Months",col="black",lwd=1) legend(400,3,c("Seasonal component"),col=c("black"),lwd=1) @ <>= 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(5,25),pch=19,cex=0.4, ylab="Level",xlab="Months",lwd=1) lines(Ypred,col='red') lines(Ysmooth,col='blue') legend(0,25,c("Data","Predicted", "Smoothed"),lty=c(NA,1,1), col=c("black","red","blue"),lwd=1,pch=c(19,NA,NA)) @ <>= 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) @ <>= 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{enumerate} \item \textit{Consider the following state space model for scalar observations $Y_{1:n}$ and states $X_{1:n}$} \begin{eqnarray*} Y_{t} & = & \beta X_t + V_t\\ X_{t+1} & = & \phi X_t + W_t \end{eqnarray*} \textit{where $\{V_t\}$ and $\{W_t\}$ are independent Gaussian white noise series with variances $\sigma_V^2$ and $\sigma_W^2$ respectively. To initiate the model, suppose that $X_1 \sim \text{Normal}(0,\sigma_W^2/(1-\phi^2))$. } \begin{enumerate}[label=(\alph*)] \item \textit{Under what conditions is $\{Y_t\}$ stationary ?} \textbf{Solution:} Provided $\phi \neq 1$, we have that $\E[Y_t] = \E[X_t] = 0$, \[ \Var[Y_t] = \beta^2 \Var[X_t] + \sigma_V^2 = \frac{\beta^2}{1-\phi^2} \sigma_W^2 + \sigma_V^2 \] which does not depend on $t$, and that for $h\neq0$ \begin{align*} \Cov[Y_t,Y_{t+h}] & = \Cov[(\beta X_t + V_t),(\beta X_{t+h} + V_{t+h})]\\[6pt] & = \beta^2 \Cov[X_t,X_{t+h}] & \textrm{other terms zero}. \end{align*} Thus $\{Y_t\}$ is stationary if $\{X_t\}$ is stationary, which holds provided $\phi \neq 1$. \item \textit{Assuming $\{Y_t\}$ is stationary, derive its ACVF, $\gamma_Y(h)$}. \medskip \textbf{Solution:} In this case, from properties of the stationary $\AR(1)$ process, we have from part (a) that \[ \gamma_Y(h) = \frac{\beta^2 \phi^h}{1-\phi^2} \sigma_W^2 \] for $h=0,\pm 1,\pm 2,\ldots$. \item\textit{ Derive the Kalman filtering recursions for this model, based on the filtering and prediction distributions $p(x_t|y_{1:t})$ and $p(x_{t+1}|y_{1:t})$. Compute also the prediction distribution $p(y_{t+1}|y_{1:t})$.} \medskip \textbf{Solutions:} We have two components of the recursion. Using the notation and results from lectures, we have for the scalar state process $\{X_t\}$ the filtering and prediction distributions \[ p(x_t|y_{1:t}) \equiv Normal(a_{t|t}, P_{t|t}) \qquad p(x_{t+1}|y_{1:t}) \equiv Normal(a_{t+1|t}, P_{t+1|t}). \] \begin{itemize} \item \textbf{\textit{Filtering:}} we have that \[ p(x_t|y_{1:t}) \propto p(y_t|x_t) p(x_t|y_{1:(t-1)}) \] where \[ p(y_t|x_t) \equiv Normal(\beta x_t, \sigma_V^2) \qquad p(x_t|y_{1:(t-1)}) \equiv Normal(a_{t|t-1}, P_{t|t-1}) \] In the relevant quadratic component of the calculation, we have \[ \frac{1}{\sigma_V^2} (y_t - \beta x_t)^2 + \frac{1}{P_{t|t-1}} (x_t - a_{t|t-1})^2 \] and rearranging, this yields \[ \frac{\beta^2}{\sigma_V^2} (x_t - y_t/\beta)^2 + \frac{1}{P_{t|t-1}} (x_t - a_{t|t-1})^2. \] Using the standard complete-the-square formula \[ A(x-a)^2 + B(x-b)^2 = (A+B)\left(x - \frac{Aa+Bb}{A+B} \right)^2 + \frac{AB}{A+B}(a-b)^2 \] may conclude that \[ a_{t|t} = \frac{(\beta y_t /\sigma_V^2) + (a_{t|t-1}/P_{t|t-1})}{(\beta^2/\sigma_V^2) + (1/P_{t|t-1})} \] and \[ P_{t|t} = \left( \frac{\beta^2}{\sigma_V^2} + \frac{1}{P_{t|t-1}} \right)^{-1} \] \medskip \item \textbf{\textit{Prediction:}} We have that \[ p(x_{t+1}|y_{1:t}) = \int p(x_{t+1}|x_t) p(x_t|y_{1:t}) \ d x_t. \] Taking expectations in the state equation, we have that \[ a_{t+1|t} = \E_{X_t|Y_{1:t}}[\phi X_t + W_t|y_{1:t}] = \phi \E_{X_t|Y_{1:t}}[X_t|y_{1:t}] = \phi a_{t|t} \] and \[ P_{t+1|t} = \Var_{X_t|Y_{1:t}}[\phi X_t + W_t|y_{1:t}] = \phi^2 \Var_{X_t|Y_{1:t}}[X_t|y_{1:t}] + \sigma_W^2 = \phi^2 P_{t|t-1} + \sigma_W^2. \] \medskip Finally, from the observation equation, we have that \[ \E_{Y_{t+1}|Y_{1:t}}[Y_{t+1}|y_{1:t}] = \E_{Y_{t+1}|Y_{1:t}}[\beta X_{t+1} + V_t|y_{1:t}] = \beta a_{t+1|t} \] and \[ \Var_{Y_{t+1}|Y_{1:t}}[Y_{t+1}|y_{1:t}] = \Var_{Y_{t+1}|Y_{1:t}}[\beta X_{t+1} + V_t|y_{1:t}] = \beta^2 \Var_{Y_{t+1}|Y_{1:t}}[X_{t+1}|y_{1:t}] + \sigma_V^2 = \beta^2 P_{t+1|t} + \sigma_W^2 \] and so \[ p(y_{t+1}|y_{1:t}) \equiv Normal(\beta a_{t+1|t}, \beta^2 P_{t+1|t} + \sigma_W^2) \] \end{itemize} \item \textit{Using code similar to that provided in \texttt{knitr} sheets, or otherwise, perform maximum likelihood estimation of the system parameters for the data in the \texttt{Ex3-Q2data.txt}} \begin{center} \href{https://www.math.mcgill.ca/dstephens/545/Data/Ex3-Q2data.txt}{\texttt{https://www.math.mcgill.ca/dstephens/545/Data/Ex3-Q2data.txt}} \end{center} \end{enumerate} \textbf{Solution:} <>= Q2.url<-'https://www.math.mcgill.ca/dstephens/545/Data/Ex3-Q2data.txt' y<-scan(Q2.url) par(mar=c(4,4,2,2)) plot(y,type='l',xlab='t') points(y,pch=19,cex=0.7) n<-length(y) K<-1;L<-1 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)) @ <>= #Perform the optimization ml.filter.new<-function(x,nv,Lval,Kval,yv){ Yt<-yv Sigmat<-exp(x[2]) Omegat<-exp(x[1]) 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] Ppt[1,,]<-Omegat Fm<-matrix(x[3],1,1) Gm<-matrix(x[4],1,1) #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)) } @ <>= ##Compute Estimates xv<-c(0.1,0.1,0,0.5,0.5) Q2.ml<-optim(xv,ml.filter.new,nv=n,Kval=K,Lval=L,yv=y, method="BFGS",control=list(trace=10,fnscale=-500,maxit=500,abstol=1e-6)) xv<-Q2.ml$par Q2.ml X.t<-matrix(0,nrow=length(y),ncol=1) X.t[1,]<-Q2.ml$par[5] sig.W<-exp(Q2.ml$par[1]) sig.V<-exp(Q2.ml$par[2]) beta.hat<-Q2.ml$par[3] phi.hat<-Q2.ml$par[4] @ For this model the estimates are $\beta = \Sexpr{beta.hat}$, $\phi = \Sexpr{phi.hat}$ \[ \sigma_W^2 = \Sexpr{sig.W} \qquad \sigma_2^2 = \Sexpr{sig.V} \] with initial value $X_1 = \Sexpr{X.t[1,1]}$. <>= Y.t<-array(0,c(n,K)) for(i in 1:n){ F.t[i,1,1]<-beta.hat G.t[i,,]<-phi.hat Sigma.t[i,,]<-sig.V Omega.t[i,,]<-sig.W 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,,]<-sig.W new.a.pt[1,]<-X.t[1,] new.P.pt[1,,]<-sig.W #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,,]<-phi.hat 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(-8,8), ylab=expression(X[t]),xlab="t",col="red",lwd=1) @ <>= par(mar=c(4,4,2,2)) Ypred<-as.vector(matrix(beta.hat,nrow=1) %*% t(a.pt)) Ysmooth<-as.vector(matrix(beta.hat,nrow=1) %*% t(a.tT)) plot(y,type="p",ylim=range(-8,8),pch=19,cex=0.4, ylab="Y",xlab="t",lwd=1) lines(Ypred,col='red') lines(Ysmooth,col='blue') legend(0,25,c("Data","Predicted", "Smoothed"),lty=c(NA,1,1), col=c("black","red","blue"),lwd=1,pch=c(19,NA,NA)) @ <>= Yresids<-Y.t-Ypred[-n] par(mar=c(4,4,2,2)) plot(Yresids,type="l",ylim=range(-8,8), ylab="Prediction Residuals",xlab="t",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') @ \pagebreak \item \textit{Formulate the stationary $\ARMA(2,2)$ model \[ (1-0.5 B) (1-0.9B) X_t = (1+B) (1+2B) Z_t \] as a state space model.} \medskip \textit{If it was required to add a \textbf{deterministic} seasonal component, $\{S_t\}$, with period 4 to this model, so that say \[ Y_t = X_t + S_t \] how would the state space model be amended? Identify the corresponding $G$ and $F$ matrices, and the other system components.} \medskip \textbf{Solution:} We can write the model \[ X_t -1.4 X_{t-1} + 0.45 X_{t-2} = Z_t + 3 Z_{t-1} + 2 Z_{t-2} \] The formulation for $\ARMA(2,2)$ is given as follows: let $m=\max\{p,q+1\} = 3$, and \begin{align*} Y_t &= \bF_t \bX_t + V_t \\[6pt] \bX_{t+1} &= \bG_t \bX_t + \bW_t \end{align*} where we define the state $\bX_t = (X_{1,t-2},X_{1,t-1},\ldots,X_{1,t})^\top$, and $\bW_t = (0,0,0,\ldots,Z_{t+1})^\top$, so that \[ \bG_t = \begin{bmatrix} 0 & 1 & 0 \\[3pt] 0 & 0 & 1 \\[3pt] \phi_3 & \phi_{2} & \phi_{1} \end{bmatrix} \qquad \bF_t = \begin{bmatrix} \theta_{2} & \theta_{1} & \theta_0 \end{bmatrix}. \] with $V_t \equiv 0$ and $Z_{t} \sim \text{Normal}(0,\sigma_{Z}^2)$ and where \[ (\phi_1,\phi_2,\phi_3) = (1.4,-0.45,0) \qquad (\theta_0,\theta_1,\theta_2) = (1,3,2). \] The following code verifies the simulation <>= set.seed(3834);n<-10000;p<-2;q<-2;m<-max(p,q+1) phi.vec<-theta.vec<-rep(0,m) phi.vec[1:p]<-c(1.4,-0.45) theta.vec[1:q]<-c(3,2) X<-matrix(0,ncol=n,nrow=m);Y<-mt<-rep(0,n) F<-c(rev(theta.vec[1:(m-1)]),1) G<-matrix(0,m,m);for(j in 1:(m-1)){G[j,j+1]<-1} G[m,]<-rev(phi.vec) sigV<-0;sigZ<-1 V<-rnorm(n)*sigV;Z<-rnorm(n)*sigZ X[,1]<-rep(0,m) Y[1]<-(F%*%X[,1])[1,1]+V[1] for(i in 2:n){ Wi<-c(rep(0,m-1),Z[i]) 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[1:500],type='l',ylim=range(Y)); points(Y[1:500],pch=19,cex=0.3) legend(0,max(Y),c('X','Y'),col=c('black','black'),pch=c(-1,19),lty=c(1,-1)) @ <>= psi.vec<-ARMAtoMA(ar=phi.vec,ma=theta.vec,lag.max=30) varY<-1+sum(psi.vec^2) true.acf<-ARMAacf(ar=phi.vec,ma=theta.vec,lag.max=30) true.acv<-true.acf*varY varY;var(Y) par(mar=c(4,4,0,0)) acf(Y,lag.max=30,main='ACF of Y');points(0:30,true.acf,col='red',pch=19) @ The simulation using the state space form matches the target ACF. \medskip To add the seasonal component, merely add the additional components to the state vector. In the deterministic seasonal model, we have that \[ S_{t+1} = -S_t - S_{t-1} - S_{t-2} \] so the state equation becomes \[ \begin{bmatrix} X_{1,t-1} \\[3pt] X_{1,t} \\[3pt] X_{1,t+1} \\[3pt] S_{t+1} \\[3pt] S_t \\[3pt] S_{t-1} \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 \\[3pt] 0 & 0 & 1 & 0 & 0 & 0 \\[3pt] 0 & \phi_2 & \phi_1 & 0 & 0 & 0 \\[3pt] 0 & 0 & 0 & -1 & -1 & -1 \\[3pt] 0 & 0 & 0 & 1 & 0 & 0 \\[3pt] 0 & 0 & 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} X_{1,t-2} \\[3pt] X_{1,t-1} \\[3pt] X_{1,t} \\[3pt] S_{t} \\[3pt] S_{t-1} \\[3pt] S_{t-2} \end{bmatrix} + \begin{bmatrix} 0 \\[3pt] 0 \\[3pt] Z_t \\[3pt] 0 \\[3pt] 0 \\[3pt] 0 \end{bmatrix} \] with $\bF = \begin{bmatrix} \theta_{2} & \theta_{1} & \theta_0 & 1 & 0 & 0\end{bmatrix}$ <>= set.seed(3834);n<-10000;p<-2;q<-2;m<-max(p,q+1) phi.vec<-theta.vec<-rep(0,m) phi.vec[1:p]<-c(1.4,-0.45) theta.vec[1:q]<-c(3,2) X<-matrix(0,ncol=n,nrow=m+3);Y<-mt<-rep(0,n) F<-c(rev(theta.vec[1:q]),1,1,0,0) G<-matrix(0,m+3,m+3);for(j in 1:(m-1)){G[j,j+1]<-1} G[m,]<-rev(phi.vec) G[m+1,m+c(1:3)]<--1 G[m+2,m+1]<-1 G[m+3,m+2]<-1 sigV<-0;sigZ<-1 V<-rnorm(n)*sigV;Z<-rnorm(n)*sigZ X[,1]<-c(rep(0,m),5,-5,10) Y[1]<-(F%*%X[,1])[1,1]+V[1] for(i in 2:n){ Wi<-c(rep(0,m-1),Z[i],rep(0,3)) 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[1:500],type='l',ylim=range(Y)); points(Y[1:500],pch=19,cex=0.3) legend(100,max(Y),c('X','Y'),col=c('black','black'),pch=c(-1,19),lty=c(1,-1)) @ <>= par(mar=c(2,3,3,0)) acf(Y,lag.max=30,main='ACF of Y') @ \item \textit{The following data from the \texttt{astsa} library comprise} \begin{itemize} \item \texttt{cmort}: \textit{average weekly cardiovascular mortality in Los Angeles County, US. The data are 508 six-day averages obtained from daily values over the period 1970-1979;} \item \texttt{tempr}: \textit{the corresponding weekly temperature data for the same period;} \item \texttt{part}: \textit{a corresponding particulate matter series for the same period.} \end{itemize} \textit{studied in the book by Shumway et al.} <>= suppressMessages(library(vars));library(astsa) x = cbind(cmort, tempr, part) par(mar=c(4,4,0,2)) plot.ts(x , main = "", xlab = "") @ \begin{enumerate}[label=(\alph*)] \item \textit{Using the \texttt{VAR} function from the \texttt{vars} library, fit $\VAR(p)$ models for various $p$, and compare them using the AIC.} \medskip \textbf{Solution: } The ACF of the series indicates non-stationarity. <>= par(mar=c(4,4,0,2)) acf(x,lag.max=30,main='ACF of original series') @ The following code performs automatic selection of the model order <>= par(mar=c(4,4,0,2)) acf(x,lag.max=30,main='ACF of original series') fit1<-VAR(x,lag.max=12,type='const',ic='AIC') summary(fit1) @ The selected model under AIC is $p=9$. Under BIC (or Schwarz criterion), the selected model order is $p=2$. <>= fit2<-VAR(x,lag.max=12,type='const',ic='SC') summary(fit2) @ \item \textit{Interpret the model selected using AIC, and check the fit of this model by inspecting the residual series.} The $p=9$ selected by AIC suggests a complex pattern of cross-correlation. <>= acf(residuals(fit1),lag.max=30) @ The $p=2$ selected by BIC is simpler, but apparently still gives a reasonable fit, and residuals that are largely uncorrelated. <>= acf(residuals(fit2),lag.max=30) @ \item\textit{ Assess whether seasonal differencing of the series improves the quality of fit.} <>= xd<-diff(x,12) par(mar=c(4,4,0,2)) plot.ts(xd , main = "", xlab = "") acf(xd,lag.max=30,main='ACF differencing at lag 12') fit3<-VAR(xd,lag.max=12,type='const',ic='AIC') fit3$p fit4<-VAR(xd,lag.max=12,type='const',ic='SC') fit4$p @ Thus it seems seasonal difference does not change the complexity of the inferred models, in part because the ACF indicates remaining non-stationarity. First differencing, however, does seem to have an impact <>= xd1<-diff(x,1) par(mar=c(4,4,0,2)) plot.ts(xd1 , main = "", xlab = "") acf(xd1,lag.max=30,main='ACF differencing at lag 1') fit5<-VAR(xd1,lag.max=12,type='const',ic='AIC') fit5$p fit6<-VAR(xd1,lag.max=12,type='const',ic='SC') fit6$p @ The conclusion from this analysis is that the series should be first differenced, and this effectively removes the apparent seasonal structure, but that seasonal differencing does not really simplify the data structure. \end{enumerate} \end{enumerate} \end{document}