\documentclass[notitlepage]{article} \usepackage{../../Math545} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{bbm} \usepackage{bm} \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}} \newcommand{\bs}[1]{\bm{#1}} \newcommand{\Rho}{\mathrm{P}} \newcommand{\cif}{\text{ if }} \newcommand{\Ra}{\Longrightarrow} \newcommand{\ra}{\longrightarrow} \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{ARMA processes} \begin{document} \maketitle An $\ARMA(p,q)$ process $\{X_t\}$ is defined for $\{Z_t\} \sim \WN(0,\sigma_Z^2)$ and integers $p,q \geq 0$ as the process defined implicitly by \begin{equation} X_t = \sum_{j=1}^p \phi_j X_{t-j} + \sum_{j=0}^q \theta_{j} Z_{t-j} \qquad \forall t \label{eq:ARMApq} \end{equation} with $\theta_0 \equiv 1$, or equivalently \begin{equation} \Phi(B) X_t = \Theta(B) Z_t \label{eq:ARMAdef} \end{equation} where \[ \Phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p \qquad \qquad \Theta(B) = 1 + \theta_1 B + \theta B^2 + \cdots + \theta_q B^q. \] Using the factorization method, we write \[ \Phi(B) = \prod_{j=1}^p (1-\xi_j B) \qquad \qquad \Theta(B) = \prod_{j=1}^q (1-\omega_j B) \] and assume that there are no common factors in the two products (otherwise we could cancel terms on both sides, and reduce the orders of the two polynomials). We have shown that the process is stationary and causal if and only if $|\xi_j| < 1, j=1,\ldots,p$, in which case we know that the linear process representation \[ X_t = \Psi(B) Z_t = \left(\sum_{j=0}^\infty \psi_j B^j \right) Z_t = \sum_{j=0}^\infty \psi_j Z_{t-j} \] exists with $\sum\limits_{j=0}^\infty |\psi_j| < \infty$. \medskip We have also seen some general methods of computation. To identify the $\psi_j$, we can plug the proposed solution into \eqref{eq:ARMAdef} and equate coefficients: \[ \Phi(B) \Psi(B) = \Theta(B) \] or \[ (1 + \phi_1 B + \cdots + \phi_p B^p) (\psi_0 + \psi_1 B + \psi_2 B^2 + \cdots) = 1 + \theta_1 B + \theta_2 B^2 + \cdots + \theta_q B^q. \] Equating the coefficient of $B^j$, we have \[ \psi_j - \phi_1 \psi_{j-1} - \cdots - \psi_{j-p} \phi_p = \theta_j \qquad j = 0,1,2,\ldots \] where $\theta_j = 0$, $j > q$, and $\psi_j = 0$, $j < 0$. Rearranging, we get the recursive formula \[ \psi_j = \theta_j + \sum_{k=1}^{p} \phi_k \psi_{j-k} \qquad j = 1,2,\ldots \] with $\psi_0 = 1$, or more precisely \[ \psi_j = \left\{ \begin{array}{rll} \theta_j + \sum\limits_{k=1}^{j} \phi_k \psi_{j-k} & 0 \leq j < \max\{p,q+1\}\\[12pt] \sum\limits_{k=1}^{p} \phi_k \psi_{j-k} & j \geq \max\{p,q+1\} \end{array} \right. \] Note that this is merely a system of linear recurrence relations for which we have seen methods for constructing the exact solution using eigen decompositions. \pagebreak To illustrate this recursion, consider the $\ARMA(3,1)$ model \[ X_t = 0.5 X_{t-1} + 0.2 X_{t-2} - 0.2 X_{t-3} + Z_t + 2 Z_{t-1} \] <>= p<-3;q<-1 lmax<-200 phi.vec<-theta.vec<-rep(0,lmax+1) phi.vec[1:p]<-c(0.5,0.2,-0.2) theta.vec[1:q]<-c(2) true.psi.vec<-ARMAtoMA(ar=phi.vec[1:p],ma=theta.vec[1:q],lag.max=lmax) #True values psi.vec<-rep(0,lmax+1) psi.vec[0+1]<-1 for(j in 1:lmax){ kvec<-1:min(j,p) psi.vec[j+1]<-theta.vec[j]+sum(phi.vec[kvec]*psi.vec[j-kvec+1]) } cbind(1:10,psi.vec[1:10+1],true.psi.vec[1:10]) @ The other key computation is for the autocovariance sequence. The technique here is to multiply both sides of \eqref{eq:ARMApq} by $X_{t-h}$, and then take expectations: \begin{equation} \E[X_t X_{t-h}] = \sum_{j=1}^p \phi_j \E[X_{t-j} X_{t-h}] + \sum_{j=0}^q \theta_{j} \E[Z_{t-j} X_{t-h}]. \label{eq:ARMAacf1} \end{equation} Now \[ \E[Z_{t-j} X_{t-h}] = \sum_{k=0}^\infty \psi_k \E[Z_{t-j} Z_{t-h-k}] = \sigma_Z^2 \psi_{j-h} \] Thus we have for $h \geq 0$ \begin{equation} \gamma_X(h) - \sum_{j=1}^p \phi_j \gamma_X(h-j) = \sigma_Z^2 \sum_{j=h}^q \theta_{j} \psi_{j-h} \equiv \sigma_Z^2 \sum_{j=0}^{q-h} \theta_{j+h} \psi_{j} \label{eq:ARMAacf2} \end{equation} as, recall, $\psi_{j-h} = 0$ if $j < h$. From \eqref{eq:ARMAacf2}, we can obtain a linear system of equations, and then deduce the remaining values in the sequence from the fact that \begin{equation} \gamma_X(h)= \sum_{j=1}^p \phi_j \gamma_X(h-j) \label{eq:ARMAacf3}. \end{equation} whenever $h \geq m = \max\{p,q+1\}$. A rule of thumb is \begin{itemize} \item if $p \geq q+1$, inspect equations for lags $h=0,1,\ldots,p$; \item if $p < q+1$, inspect equations for lags $h=0,1,\ldots,q$; \end{itemize} \pagebreak \textbf{EXAMPLE:} For the $\ARMA(3,1)$ model \[ X_t = 0.5 X_{t-1} + 0.2 X_{t-2} - 0.2 X_{t-3} + Z_t + 2 Z_{t-1} \] we have $m=\max\{p,q+1\} = 3$ with $p > q+1$, and therefore inspect $h=0,1,2,3$. We have the system \[ \begin{array}{lccl} h = 0: \qquad & \gamma_X(0) & = & \phi_1 \gamma_X(1) + \phi_2 \gamma_X(2) + \phi_3 \gamma_X(3) + \sigma_Z^2 (1 + \theta_1 \psi_1)\\[6pt] h = 1: & \gamma_X(1) & = & \phi_1 \gamma_X(0) + \phi_2 \gamma_X(1) + \phi_3 \gamma_X(2) + \sigma_Z^2 \theta_1\\[6pt] h = 2: & \gamma_X(2) & = & \phi_1 \gamma_X(1) + \phi_2 \gamma_X(0) + \phi_3 \gamma_X(1) \\[6pt] h = 3: & \gamma_X(3) & = & \phi_1 \gamma_X(2) + \phi_2 \gamma_X(1) + \phi_3 \gamma_X(0) \\[6pt] \end{array} \] with $\phi_1 = 0.5$, $\phi_1=0.2$, $\phi_3 = -0.2$ and $\theta_1 = 2$. We re-write this as \[ \begin{bmatrix} 1 & -\phi_1 & - \phi_2 & -\phi_3 \\[6pt] - \phi_1 & 1-\phi_2 & - \phi_3 & 0 \\[6pt] - \phi_2 & -\phi_1-\phi_3 & 1 & 0 \\[6pt] - \phi_3 & -\phi_2 & -\phi_1 & 1 \end{bmatrix} \begin{bmatrix} \gamma_X(0) \\[6pt] \gamma_X(1) \\[6pt] \gamma_X(2) \\[6pt] \gamma_X(3) \end{bmatrix} = \begin{bmatrix}\sigma_Z^2 (1 + \theta_1 \psi_1) \\[6pt] \sigma_Z^2 \theta_1 \\[6pt] 0 \\[6pt] 0 \end{bmatrix} \] or $\bA \mathbf{x} = \mathbf{b}$, which we want to solve for $\mathbf{x}$. <>= p<-3;q<-1;m<-max(p,q+1) phi.vec<-theta.vec<-rep(0,lmax+1) phi.vec[1:p]<-c(0.5,0.2,-0.2) theta.vec[1:q]<-c(2) true.gammaX<-ARMAacf(ar=phi.vec[1:p],ma=theta.vec[1:q],lag.max=lmax) #True values Amat<-matrix(0,4,4) Amat[1,]<-c(1,-phi.vec[1:p]) Amat[2,]<-c(-phi.vec[1],1-phi.vec[2],-phi.vec[3],0) Amat[3,]<-c(-phi.vec[2],-phi.vec[1]-phi.vec[3],1,0) Amat[4,]<-c(-phi.vec[3],-phi.vec[2],-phi.vec[1],1) bvec<-c(1+theta.vec[1]*psi.vec[1+1], theta.vec[1],0,0) gammaX<-solve(Amat,bvec) gammaX/gammaX[1] true.gammaX[1:4] @ Then using \eqref{eq:ARMAacf3}, we can compute the remaining values (out to $h=9$ say): <>= for(h in (m+1):lmax){ gammaX<-c(gammaX,sum(phi.vec*gammaX[h-(1:p)+1])) } cbind(gammaX[1:10]/gammaX[1],true.gammaX[1:10]) #Comparison of equivalent computations @ \textbf{EXAMPLE:} For the $\ARMA(2,4)$ model \[ X_t = X_{t-1} - 0.25 X_{t-2} + Z_t + Z_{t-1} + 3 Z_{t-2} - Z_{t-3} + 2 Z_{t-4} \] we have $m=\max\{p,q+1\} = 5$ with $p < q+1$, and therefore inspect $h=0,1,2,3,4$. We have the system \[ \begin{array}{lccl} h = 0: \qquad & \gamma_X(0) & = & \phi_1 \gamma_X(1) + \phi_2 \gamma_X(2) + \sigma_Z^2 (1 + \theta_1 \psi_1+ \theta_2 \psi_2 + \theta_3 \psi_3+ \theta_4 \psi_4)\\[6pt] h = 1: & \gamma_X(1) & = & \phi_1 \gamma_X(0) + \phi_2 \gamma_X(1) + \sigma_Z^2 (\theta_1 + \theta_2 \psi_1 + \theta_3 \psi_2+ \theta_4 \psi_3)\\[6pt] h = 2: & \gamma_X(2) & = & \phi_1 \gamma_X(1) + \phi_2 \gamma_X(0) + \sigma_Z^2 (\theta_2 + \theta_3 \psi_1 + \theta_4 \psi_2) \\[6pt] h = 3: & \gamma_X(3) & = & \phi_1 \gamma_X(2) + \phi_2 \gamma_X(1) + \sigma_Z^2 (\theta_3 + \theta_4 \psi_1 )\\[6pt] h = 4: & \gamma_X(4) & = & \phi_1 \gamma_X(3) + \phi_2 \gamma_X(2) + \sigma_Z^2 \theta_4 \\[6pt] \end{array} \] with $\phi_1 = 1$, $\phi_2 = -0.25$, and $\theta_1 = 1, \theta_2 = 3, \theta_3 = -1, \theta_4 = 2$. We re-write this as \[ \begin{bmatrix} 1 & -\phi_1 & - \phi_2 & 0 & 0 \\[6pt] - \phi_1 & 1-\phi_2 & 0 & 0 & 0\\[6pt] - \phi_2 & -\phi_1 & 1 & 0 & 0\\[6pt] 0 & -\phi_2 & -\phi_1 & 1 & 0\\[6pt] 0 & 0 & -\phi_2 & -\phi_1 & 1 \end{bmatrix} \begin{bmatrix} \gamma_X(0) \\[6pt] \gamma_X(1) \\[6pt] \gamma_X(2) \\[6pt] \gamma_X(3) \\[6pt] \gamma_X(4) \end{bmatrix} = \begin{bmatrix}\sigma_Z^2 (1 + \theta_1 \psi_1+ \theta_2 \psi_2 + \theta_3 \psi_3+ \theta_4 \psi_4) \\[6pt] \sigma_Z^2 (\theta_1 + \theta_2 \psi_1 + \theta_3 \psi_2+ \theta_4 \psi_3) \\[6pt] \sigma_Z^2 (\theta_2 + \theta_3 \psi_1 + \theta_4 \psi_2) \\[6pt] \sigma_Z^2 (\theta_3 + \theta_4 \psi_1 ) \\[6pt] \sigma_Z^2 \theta_4 \end{bmatrix} \] or $\bA \mathbf{x} = \mathbf{b}$, which we want to solve for $\mathbf{x}$. <>= p<-2;q<-4;m<-max(p,q+1) phi.vec<-theta.vec<-rep(0,lmax+1) phi.vec[1:p]<-c(1,-0.25) theta.vec[1:q]<-c(1,3,-1,2) true.psi.vec<-ARMAtoMA(ar=phi.vec[1:p],ma=theta.vec[1:q],lag.max=lmax) #True values psi.vec<-rep(0,lmax+1) psi.vec[0+1]<-1 for(j in 1:lmax){ kvec<-1:min(j,p) psi.vec[j+1]<-theta.vec[j]+sum(phi.vec[kvec]*psi.vec[j-kvec+1]) } cbind(1:10,psi.vec[1:10+1],true.psi.vec[1:10]) #Comparison of equivalent computations true.gammaX<-ARMAacf(ar=phi.vec[1:p],ma=theta.vec[1:q],lag.max=lmax) #True values Amat<-matrix(0,5,5) Amat[1,]<-c(1,-phi.vec[1:p],0,0) Amat[2,]<-c(-phi.vec[1],1-phi.vec[2],0,0,0) Amat[3,]<-c(-phi.vec[2],-phi.vec[1],1,0,0) Amat[4,]<-c(0,-phi.vec[2],-phi.vec[1],1,0) Amat[5,]<-c(0,0,-phi.vec[2],-phi.vec[1],1) bvec<-c(sum(c(1,theta.vec[1:4])*psi.vec[1:5]),sum(theta.vec[1:4]*psi.vec[1:4]), sum(theta.vec[2:4]*psi.vec[1:3]),sum(theta.vec[3:4]*psi.vec[1:2]), sum(theta.vec[4]*psi.vec[1])) gammaX<-solve(Amat,bvec) gammaX/gammaX[1] true.gammaX[1:5] @ Then using \eqref{eq:ARMAacf3}, we can compute the remaining values (out to $h=9$ say): <>= for(h in m:lmax){ gammaX<-c(gammaX,sum(phi.vec[1:p]*gammaX[h-(1:p)+1])) } cbind(gammaX[1:10]/gammaX[1],true.gammaX[1:10]) #Comparison of equivalent computations @ For a simulation from this process: <>= set.seed(234) X<-arima.sim(n=200,model=list(ar=phi.vec[1:p],ma=theta.vec[1:q])) par(mar=c(4,4,0,1)) plot(X,type='l',ylim=range(-30,30));points(X,pch=19,cex=0.6) @ \end{document}