20  ARMA processes

library(dplyr)
# required for figures  
library(ggplot2)
library(gridExtra)
# required for render latex 
library(backports)
library(latex2exp)
# required for render tables
library(knitr)
library(kableExtra)
# Random seed 
set.seed(1)

20.1 ARMA(p, q)

An Auto Regressive Moving Average processes, or for short ARMA(p,q), is a combination of an MA(q) (Equation 19.1) and an AR(p) (Equation 19.3), i.e.  $$\begin{array}{}\text{(20.1)}& Y_t=\mu +\sum _{i=1}^p{\varphi }_i{Y}_{t-i}+\sum _{j=1}^q{\theta }_j{u}_{t-j}+u_t\text{,}\end{array}$$ where in general $u_t\sim \text{WN}(0,{\sigma }_{\text{u}}^2)$ (Definition 18.5).

set.seed(6) # random seed
# *************************************************
#                      Inputs 
# *************************************************
# Number of steps ahead 
t_bar <- 300 
# Intercept 
mu <- 0.5
# AR parameters
phi <- c(phi1 = 0.95)
# MA parameters
theta <- c(theta1 = 0.45)
# Variance of the residuals 
sigma2_u <- 1 
# *************************************************
#                   Simulations 
# *************************************************
# Initial value
Yt <- c(mu / (1 - phi[1])) 
# Simulated residuals
u_t <- rnorm(t_bar, 0, sqrt(sigma2_u)) 
# Simulated ARMA(1,1)
for(t in 2:t_bar){
  # AR component
  sum_AR <- phi[1] * Yt[t-1]
  # MA component
  sum_MA <- theta[1] * u_t[t-1] 
  # Update Yt 
  Yt[t] <- mu + sum_AR + sum_MA + u_t[t]
}

Figure 20.1: ARMA(1,1) simulation on the top and empirical autocovariance at the bottom.

20.1.1 Matrix form

An Autoregressive process of order p AR(p) (Equation 19.3) can be written in matrix form as $${\mathbf{Y}}_t=\mathbf{c}+\mathbf{\Phi }{\mathbf{Y}}_{t-1}+\mathbf{b}{u}_t\text{,}$$ where $$\underset{{\mathbf{Y}}_t}{\underbrace{\left(\begin{array}{c}{Y}_t\\ Y_{t-1}\\ ⋮\\ Y_{t-p}\end{array}\right)}}=\underset{\mathbf{c}}{\underbrace{\left(\begin{array}{c}\mu \\ 0\\ ⋮\\ 0\end{array}\right)}}+\underset{\mathbf{\Phi }}{\underbrace{\left(\begin{array}{ccccc}{\varphi }_1& {\varphi }_2& \dots & {\varphi }_{p-1}& {\varphi }_p\\ 1& 0& \dots & 0& 0\\ 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 1& 0\\ 0& 0& \dots & 0& 0\end{array}\right)}}\underset{{\mathbf{Y}}_{t-1}}{\underbrace{\left(\begin{array}{c}{Y}_{t-1}\\ Y_{t-2}\\ ⋮\\ Y_{t-p-1}\end{array}\right)}}+\underset{\mathbf{b}}{\underbrace{\left(\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right)}}\cdot u_t$$ where ${\mathbf{Y}}_t$, $\mathbf{c}$ and $\mathbf{b}$ have dimension $p\times 1$, while $\mathbf{\Phi }$ is $p\times p$.

How to construct the companion matrix $\mathbf{\Phi }$?

Let’s consider the vector containing the coefficients of the model, i.e. $$\underset{p\times 1}{\mathit{\gamma }}=\left(\begin{array}{ccccc}{\varphi }_1& {\varphi }_2& \dots & {\varphi }_{p-1}& {\varphi }_p\end{array}\right)\text{.}$$ If the order of the Autoregressive process is greater than 1, i.e. $p>1$, let’s consider an identity matrix (Equation 31.3) with dimension $(p-1)\times (p-1)$, i.e.  $${\mathbf{I}}_{p-1}=\left(\begin{array}{ccccc}1& 0& \dots & 0& 0\\ 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 1& 0\\ 0& 0& \dots & 0& 1\end{array}\right)\text{,}$$ and combine it with a column of zeros, i.e. $$\underset{(p-1)\times p}{\mathbf{L}}=\left(\begin{array}{cc}{\mathbf{I}}_{p-1}& {\mathbf{0}}_{1\times (p-1)}\end{array}\right)\text{.}$$ Finally, add on the top of $\mathbf{L}$ with the AR parameters, i.e.  $$\underset{p\times p}{\mathbf{\Phi }}=\left(\begin{array}{c}\mathit{\gamma }\\ \mathbf{L}\end{array}\right)\text{.}$$

