\documentclass[notitlepage]{article}
\usepackage{../../Math545}
\usepackage{listings}
\usepackage{numprint}
\usepackage{enumerate}
\usepackage{bbm}

\usepackage{amsfonts,amsmath}

\setlength{\parindent}{0in}

\def\Var{\text{Var}}
\def\Cov{\text{Cov}}
\def\Corr{\text{Corr}}



\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{Stationary, Causal and Invertible \ARMA~solutions}


\begin{document}

\maketitle

$\{X_t\}$ is an $\ARMA(p,q)$ process if it is stationary and satisfies
\begin{equation}
\Phi(B) X_t = \Theta(B) Z_t \label{eq:ARMA}
\end{equation}
where $\{Z_t\} \sim \WN(0,\sigma_Z^2)$, and
\begin{eqnarray*}
\Phi(z) & = & 1 - \phi_1 z - \cdots - \phi_p z^p = \prod_{j=1}^p (1-\xi_j z)\\
\Theta(z) & = & 1 + \theta_1 z + \cdots + \theta_q z^q = \prod_{j=1}^q (1-\omega_j z)\\
\end{eqnarray*}
That is
\[
X_t - \phi_1 X_{t-1} - \phi_2 X_{t-2} - \cdots - \phi_p X_{t-p} = Z_t + \theta_1 Z_{t-1} + \theta_2 Z_{t-2} + \cdots \theta_q Z_{t-q}
\]
or equivalently
\[
X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \cdots + \phi_p X_{t-p} + Z_t + \theta_1 Z_{t-1} + \theta_2 Z_{t-2} + \cdots \theta_q Z_{t-q}
\]

\textbf{Reciprocal roots:}  For the AR component, the parameters $\xi_1,\ldots,\xi_p$ are the reciprocal roots of the equation
\[
\Phi(z) = 1 - \phi_1 z - \cdots - \phi_p z^p = 0
\]
that is, the reciprocals of the roots, $\eta_j$, of the equation
\[
\prod_{j=1}^p (\eta_j - z) = 0
\]
To see this, note that $(1-\xi z) = \xi (\xi^{-1} - z) = \xi (\eta - z)$ so therefore
\[
\prod_{j=1}^p (1-\xi_j z) = \left( \prod_{j=1}^p \xi_j \right) \prod_{j=1}^p (\xi_j^{-1} - z).
\]
The function \texttt{polyroot} in \texttt{R} can be used to find the reciprocal roots from an arbitrary polynomial.

\medskip
\textbf{Example:} Consider the polynomial
\[
\Phi(z) = 1 - 2 z + \frac{1}{2} z^2 + \frac{2}{3} z^3
\]
to define the $\AR(3)$ model
\[
X_t - 2 X_{t-1} + \frac{1}{2} X_{t-2} + \frac{2}{3} X_{t-3} = Z_t
\]
or equivalently
\[
X_t = 2 X_{t-1} - \frac{1}{2} X_{t-2} - \frac{2}{3} X_{t-3} + Z_t.
\]

<<P1>>=
covec<-c(1,-2,1/2,2/3)  #Enter the coefficient of 1, z, z^2, z^3
polyroot(covec)
1/polyroot(covec)
@
The three reciprocal roots are
\[
\Sexpr{round(1/polyroot(covec),4)}
\]
that is, there is one real reciprocal root, \Sexpr{Re(round(1/polyroot(covec),4)[3])}, and two complex reciprocal roots, that appear in a complex conjugate pair.  We can check the modulus of the three roots using the function \texttt{Mod}:
<<P2>>=
Mod(1/polyroot(covec))
@
This tells us that the complex-valued reciprocal roots have modulus \Sexpr{Mod(1/polyroot(covec))[1]}, and so no stationary AR(3) component can be defined with these coefficients.  We can example the AR(3) process defined in this way.
<<P3,fig.height=4,fig.width=8>>=
set.seed(23)
n<-20
Z<-rnorm(n)
X<-rep(0,n)
for(i in 4:n){
    X[i]<-2*X[i-1]-(1/2)*X[i-2]-(2/3)*X[i-3]+Z[i]
}
par(mar=c(4,4,2,2))
plot(X,type='l',xlab='t')
points(X,pch=19,cex=0.6)
@

The system exhibits explosive behaviour.

