We have a strong definition of our ARMA process now, but we’ll also want to quantify its statistics and parameters given data. First, let’s talk about how to estimate our two main statistics of interest: the mean and the covariance function. Throughout we will assume the process is stationary. The sample mean is the easiest, it can be estimated by:

where T is the length of our sample. The sample covariance is a bit more involved:

Note that the effective size of our sample is smaller when we calculate the sample covariance because we require an offset of h. Regardless we still divide by T because this estimator has a few nicer properties than the estimator we get if we divide by T-h (see Shumway and Stoffer). We can also calculte the sample correlation from the sample covariance:

\begin{equation*} \hat{\rho}(h)=\frac{\hat{\gamma}(h)}{\hat{\gamma}(0)} \tag{3.3} \end{equation*}

These estimators are great, but unless T \rightarrow \infty they are not going to give us the true values. It’s useful to understand the expected error of these estimators to have a handle on the significance of the measurement. For our correlation function, the best metric is a comparison to white noise (which we know has a correlation coefficient of 0 for all h>0 ). For a large sample and finite h, the sample ACF of white noise is normally distributed with mean 0 and variance:

\begin{equation*} \sigma_{\hat{\rho}(h)}^{2}=\frac{1}{T} \tag{3.4} \end{equation*}

It is common to declare a measurement of \hat{\rho}(h) significant if:

\begin{equation*} |\hat{\rho}(h)|>\frac{2}{\sqrt{T}} \tag{3.5} \end{equation*}

This is equivalent to the 95^{\text {th }} percentile.

Let’s combine two concepts and see how the backshift operator allows us calculate the correlation coefficient for an AR(2) process. Let’s start with our equation:

\begin{equation*} X_{t+h}=\phi_{1} X_{t+h-1}+\phi_{2} X_{t+h-2}+W_{t+h} \tag{3.6} \end{equation*}

Multiply through by X_{t} and taking an expectation:

For a stationary process, we can substitute the covariance terms, recalling that \mathbb{E}\left[X_{t+h} X_t\right]=\gamma(h). The final term, \mathbb{E}\left[W_{t+h} X_t\right], is zero for h>0 because the white noise term W_{t+h} is uncorrelated with all past values of X_t. This gives us the key result:

This is known as a difference equation. We can rewrite this in terms of our backshift operator:

\begin{align*} P(B) \gamma(h) & =0 \tag{3.9} \end{align*}

The roots of our backshift operator are the same roots of our difference equation. The solutions of a difference equation of order p:

\begin{equation*} u_{n}-\alpha_{1} u_{n-1}-\ldots-\alpha_{p} u_{n-p}=0 \tag{3.10} \end{equation*}

are given by:

\begin{equation*} u_{n}=z_{1}^{-n} P_{1}(n)+z_{2}^{-n} P_{2}(n)+\ldots+z_{r}^{-n} P_{r}(n), \tag{3.11} \end{equation*}

where r is the number of distinct roots and P_{j}(n) is a polynomial in n of degree m_{j}-1 where m_{j} is the multiplicity of the root z_{j}. For \operatorname{AR}(2) we have a difference equation of order 2. If there are two distinct roots, z_{1} \neq z_{2}, then the solution has the form:

\begin{equation*} \rho(h)=c_{1} z_{1}^{-|h|}+c_{2} z_{2}^{-|h|} \tag{3.12} \end{equation*}

If instead the two roots are equal then the solution has the form:

\begin{equation*} z_{1}^{-|h|}\left(c_{1}+c_{2}|h|\right) . \tag{3.13} \end{equation*}

For higher orders there is a wider variety of solutions.

Let’s imagine we have a process we suspect is a MA process. If we have a reliable estimate of the ACF we can immediately tell its order by the largest non-zero value. For an AR process that’s a lot less clear. We have already shown that for an AR(2) process the ACF is non-zero for all separations. More generically, even if X_{t+h} does not directly depend on X_{t} in the equation that defines the process, it can still be correlated through some intermediary X_{t+k} with k<h. What we would like to know is the correlation between X_{t+h} and X_{t} conditioned on X_{t+1: t+h-1} :

\begin{equation*} \gamma_{X_{t+h}, X_{t} \mid X_{t+1: t+h-1}} \tag{3.14} \end{equation*}

If we assume that our process is Gaussian, between X_{t+h} and X_{t} after having subtracted the optimal linear predictor using X_{t+1: t+h-1}. We will call these linear predictor \hat{X}_{t+h} and \hat{X}_{t}. This is known as the partial auto-correlation function:

For now we will assert that the PACF of an AR(p) process is 0 for h>p.

Once we’ve selected what model to use on our data, we will want to estimate its parameters. For now, let’s introduce the simplest choice we can make in Bayesian inference: the maximum a posteriori (MAP) estimate. Assume we have some data X and some model parameters \theta. The MAP estimate is given by:

We already have all the tools we need to evaluate this on our ARMA processes. We will see an example of prior choices in lab, and we can write the likelihood from the equations that define the process. For an AR(1) process we have two parameters, \phi_{1} and \sigma_{w}. If we set the boundary condition X_{1}=W_{1}, we can write:

\begin{align*} & p\left(X_{1: t} \mid \phi_{1}, \sigma_{w}\right)=p\left(X_{1} \mid \phi_{1}, \sigma_{w}\right) p\left(X_{2} \mid X_{1}, \phi_{1}, \sigma_{w}\right) \times \\ & p\left(X_{3} \mid X_{1: 2}, \phi_{1}, \sigma_{w}\right) \ldots p\left(X_{t} \mid X_{1: t-1}, \phi_{1}, \sigma_{w}\right) . \tag{3.28} \end{align*}

But we know that for an \mathrm{AR}(1) process X_{t} is conditionally independent of X_{1: t-2} given X_{t-1} so we can simplify this equation to:

The likelihood of X_{t} given X_{t-1} and the parameters is just a Gaussian distribution with mean \phi_{1} X_{t-1} and variance \sigma_{w}^{2} :

\begin{align*} p\left(X_{t} \mid X_{t-1}, \phi_{1}, \sigma_{w}\right) & =\mathcal{N}\left(\phi_{1} X_{t-1}, \sigma_{w}\right) \tag{3.30}\\ & =\frac{1}{\sqrt{2 \pi} \sigma_{w}} \exp \left[\frac{-\left(X_{t}-\phi_{1} X_{t-1}\right)^{2}}{2 \sigma_{w}^{2}}\right] \tag{3.31} \end{align*}

You will need to extend this to \operatorname{AR}(2) for lab this week, but I will spoil the answer and tell you that it’s just a matter of adding \phi_{2} X_{t-2} to the mean.

One final note - when we take the product of probability distributions we quickly get vanishingly small numbers. From an implementation standpoint, this will run into the numerical precision of our computer. To avoid this, we often operate with log pdfs instead of the pdf itself and take advantage of the property that:

\begin{equation*} \log [p(x) p(y)]=\log [p(x)]+\log [p(y)] \tag{3.32} \end{equation*}