Example: AR(4) in matrix form

Example 20.1 Let’s consider an AR(4) process of the form $$Y_t=\mu +{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+{\varphi }_3{Y}_{t-3}+{\varphi }_4{Y}_{t-4}+u_t\text{.}$$ Then, the vector containing the coefficients of the model reads $$\underset{4\times 1}{\mathit{\gamma }}=\left(\begin{array}{cccc}{\varphi }_1& {\varphi }_2& {\varphi }_3& {\varphi }_4\end{array}\right)\text{.}$$

For an AR(4) we start with a diagonal matrix with dimension $3\times 3$, i.e.  $${\mathbf{I}}_3=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\text{,}$$ and we add a column of zeros, i.e. $$\underset{3\times 4}{{\mathbf{L}}_3}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right)\text{.}$$ and the parameters $$\underset{4\times 4}{\mathbf{\Phi }}=\left(\begin{array}{cccc}{\varphi }_1& {\varphi }_2& {\varphi }_3& {\varphi }_4\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right)\text{.}$$ In order to verify the above result, let’s consider the iteration from $t-1$ to time $t$, in matrix form one obtain $$\left(\begin{array}{c}{Y}_t\\ Y_{t-1}\\ Y_{t-2}\\ Y_{t-3}\end{array}\right)=\left(\begin{array}{c}\mu \\ 0\\ 0\\ 0\end{array}\right)+\left(\begin{array}{cccc}{\varphi }_1& {\varphi }_2& {\varphi }_3& {\varphi }_4\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right)\left(\begin{array}{c}{Y}_{t-1}\\ Y_{t-2}\\ Y_{t-3}\\ Y_{t-4}\end{array}\right)+\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right)u_t\text{.}$$ After an explicit computation, one can verify that it lead to the classic AR(4) recursion, i.e. $$\begin{array}{rl}\left(\begin{array}{c}{Y}_t\\ Y_{t-1}\\ Y_{t-2}\\ Y_{t-3}\end{array}\right)& =\left(\begin{array}{c}\mu \\ 0\\ 0\\ 0\end{array}\right)+\left(\begin{array}{c}{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+{\varphi }_3{Y}_{t-3}+{\varphi }_4{Y}_{t-4}\\ Y_{t-1}\\ Y_{t-2}\\ Y_{t-3}\end{array}\right)+\left(\begin{array}{c}{u}_t\\ 0\\ 0\\ 0\end{array}\right)\end{array}$$

Let’s now consider the matrix form an ARMA (p,q) process (Equation 20.1), i.e.  $${\mathbf{X}}_t=\mathbf{c}+\mathbf{A}{\mathbf{X}}_{t-1}+\mathbf{b}{u}_t\text{,}$$ where, the first component, namely $Y_t={\mathbf{e}}_{p+q}^{\mathrm{\top }}{\mathbf{X}}_t$, is extracted as: $$Y_t={\mathbf{e}}_{p+q}^{\mathrm{\top }}\mathbf{c}+{\mathbf{e}}_{p+q}^{\mathrm{\top }}\mathbf{A}{\mathbf{X}}_{t-1}+{\mathbf{e}}_{p+q}^{\mathrm{\top }}\mathbf{b}{u}_t\text{.}$$ where ${\mathbf{e}}_{p+q}^{\mathrm{\top }}$ is a basis vector (Equation 31.1) with dimension $p+q$.

$$\underset{{\mathbf{X}}_t}{\underbrace{\left(\begin{array}{c}{Y}_t\\ Y_{t-1}\\ ⋮\\ Y_{t-p}\\ u_t\\ u_{t-1}\\ ⋮\\ u_{t-q}\end{array}\right)}}=\underset{\mathbf{c}}{\underbrace{\left(\begin{array}{c}\mu \\ 0\\ 0\\ ⋮\\ 0\\ 0\end{array}\right)}}+\underset{\mathbf{A}}{\underbrace{\left(\begin{array}{cccccccccc}{\varphi }_1& {\varphi }_2& \dots & {\varphi }_{p-1}& {\varphi }_p& {\theta }_1& {\theta }_2& \dots & {\theta }_{q-1}& {\theta }_q\\ 1& 0& \dots & 0& 0& 0& 0& \dots & 0& 0\\ 0& 1& \dots & 0& 0& 0& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 1& 0& 0& 0& \dots & 0& 0\\ 0& 0& \dots & 0& 0& 0& 0& \dots & 1& 0\end{array}\right)}}\underset{{\mathbf{X}}_{t-1}}{\underbrace{\left(\begin{array}{c}{Y}_{t-1}\\ Y_{t-2}\\ ⋮\\ Y_{t-p-1}\\ u_{t-1}\\ u_{t-2}\\ ⋮\\ u_{t-p-1}\end{array}\right)}}+\underset{\mathbf{b}}{\underbrace{\left(\begin{array}{c}1\\ 0\\ ⋮\\ 0\\ 1\\ 0\\ ⋮\\ 0\end{array}\right)}}\cdot u_t$$

How to construct the companion matrix for an ARMA?

Let’s consider the vector $\gamma$ containing the coefficients of the model, i.e. $$\underset{1\times (q+p)}{\mathit{\gamma }}=\left(\begin{array}{cccccccccc}{\varphi }_1& {\varphi }_2& \dots & {\varphi }_{p-1}& {\varphi }_p& {\theta }_1& {\theta }_2& \dots & {\theta }_{q-1}& {\theta }_q\end{array}\right)\text{.}$$ For an ARMA(p,q) we consider two distinct matrices, i.e. a matrix for the AR part we consider the identity matrix (Equation 31.3) with order $p-1$, i.e. $$\underset{(p-1)\times (p-1)}{{\mathbf{I}}_{p-1}}=\left(\begin{array}{ccccc}1& 0& \dots & 0& 0\\ 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 1& 0\\ 0& 0& \dots & 0& 1\end{array}\right)\text{,}$$ and we combine it with a column of zeros, i.e. $$\underset{(p-1)\times p}{{\mathbf{L}}^{\text{AR}}}=\left(\begin{array}{cc}{\mathbf{I}}_{p-1}& {\mathbf{0}}_{(p-1)\times 1}\end{array}\right)=\left(\begin{array}{cccccc}1& 0& \dots & 0& 0& 0\\ 0& 1& \dots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& 0\\ 0& 0& \dots & 1& 0& 0\\ 0& 0& \dots & 0& 1& 0\end{array}\right)\text{.}$$ For the MA part we consider the identity matrix with order $q-1$, i.e. $${\mathbf{I}}_{q-1}=\left(\begin{array}{ccccc}1& 0& \dots & 0& 0\\ 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 1& 0\\ 0& 0& \dots & 0& 1\end{array}\right)\text{,}$$ and we combine it with a column of zeros, i.e. $$\underset{(q-1)\times q}{{\mathbf{L}}^{\text{MA}}}=\left(\begin{array}{cc}{\mathbf{I}}_{q-1}& {\mathbf{0}}_{(q-1)\times 1}\end{array}\right)=\left(\begin{array}{cccccc}1& 0& \dots & 0& 0& 0\\ 0& 1& \dots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& 0\\ 0& 0& \dots & 1& 0& 0\\ 0& 0& \dots & 0& 1& 0\end{array}\right)\text{.}$$ To combine the matrices ${\mathbf{L}}_{p-1}$ and ${\mathbf{L}}_{q-1}$, we need add a matrix of zeros. More precisely: $$\underset{(q+p-1)\times (q+p)}{\mathbf{L}}=\left(\begin{array}{cc}{\mathbf{L}}^{\text{AR}}& {\mathbf{0}}_{(p-1)\times q}\\ {\mathbf{0}}_{1\times p}& {\mathbf{0}}_{1\times q}\\ {\mathbf{0}}_{(q-1)\times p}& {\mathbf{L}}^{\text{MA}}\end{array}\right)\text{.}$$ Finally we combine $\mathbf{L}$ with the parameters, i.e. $$\underset{(q+p)\times (q+p)}{\mathbf{A}}=\left(\begin{array}{c}\mathit{\gamma }\\ \mathbf{L}\end{array}\right)\text{.}$$ Then, the vector $\mathbf{b}$ is constructed by combining two basis vectors (Equation 31.1), to ensure it is equal to one in the position 1 and in the position p+1, i.e. $$\mathbf{b}=\left(\begin{array}{c}{\mathbf{e}}_p\\ {\mathbf{e}}_q\end{array}\right)\text{.}$$

Example: matrix form of an ARMA(2,3)

Example 20.2 For example let’s consider an ARMA(2,3): $$\mathit{\gamma }=\left(\begin{array}{ccccc}{\varphi }_1& {\varphi }_2& {\theta }_1& {\theta }_2& {\theta }_3\end{array}\right)\text{.}$$ The matrix for the AR part that is equal to $\mathbf{\Phi }$ combined with a matrix of zeros, i.e. $$\underset{1\times 2}{{\mathbf{L}}^{\text{AR}}}=\left(\begin{array}{cc}1& 0\end{array}\right)\text{,}$$ For the MA part we consider a similar matrix as in the single case but with first row equal to zero, i.e. $$\underset{2\times 3}{{\mathbf{L}}^{\text{MA}}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right)\text{,}$$ and we combine $$\underset{4\times 5}{\mathbf{L}}=\left(\begin{array}{cc}{\mathbf{L}}^{\text{AR}}& {\mathbf{0}}_{2\times 3}\\ {\mathbf{0}}_{1\times 2}& {\mathbf{0}}_{1\times 3}\\ {\mathbf{0}}_{2\times 2}& {\mathbf{L}}^{\text{MA}}\end{array}\right)=\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\end{array}\right)\text{.}$$ Finally we combine $$\underset{5\times 5}{\mathbf{A}}=\left(\begin{array}{c}\mathit{\gamma }\\ \mathbf{L}\end{array}\right)=\left(\begin{array}{ccccc}{\varphi }_1& {\varphi }_2& {\theta }_1& {\theta }_2& {\theta }_3\\ 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\end{array}\right)\text{.}$$ Then, the vector $\mathbf{b}$ is constructed by combining two vectors $$\mathbf{b}=\left(\begin{array}{c}{\mathbf{e}}_2\\ {\mathbf{e}}_3\end{array}\right)={\left(\begin{array}{ccccc}1& 0& 1& 0& 0\end{array}\right)}^{\mathrm{\top }}\text{.}$$ Let’s check it would lead to a classic ARMA(2,3): $$\begin{array}{rl}\left(\begin{array}{c}{Y}_t\\ Y_{t-1}\\ u_t\\ u_{t-1}\\ u_{t-2}\end{array}\right)& =\left(\begin{array}{c}{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+{\theta }_1{u}_{t-1}+{\theta }_2{u}_{t-1}+{\theta }_2{u}_{t-2}\\ Y_{t-1}\\ 0\\ u_{t-1}\\ u_{t-2}\end{array}\right)+\left(\begin{array}{c}{u}_t\\ 0\\ u_t\\ 0\\ 0\end{array}\right)\\ & =\left(\begin{array}{c}{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+{\theta }_1{u}_{t-1}+{\theta }_2{u}_{t-1}+{\theta }_3{u}_{t-2}+u_t\\ Y_{t-1}\\ u_t\\ u_{t-1}\\ u_{t-2}\end{array}\right)\end{array}$$

20.2 Moments

Proposition 20.1 Given the information at time $t$, the forecasted value at time $t+h$ for an AR(p) process reads: $$\begin{array}{}\text{(20.2)}& {\mathbf{X}}_{t+h}=\sum _{j=0}^{h-1}{\mathbf{A}}^j\mathbf{c}+{\mathbf{A}}^h{\mathbf{X}}_t+\sum _{j=0}^{h-1}{\mathbf{A}}^j\mathbf{b}{u}_{t+h-j}\text{.}\end{array}$$

Proof of Proposition 20.1

Proof. Let’s start developing the Equation 20.2 with $h=1$, i.e.  $${\mathbf{X}}_{t+1}=\mathbf{c}+\mathbf{A}{\mathbf{X}}_t+\mathbf{b}{u}_{t+1}\text{,}$$ with $h=2$: $$\begin{array}{rl}{\mathbf{X}}_{t+2}& =\mathbf{c}+\mathbf{A}{\mathbf{X}}_{t+1}+\mathbf{b}{u}_{t+2}=\\ & =\mathbf{A}(\mathbf{c}+\mathbf{A}{\mathbf{X}}_t+\mathbf{b}{u}_{t+1})+\mathbf{b}\cdot u_{t+2}=\\ & =\mathbf{c}+\mathbf{A}\mathbf{c}+{\mathbf{A}}^2{\mathbf{X}}_t+\mathbf{b}{u}_{t+2}+\mathbf{A}\mathbf{b}{u}_{t+1}\end{array}$$ with $h=3$: $$\begin{array}{rl}{\mathbf{X}}_{t+3}& =\mathbf{c}+\mathbf{A}{\mathbf{X}}_{t+2}+\mathbf{b}{u}_{t+3}=\\ & =\mathbf{c}+\mathbf{A}(\mathbf{c}+\mathbf{A}\mathbf{c}+{\mathbf{A}}^2{\mathbf{X}}_t+\mathbf{b}{u}_{t+2}+\mathbf{A}\mathbf{b}{u}_{t+1})+\mathbf{b}{u}_{t+3}=\\ & =\mathbf{c}+\mathbf{A}\mathbf{c}+{\mathbf{A}}^2\mathbf{c}+{\mathbf{A}}^3{\mathbf{X}}_t+\mathbf{b}{u}_{t+3}+\mathbf{A}\mathbf{b}{u}_{t+2}+{\mathbf{A}}^2\mathbf{b}{u}_{t+1}\end{array}$$ and so on. Hence the impact of the shocks decrease exponentially over time.

Proposition 20.2 ($\mathbf{\text{Expectation of an ARMA}}$)
The conditional expectation at time $t+h$ of ${\mathbf{X}}_{t+h}$ given the information up to time $t$ can be easily computed from Equation 20.2, i.e.  $$\mathbb{E}\{{\mathbf{X}}_{t+h}\mid {\mathcal{F}}_t\}=\mathbf{c}({\mathbf{I}}_p-{\mathbf{A}}^h)({\mathbf{I}}_p-\mathbf{A}{)}^{-1}+{\mathbf{A}}^h{\mathbf{X}}_t\text{.}$$

Proof of Proposition 20.2

Proof. The expectation of Equation 20.2 reads: $$\mathbb{E}\{{\mathbf{X}}_{t+h}\mid {\mathcal{F}}_t\}=\sum _{j=0}^{h-1}{\mathbf{A}}^j\mathbf{c}+{\mathbf{A}}^h{\mathbf{X}}_t+\sum _{j=0}^{h-1}{\mathbf{A}}^j\mathbf{b}\cdot \mathbb{E}\{u_{t+h-j}\mid {\mathcal{F}}_t\}\text{.}$$ Under the assumption that $u_{t+h-j}\sim \text{MDS}(0,1)$, we have that for all $t$ $$\mathbb{E}\{u_{t+1}\mid {\mathcal{F}}_t\}=0\text{.}$$ Therefore $$\mathbb{E}\{{\mathbf{X}}_{t+h}\mid {\mathcal{F}}_t\}=\sum _{j=0}^{h-1}{\mathbf{A}}^j\mathbf{c}+{\mathbf{A}}^h{\mathbf{X}}_t\text{,}$$ where for constant $\mathbf{c}$, the series further simplifies in $$\sum _{j=0}^{h-1}{\mathbf{A}}^j\mathbf{c}=({\mathbf{I}}_p-{\mathbf{A}}^h)({\mathbf{I}}_p-\mathbf{A}{)}^{-1}\mathbf{c}\text{,}$$ with ${\mathbf{I}}_p$ the identity matrix (Equation 31.3).

Proposition 20.3 ($\mathbf{\text{Variance of an ARMA}}$)
Taking the variance of Equation 20.2 on both sides gives $$\mathbb{V}\{{\mathbf{X}}_{t+h}\mid {\mathcal{F}}_t\}={\sigma }_{\text{u}}^2\sum _{j=0}^{h-1}{\mathbf{A}}^j\mathbf{b}{\mathbf{b}}^{\mathrm{\top }}({\mathbf{A}}^j{)}^{\mathrm{\top }}\text{.}$$ Moreover, Assuming $u_t\sim \text{WN}(0,{\sigma }_{\text{u}}^2)$, and independence across time, the conditional covariance is: $$\mathbb{C}v\{{\mathbf{X}}_{t+h},{\mathbf{X}}_{t+k}\mid {\mathcal{F}}_t\}={\sigma }_{\text{u}}^2\sum _{j=0}^{min(h,k)-1}{\mathbf{A}}^{h-1-j}\mathbf{b}{\mathbf{b}}^{\mathrm{\top }}({\mathbf{A}}^{k-1-j}{)}^{\mathrm{\top }}.$$

Proof of Proposition 20.3

Proof. Taking the conditional variance gives $$\mathbb{V}\{{\mathbf{X}}_{t+h}\mid {\mathcal{F}}_t\}=\sum _{j=0}^{h-1}\mathbb{V}\{{\mathbf{A}}^j\mathbf{b}{u}_{t+h-j}\mid {\mathcal{F}}_t\}=\sum _{j=0}^{h-1}{\mathbf{A}}^j\mathbf{b}\mathbb{V}\{u_{t+h-j}\mid {\mathcal{F}}_t\}{\mathbf{b}}^{\mathrm{\top }}({\mathbf{A}}^j{)}^{\mathrm{\top }}\text{.}$$ Then, if $u_t$ is White Noise process, then its variance is constant, hence independent from $t$, and equal to ${\sigma }_{\text{u}}^2$. For the covariance formula, each ${\mathbf{X}}_{t+h}$ is influenced by past shocks $u_{t+h-j}$. Since the shocks are uncorrelated across time, only shared shocks affect both ${\mathbf{X}}_{t+h}$ and ${\mathbf{X}}_{t+k}$. The number of common shocks is $min(h,k)$, and each contributes: $$\mathbb{C}v\{{\mathbf{A}}^{h-1-j}\mathbf{b}{u}_{t+1+j},{\mathbf{A}}^{k-1-j}\mathbf{b}{u}_{t+1+j}\}={\sigma }_{\text{u}}^2{\mathbf{A}}^{h-1-j}\mathbf{b}{\mathbf{b}}^{\mathrm{\top }}({\mathbf{A}}^{k-1-j}{)}^{\mathrm{\top }}\text{.}$$

Example: ARMA(2,3) iterative forecast

Example 20.3 Let’s use Monte Carlo simulations to establish if the results are accurate.

# Simulate an ARMA(p, q) process
simulate_ARMA <- function(t_bar, nsim = 1, mu = 0, phi = c(0.5, 0.2), theta = c(0.4, 0.1), sigma2 = 1, y0, eps0){
  scenarios <- list()
  # MA order
  q <- length(theta)
  # AR order 
  p <- length(phi)
  # Initial lag to start 
  i_start <- max(c(p, q)) + 1
  for(sim in 1:nsim){
    # Simulated residuals 
    if (!missing(eps0)){
      eps <- c(eps0, rnorm(t_bar - q, sd = sqrt(sigma2)))
    } else {
      eps <- rnorm(t_bar, sd = sqrt(sigma2))
    }
    # Initialize the simulation 
    if (!missing(y0)){
      x <- c(y0)
    } else {
      x <- eps
    }
    for(t in i_start:t_bar){
      # AR part 
      ar_part <- 0
      for(i in 1:length(phi)){
        ar_part <- ar_part + phi[i] * x[t-i]
      }
      # MA part 
      ma_part <- 0
      for(i in 1:length(theta)){
        ma_part <- ma_part + theta[i] * eps[t-i]
      }
      # Next step simulation 
      x[t] <- mu + ar_part + ma_part + eps[t]
    }
    scenarios[[sim]] <- x
  }
  scenarios
}

# Compute the power of a matrix 
pow_matrix <- function(x, n = 0){
  M <- x
  # Case power == 0
  if (n == 0) {
    M <- diag(1, nrow(M), ncol(M))
  } else if(n > 1) {
    M_pow <- M
    for(i in 2:n){
      M_pow <- M_pow %*% M
    }
    M <- M_pow
  }
  return(M)
}

# Compute the expected value with the iterative formula 
expectation_ARMA <- function(X0, h = 10, A, b, c){
  A_h <- pow_matrix(A, h)
  I_p <- diag(1, nrow = nrow(A_h))
  (I_p - A_h) %*% solve(I_p - A) %*% c  + A_h %*% X0
}

# Compute the covariance with the iterative formula 
covariance_ARMA <- function(h, k, A, b, sigma2 = 1){
  cv_x <- 0
  for(j in 0:(min(h, k)-1)){
    A_hj <- pow_matrix(A, h-1-j)
    A_kj <- pow_matrix(A, k-1-j)
    cv_x <- cv_x + A_hj %*% b %*% t(b) %*% t(A_kj)
  }
  sigma2 * cv_x
}

# *************************************************
#                      Inputs 
# *************************************************
set.seed(6) # random seed
nsim <- 200 # Number of simulations
t_bar <- 100000 # Steps for each simulation 
phi <- c(0.3, 0.15) # AR parameters
theta <- c(0.3, 0.15, 0.1) # MA parameters
mu <- 0.3 # Intercept 
sigma2 <- 1 # Variance of the residuals 
y0 <- c(0.9, 0.2, 0) # Initial values for y0
eps0 <- c(0.3, -0.2, 0.1) # Initial values for eps0
h <- 200 # Number of steps ahead
# *************************************************
#                   Simulations 
# *************************************************
scenarios <- simulate_ARMA(t_bar, nsim, mu, phi, theta, 
                           sigma2, y0, eps0)
# Build companion matrix 
L_AR <- matrix(c(1, 0, 0, 0, 0), nrow = 1)
L_MA <- cbind(c(0, 0), c(0, 0), diag(2), c(0, 0))
A <- rbind(c(phi, theta), L_AR, rep(0, 5), L_MA)
# Build vector b
b <- matrix(c(1, 0, 1, 0, 0), ncol = 1)
# Build vector c
c <- matrix(c(mu, 0, 0, 0, 0), ncol = 1)
# *************************************************
#                      Moments  
# *************************************************
# Compute expectation and variances with Monte Carlo
e_mc <- mean(purrr::map_dbl(scenarios, ~mean(.x)))
v_mc <- mean(purrr::map_dbl(scenarios, ~var(.x)))
cv_mc_L1 <- mean(purrr::map_dbl(scenarios, ~cov(.x[-1], lag(.x, 1)[-1])))
cv_mc_L2 <- mean(purrr::map_dbl(scenarios, ~cov(.x[-c(1,2)], lag(.x, 2)[-c(1,2)])))
cv_mc_L3 <- mean(purrr::map_dbl(scenarios, ~cov(.x[-c(1:3)], lag(.x, 3)[-c(1:3)])))
cv_mc_L4 <- mean(purrr::map_dbl(scenarios, ~cov(.x[-c(1:4)], lag(.x, 4)[-c(1:4)])))
cv_mc_L5 <- mean(purrr::map_dbl(scenarios, ~cov(.x[-c(1:5)], lag(.x, 5)[-c(1:5)])))
# Compute expectation and variances with Formula
# State vector for forecasting 
X0 = c(y0[-3][c(2,1)], eps0[c(3:1)])
e_mod <- expectation_ARMA(X0, h = h, A, b, c)[1,1]
v_mod <- covariance_ARMA(h, h, A, b, sigma2)[1,1]
cv_mod_L1 <- covariance_ARMA(h, h-1, A, b, sigma2 = 1)[1,1]
cv_mod_L2 <- covariance_ARMA(h, h-2, A, b, sigma2 = 1)[1,1]
cv_mod_L3 <- covariance_ARMA(h, h-3, A, b, sigma2 = 1)[1,1]
cv_mod_L4 <- covariance_ARMA(h, h-4, A, b, sigma2 = 1)[1,1]
cv_mod_L5 <- covariance_ARMA(h, h-5, A, b, sigma2 = 1)[1,1]

Covariance	Formula	MonteCarlo
$\mathbb{E}\{Y_t\}$	0.5454545	0.5451268
$\mathbb{V}\{Y_t\}$	1.7466575	1.7472189
$\mathbb{C}v\{Y_t,Y_{t-1}\}$	1.1317615	1.1322251
$\mathbb{C}v\{Y_t,Y_{t-2}\}$	0.8115271	0.8118819
$\mathbb{C}v\{Y_t,Y_{t-3}\}$	0.5132223	0.5133504
$\mathbb{C}v\{Y_t,Y_{t-4}\}$	0.2756958	0.2758340
$\mathbb{C}v\{Y_t,Y_{t-5}\}$	0.1596921	0.1596755

Table 20.1: Theoric long term moments and moments computed on 200 Monte Carlo simulations (t = 100000).

Figure 20.2: ARMA(2,3) simulations with expected value (red) and confidence intervals with $\alpha =0.1$ (green), $\alpha =0.05$ (magenta) and $\alpha =0.01$ (purple).