\medskip
We can construct a stationary process by finding the coefficients corresponding to a particular set of reciprocal roots.  For example, if we consider entirely real reciprocal roots
\[
\xi_1  = -0.1, \quad \xi_2 = 0.25, \quad \xi_3 = 0.75
\qquad \Longrightarrow \qquad
\eta_1 = -10, \quad \eta_2 = 4, \quad \eta_3 = \frac{4}{3}
\]
so that the $\AR(3)$ polynomial
\[
(z-\eta_1) (z-\eta_2) (z-\eta_3)
\]
can be computed directly using the \texttt{polynom} package.  The \texttt{polynom} package allows multiplication of polynomials, in this case, to perform the calculation
\[
(10+z)(-4+z)(-4/3+z)
\]
<<c4>>=
library(polynom)
t1 <- polynomial(c(10,1))    #First term:  10+z
t2 <- polynomial(c(-4,1))    #Second term: -4+z
t3 <- polynomial(c(-4/3,1))  #Third term: -4/3+z
(t1*t2)*t3
avec<-coefficients((t1*t2)*t3)
polyroot(avec)               #Roots
polyroot(rev(avec))          #Reciprocal roots
covec<-avec/avec[1]
covec
@
Thus the $\AR(3)$ model
\[
X_t = \Sexpr{covec[2]} X_{t-1} + \Sexpr{covec[3]} X_{t-2} + \Sexpr{covec[4]} X_{t-3} + Z_t
\]
should be stationary.  We use the \texttt{arima.sim} to check this:
<<c5,fig.height=4.5,fig.width=7>>=
n<-500
X<-arima.sim(n,model=list(ar=covec[2:4]))
par(mar=c(4,4,1,2))
plot(X,type='l'); points(X,pch=19,cex=0.5)
acf(X,lag.max=50)
ARMAacf(ar=covec[2:4],lag.max=10)
@
The process is stationary, and has ACF alternating in sign, with odd lags exhibiting negative correlation, and even lags exhibiting positive dependence, both decaying to negligible values before lag 25.

\bigskip

\textbf{Invertibility:}  In the $\ARMA(p,q)$ process in \eqref{eq:ARMA}, the existence of a stationary and causal solution depends only on the AR polynomial $\Phi(z)$.  However, it is sometimes useful to be able to re-write \eqref{eq:ARMA} in the form
\[
\{\Theta(B)\}^{-1} \Phi(B) X_t =  Z_t
\]
filtering both sides of \eqref{eq:ARMA} by the inverse operator filter $\{\Theta(B)\}^{-1}$ to give the series expansion
\begin{equation}
Z_t = \pi(B) X_t = \sum_{j=0}^\infty \pi_j X_{t-j} \label{eq:ARInfinity}
\end{equation}
where
\[
\pi(B) = \{\Theta(B)\}^{-1} \Phi(B)
\]
has the series expansion
\[
\pi(B) =  \sum_{j=0}^\infty \pi_j B^j.
\]
If such a convergent expansion given in \eqref{eq:ARInfinity} exists, then $\{X_t\}$ is termed an \textit{invertible} process.  By identical arguments to those made for the AR model, the model is an invertible if
\[
\sum_{j=0}^\infty |\pi_j| < \infty
\]
and in the ARMA case, if the AR component is stationary and causal, the model is invertible if
\[
\Theta(z) = 1 + \theta_1 z + \theta_2 z^2 + \cdots + \theta_q z^2 = \prod_{j=1}^q (1-\omega_j z)
\]
gives
\[
\beta(z) = \{\Theta(z)\}^{-1} = \prod_{j=1}^q (1-\omega_j z)^{-1}
\]
which has a convergent series expansion when $|z| = 1$.  This is equivalent to requiring that $|\omega_j| < 1$ for all $j$.

\bigskip

\textbf{Example:}  Consider the MA(2) process
\[
X_t = Z_t - 5 Z_{t-1} + 6 Z_{t-2} = (1-2B)(1-3B) Z_t
\]
This process is not invertible, as the expansions of $(1-2z)^{-1}$ and $(1-3z)^{-1}$ are divergent when $|z| = 1$.  Despite this, the process is stationary, and well behaved.  If we take $\{Z_t\}$ to be an IID Gaussian sequence with $\sigma_Z^2=1$, the process can be simulated.

<<c6,fig.height=4.5,fig.width=7>>=
Z<-rnorm(502)
X<-Z[3:502]-5*Z[2:501]+6*Z[1:500]
par(mar=c(4,4,2,2))
plot(X,type='l',xlab='Time')
points(X,pch=19,cex=0.5)
@
\bigskip

