\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}}} \def\ctitle{MATH 545: Introduction to Time Series Analysis} \def\ttitle{Winter 2025} \def\stitle{The Kalman Filter for ARMA 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 The Kalman filter can be used for inference for $\ARMA(p,q)$ models if we specify the two equations as follows: let $m=\max\{p,q+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 = (X_{1,t-m+1},X_{1,t-m+2},\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 & \cdots & 0 &0\\[3pt] 0 & 0 & 1 & \cdots & 0 &0\\[3pt] \vdots & \vdots & \vdots & \ddots & \vdots &\vdots\\[3pt] 0 & 0 & 0 & \cdots & 1 &0\\[3pt] 0 & 0 & 0 & \cdots & 0 &1\\[3pt] \phi_m & \phi_{m-1} & \phi_{m-3} & \cdots & \phi_2 &\phi_1 \end{bmatrix} \qquad \bF_t = \begin{bmatrix} \theta_{m-1} & \theta_{m-2} & \cdots & \theta_{1} & \theta_0 \end{bmatrix}. \] with $V_t \equiv 0$ and $Z_{t} \sim \text{Normal}(0,\sigma_{Z}^2)$. \bigskip \textbf{EXAMPLE:} Consider the $\ARMA(1,1)$ case \[ Y_t - \phi Y_{t-1} = Z_t + \theta Z_{t-1} \qquad t \in \R \] so that $m=\max\{1,2 \} = 2$, and \[ \bG_t = \begin{bmatrix} 0 & 1 \\[3pt] 0 &\phi\\[3pt] \end{bmatrix} \qquad \bF_t = \begin{bmatrix} \theta & 1 \end{bmatrix} . \] as $\phi_2 \equiv 0$, and $\theta_0 = 1$. Therefore \[ \begin{array}{rccll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& \begin{bmatrix} \theta & 1 \end{bmatrix} \begin{bmatrix} X_{1,t-1} \\ X_{1,t} \end{bmatrix} + V_t \\[12pt] \\ \textrm{STATE EQUATION} \qquad\qquad & \begin{bmatrix} X_{1,t} \\[3pt] X_{1,t+1} \end{bmatrix} &=& \begin{bmatrix} 0 & 1 \\[3pt] 0 &\phi_1 \end{bmatrix} \begin{bmatrix} X_{1,t-1} \\[3pt] X_{1,t} \end{bmatrix}+ \begin{bmatrix} 0 \\[3pt] Z_{t+1} \end{bmatrix} \end{array} \] From the first equation, we have \[ Y_t = X_{1,t} + \theta X_{1,t-1} = (1+\theta \B) X_t \] and from the second equation the identity $X_{1,t} = X_{1,t}$, and \[ X_{1,t+1} = \phi_1 X_{1,t} + Z_{t+1} \] that is $(1-\phi) X_{1,t+1} = Z_{t+1}$ or equivalently $(1-\phi \B) X_{1,t} = Z_{t}$, so that \[ X_{1,t} = (1-\phi \B)^{-1} Z_{t} \] Hence \[ Y_t = (1+\theta \B) X_t = (1+\theta \B) (1-\phi \B)^{-1} Z_{t} \] as required. As a special case with $\phi = 0$, we have the $\MA(1)$ model where \[ \begin{array}{rccll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& \begin{bmatrix} \theta & 1 \end{bmatrix} \begin{bmatrix} X_{1,t-1} \\[3pt] X_{1,t} \end{bmatrix} + V_t \\[12pt] \\ \textrm{STATE EQUATION} \qquad\qquad & \begin{bmatrix} X_{1,t} \\[3pt] X_{1,t+1} \end{bmatrix} &=& \begin{bmatrix} 0 & 1 \\[3pt] 0 &0 \end{bmatrix} \begin{bmatrix} X_{1,t-1} \\[3pt] X_{1,t} \end{bmatrix}+ \begin{bmatrix} 0 \\[3pt] Z_{t+1} \end{bmatrix} \end{array} \] or equivalently \[ \begin{array}{rccll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& \begin{bmatrix}1 &\theta \end{bmatrix} \begin{bmatrix} X_{1,t} \\[3pt] X_{1,t-1}\end{bmatrix} + V_t \\[12pt] \\ \textrm{STATE EQUATION} \qquad\qquad & \begin{bmatrix} X_{1,t+1} \\[3pt] X_{1,t} \end{bmatrix} &=& \begin{bmatrix} 0 & 0 \\[3pt] 1 & 0 \end{bmatrix} \begin{bmatrix} X_{1,t} \\[3pt] X_{1,t-1}\end{bmatrix} + \begin{bmatrix} Z_{t+1} \\[3pt] 0 \end{bmatrix} \end{array} \] \textbf{EXAMPLE:} In the following simulation $\sigma_Z=1$, with $\bX_1 = (0,0,0,0)^\top$, for the $\ARMA(2,3)$ model given by \[ (1-0.8\B + 0.15 \B^2)X_t = (1+\B + \B^2 + \B^3) Z_t \] We have $m=q+1=4$, \[ \begin{array}{rccll} \textrm{OBSERVATION EQUATION} \qquad\qquad& Y_t &=& \begin{bmatrix} 1 & 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} X_{1,t-3} \\[3pt] X_{1,t-2} \\[3pt] X_{1,t-1} \\[3pt] X_{1,t} \end{bmatrix} + V_t \\[12pt] \\ \textrm{STATE EQUATION} \qquad\qquad & \begin{bmatrix} X_{1,t-2} \\[3pt] X_{1,t-1} \\[3pt] X_{1,t} \\[3pt] X_{1,t+1} \end{bmatrix} &=& \begin{bmatrix} 0 & 1 & 0 & 0\\[3pt] 0 & 0 & 1 & 0\\[3pt] 0 & 0 & 0 & 1\\[3pt] 0 & 0 & -0.15 & 0.8 \end{bmatrix} \begin{bmatrix} X_{1,t-3} \\[3pt] X_{1,t-2} \\[3pt] X_{1,t-1} \\[3pt] X_{1,t} \end{bmatrix} + \begin{bmatrix}0 \\[3pt] 0 \\[3pt] 0 \\[3pt] Z_{t+1} \end{bmatrix} \end{array} \] <>= set.seed(3834);n<-10000;p<-2;q<-3;m<-max(p,q+1) phi.vec<-theta.vec<-rep(0,m) phi.vec[1:p]<-c(0.8,-0.15) theta.vec[1:q]<-c(1,1,1) 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') @ \pagebreak To implement the Kalman Filter, we specify the various components in the usual code. <>= ##Compute Estimates using optim and a likelihoood computing function ml.armapq<-function(x,nv,yv,pv,qv){ m<-max(pv,qv+1);Kval<-1;Lval<-m phiv<-thetav<-rep(0,m) phiv[1:pv]<-x[1:pv] thetav[1:qv]<-x[(pv+1):(pv+qv)] Fm<-c(rev(thetav[1:(m-1)]),1) Gm<-matrix(0,m,m);for(j in 1:(m-1)){Gm[j,j+1]<-1} Gm[m,]<-rev(phiv) Sigmat<-0 Omegat<-diag(c(rep(0,m-1),exp(x[pv+qv+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,]<-rep(0,m) Ppt[1,,]<-Omegat #Filter Forward loglike<-0 for(i in 1:nv){ Fmat<-matrix(Fm,nrow=Kval,ncol=Lval) rt[i,]<-yv[i] - Fmat %*% apt[i,] Mt[i,,]<-Fmat %*% Ppt[i,,] %*% t(Fmat) + Sigmat if(Sigmat == Inf || Mt[i,,] == 0 || is.na(Mt[i,,])){return(Inf)} if(Mt[i,,]==Inf){ Mtinv<-matrix(0,1,1) }else{ Mtinv<-solve(Mt[i,,]) } att[i,]<-apt[i,] + Ppt[i,,] %*% t(Fmat)%*% Mtinv%*% rt[i,] Ptt[i,,]<-Ppt[i,,]-Ppt[i,,] %*% t(Fmat)%*%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)) } @ <>= #Perform the optimization using optim Y<-Y[1:1000] xv<-c(0.1,rep(0,p-1),rep(0,q+1)) ARMA.ml<-optim(xv,ml.armapq,nv=length(Y),pv=p,qv=q,yv=Y,method='BFGS',hessian=T, control=list(trace=10,fnscale=-1,maxit=500,abstol=1e-6)) @ After the optimization, the parameter estimates and standard errors are as follows: <>= phi.hat<-ARMA.ml$par[1:p] theta.hat<-ARMA.ml$par[(p+1):(p+q)] ese<-sqrt(diag(solve(-ARMA.ml$hessian))) phi.ese<-ese[1:p];theta.ese<-ese[(p+1):(p+q)] cbind(c(phi.hat,theta.hat),c(phi.ese,theta.ese)) (sigZ.hat<-exp(ARMA.ml$par[p+q+1])) @ These results resemble closely the results from \texttt{arima}: <>= arima(Y,order=c(2,0,3),method='ML',include.mean=FALSE) @ \end{document}