Under stationarity, with a causal AR component and an invertible MA component, we have that $\pi(B)$ has generating function
\begin{equation}
\pi(z) = \frac{\Phi(z)}{\Theta(z)} = \dfrac{\prod\limits_{j=1}^p (1-\xi_j z)}{\prod\limits_{k=1}^q (1-\omega_k z)} \label{eq:ARMAgf}
\end{equation}
whereas the linear process representation via $\Psi(B)$ with generating function
\begin{equation}
\Psi(z) = \frac{\Theta(z)}{\Phi(z)} = \dfrac{\prod\limits_{k=1}^q (1-\omega_k z)}{\prod\limits_{j=1}^p (1-\xi_j z)} \label{eq:ARMAgf}
\end{equation}
where $|\xi_j|< 1$ for all $j$ and $|\omega_k| < 1$ for all $k$.


\bigskip
\textbf{Non-causal case:} In the AR(1) model, if $|\phi| > 1$, it is possible to find a stationary but non-causal solution to the AR equation: recall that if $|\phi|>1$
\[
X_{t} = \phi X_{t-1} + Z_t
\]
with $\{Z_t\} \sim \WN(0,\sigma_Z^2)$, then on rearrangement we have
\[
X_{t-1} = \frac{1}{\phi} X_t - \frac{1}{\phi} Z_t
\qquad \text{or} \qquad
X_{t-1} = \phi_\ast X_t + Z_t^\ast
\]
say, where $\phi_\ast = 1/\phi$ and where $Z_t^\ast = - Z_t/\phi$ so that $\{Z_t^\ast\} \sim \WN(0,\sigma_Z^2/\phi^2)$.   Using the standard recursion approach
\begin{align*}
  X_t = \phi_\ast X_{t+1} + Z_{t+1}^\ast & = \phi_\ast^2 X_{t+2} + (\phi_\ast Z_{t+2}^\ast + Z_{t+1}^\ast)\\[6pt]
    & = \phi_\ast^3 X_{t+3} + (\phi_\ast^2 Z_{t+3}^\ast + \phi_\ast Z_{t+2}^\ast + Z_{t+1}^\ast)
\end{align*}
etc, which yields the solution
\begin{equation}
X_t = \sum_{j=0}^\infty \phi_\ast^j Z_{t+1+j}^\ast = -\sum_{j=0}^\infty \phi^{-(j+1)} Z_{t+1+j} \label{eq:NCAR}
\end{equation}

\textbf{Example:}  Consider the AR(1) process
\[
X_t = 1.5 X_{t-1} + Z_t
\]
This process is not stationary if we restrict attention to causal solutions.  However, the stationary non-causal solution based on \eqref{eq:NCAR} can be computed.  If we take $\{Z_t\}$ to be an IID Gaussian sequence with $\sigma_Z^2=1$, the process can be simulated.  The variance of this process satisfies
\[
\Var[X_{t-1}] = \phi_\ast^2 \Var[X_{t}] + \Var[Z_t^\ast] \quad \Longrightarrow \quad \Var[X_t] = \frac{\Var[Z_t^\ast]}{(1-\phi_\ast^2)} = \frac{1/\phi_\ast^2}{(1-\phi_\ast^2)}=\frac{1}{\phi^2-1}.
\]
and as $\phi_\ast = 1/1.5 = \Sexpr{round(1/1.5,4)}$, we have that $\Var[X_t]=\Sexpr{((1/1.5)^2)/(1-1/(1.5^2))}$.

<<c7,fig.height=4.0,fig.width=7>>=
phi<-1.5; phistar<-1/phi
Z<-rnorm(1001); Zstar<--Z/phi;
psivec<-phistar^c(0:500)
for(i in 1:500){
    X[i]<-sum(psivec*Zstar[(i+1):(i+1+500)])  #Generate X using (5)
}
var(X)                   #Variance
sdX<-sqrt(1/(phi^2-1))
X[500]<-rnorm(1,0,sdX)                        #Generate X using the recursion
for(i in 499:1){
    X[i]<-phistar*X[i+1]+Zstar[i+1]
}
var(X)
par(mar=c(4,4,1,2))
plot(X,type='l',xlab='Time')
points(X,pch=19,cex=0.5)
acf(X)
@

We can compute the ACVF for $h \geq 0$ as
\begin{align*}
\E[X_t X_{t+h}] & = \E\left[ \left(\sum_{j=0}^\infty \phi_\ast^j Z_{t+1+j}^\ast \right)  \left(\sum_{k=0}^\infty \phi_\ast^k Z_{t+h+1+k}^\ast \right) \right]\\[6pt]
& = \sum_{j=0}^\infty \sum_{k=0}^\infty \phi_\ast^{j+k} \E\left[ Z_{t+1+j}^\ast Z_{t+h+1+k}^\ast \right]
\end{align*}
The expectation is only non-zero when the indices match and $(t+1+j) = (t+h+1+k)$, that is, when $k=j-h$, in which case the expectation is $\Var[Z_t^\ast] = \sigma_Z^2/\phi^2$.  Therefore
\[
\E[X_t X_{t+h}] = \frac{\sigma_Z^2}{\phi^2} \sum_{j=h}^\infty \phi_\ast^{2j - h} = \frac{\sigma_Z^2}{\phi^2} \sum_{j=0}^\infty \phi_\ast^{2(j+h)-h} = \frac{\sigma_Z^2}{\phi^2} \phi_\ast^h \sum_{j=0}^\infty \phi_\ast^{2j} = \frac{\sigma_Z^2}{\phi^2} \phi_\ast^h \frac{1}{1-\phi_\ast^2}.
\]
as $\phi_\ast = 1/\phi$, after summing the geometric progression. Simplifying, we obtain for the stationary but non-causal solution that
\[
\gamma_X(h) = \frac{\sigma_Z^2}{\phi^2-1} \frac{1}{\phi^h}  \qquad h=0,1,2,\ldots.
\]

<<c8,fig.height=4.0,fig.width=7>>=
h<-1:10
gamvec<-1/(phi^2-1)*(1/phi)^h
gamvec
rhovec<-gamvec
rhovec/rhovec[1]
acf(X,lag.max=10,plot=F)
@


\end{document}

\medskip
\textbf{Example:} Consider the $\ARMA(1,1)$ model
\begin{equation}
(1-\phi B) X_t = (1+\theta B) Z_t \label{eq:ARMAZ}
\end{equation}
where $|\phi|>1$ and $|\theta| > 1$, so that the process is neither causal nor invertible.  We can construct a solution that is causal and invertible with respect to a different white noise process in the following way.  First, it is evident that the $\ARMA(1,1)$ process defined as the solution to
\begin{equation}
(1-\phi^{-1} B) X_t = (1+\theta^{-1} B) Z_t^\ast, \label{eq:ARMAZstar}
\end{equation}
where $\{Z_t^\ast\}$ is another white noise process, is stationary, causal and invertible.


\medskip
As shown previously, the linear process representation for a stationary and causal $\ARMA(1,1)$ model is determined by the coefficients
\begin{eqnarray*}
\psi_0 & = & 1\\
\psi_j & = & (\theta + \phi) \phi^{j-1}, \; j \geq 1
\end{eqnarray*}
and the autocovariance function $\gamma_X(h)$ is given by
\[
\gamma_X(h) = \sigma_Z^2 \sum_{j=0}^\infty \psi_j \psi_{j+h}.
\]
Therefore for the process defined by \eqref{eq:ARMAZstar}, we have that for $j \geq 1$, the linear process coefficients are
\[
\psi_j^\ast = (\theta^{-1} + \phi^{-1}) \phi^{-(j-1)}
\]
Now for $h \geq 2$
\begin{align*}
  \sum_{j=0}^\infty \psi_j \psi_{j+h} & = (\theta + \phi) \phi^{h-1} + \sum_{j=1}^\infty (\theta+\phi)^2 \phi^{2 (j-1) + h}\\[6pt]
  & = (\theta + \phi) \phi^{h-1} + (\theta+\phi)^2 \phi^h \sum_{j=0}^\infty \phi^{2j}\\[6pt]
  & = (\theta + \phi) \phi^{h-1} + (\theta+\phi)^2 \phi^h \frac{1}{1-\phi^2} & \text{geometric progression}\\[6pt]
  & = \left[ (\theta/\phi+1) + \frac{(\theta+\phi)^2}{1-\phi^2} \right] \phi^h
\end{align*}
Similarly
\[
\sum_{j=0}^\infty \psi_j^\ast \psi_{j+h}^\ast = (\theta^{-1} + \phi^{-1}) \phi^{1-h} + (\theta^{-1}+\phi^{-1})^2 \phi^{-h} \frac{1}{1-\phi^{-2}}
\]
But recall that the ACVF must be an even function of $h$, so we must have that $\gamma_X^\ast(h) = \gamma_X^\ast (-h)$ for $h<0$.  Note also that
\[
(\theta^{-1} + \phi^{-1}) \phi = (\phi/\theta+1)
\]
and
\[
(\theta^{-1}+\phi^{-1})^2 \frac{1}{1-\phi^{-2}} = (\phi/\theta+1)^2 \frac{\phi^2}{\phi^2-1}
\]