3 Linear models for stationary time series.

3.1.1 Example: The MA($1$) model. Let

$$x_t=\mu +e_t-\theta e_{t-1}$$

where $\mu ={\mu }_x$ and $e_t$ is white noise. The reason for the $-$ sign will beome evident.

The properties of this model are

$$\text{E}\left(x_t\right)=\mu \text{ and Var}\left(x_t\right)=\left(1+{\theta }^2\right){\sigma }_e^2$$

$${\rho }_{x,1}=-\frac{\theta }{1+{\theta }^2}\text{ and }{\rho }_{x,k}=0\text{ for }k>1.$$

Proof as an exercise.

Remarks. $\left|{\rho }_{x,1}\right|\le \frac{1}{2}$, follows from ${\left(1\pm \theta \right)}^2\ge 0\Rightarrow \left(1+{\theta }^2\right)\ge \mp 2\theta$. Also

$$\frac{\theta }{1+{\theta }^2}=\frac{{\theta }^{-1}}{1+{\left({\theta }^{-1}\right)}^2}$$

so that replacing $\theta$ by ${\theta }^{-1}$ would result in the same acf. Provided $\left|\theta \right|\ne 1$ we therefore have a choice of values for $\theta$ when modelling an acf of this form. The convention is to choose the value that satisfies .

3.1.2 Definition. A time series $x_t$ follows a moving average model of order $q$ - MA($q$) if it may be represented as

$$x_t=\mu +e_t-{\theta }_1{e}_{t-1}-\mathrm{\cdots }-{\theta }_q{e}_{t-q}$$

where $e_t$ is a white noise series.

Remarks. The model has $q$ moving average parameters ${\theta }_j$ and one variance parameter ${\sigma }_e^2$. For some purposes it is useful to have a slightly different notation and to write
3.1.3

$$x_t=\mu +{\psi }_0{e}_t+{\psi }_1{e}_{t-1}+\mathrm{\cdots }+{\psi }_q{e}_{t-q}$$

where ${\psi }_0=1$ and ${\psi }_j=-{\theta }_j$ for $j=1\mathrm{\dots }q$. The name moving average is more usually applied when one observed time series takes the place of $e_t$ in this equation, and a new series, taking the place of $x_t$, is constructed according to this equation. Our model contrasts with this in that $x_t$ is observed, and we are proposing a model which gives rise to its autocorrelation pattern.

This is a convenient point to introduce an operator notation for discrete time series.

3.1.4 Definition. The backward shift operator $B$ applied to a series $e_t$ replaces it by the shifted series with values $f_t=e_{t-1}$. Powers of $B$ are similarly defined, i.e.

$$B{e}_t=e_{t-1}\text{ and }{B}^m{e}_t=e_{t-m}.$$

Using this notation the MA($q$) model may be expressed from Definition 3.1.2 as
3.1.5

	$x_t$	$=$	$\mu +e_t-{\theta }_{1}B{e}_t-{\theta }_2{B}^2{e}_t-\mathrm{\cdots }-{\theta }_q{B}^q{e}_t$
		$=$	$\mu +\left(1-{\theta }_{1}B-{\theta }_{2}B-\mathrm{\cdots }-{\theta }_q{B}^q\right)e_t$
		$=$	$\mu +\theta \left(B\right)e_t$

For example the MA($1$) model is

$$x_t=\mu +\left(1-\theta B\right)e_t.$$

Now we can formally re-arrange this, and provided expand:

$$e_t={\left(1-\theta B\right)}^{-1}\left(x_t-\mu \right)=\left(1+\theta B+{\theta }^2{B}^2+{\theta }^3{B}^3+\mathrm{\cdots }\right)\left(x_t-\mu \right)=$$

$$\left(x_t-\mu \right)+\theta \left(x_{t-1}-\mu \right)+{\theta }^2\left(x_{t-2}-\mu \right)+{\theta }^3\left(x_{t-3}-\mu \right)+\mathrm{\cdots }.$$

We can directly verify from this expression that $e_t-\theta e_{t-1}=x_t-\mu$.

The fact that the unobserved series $e_t$ can be ‘recovered’ from the observed series $x_t$ is the reason for choosing . The general result is:

3.1.6 Definition. The MA($q$) model of 3.1.5 is said to be invertible if it is possible to express $e_t$ in terms of the present and past values of $x_t$ as

	$e_t$	$=$	$\left(x_t-\mu \right)-{\pi }_1\left(x_{t-1}-\mu \right)-{\pi }_2\left(x_{t-2}-\mu \right)-\mathrm{\cdots }$
		$=$	$\left(1-{\pi }_{1}B-{\pi }_2{B}^2-\mathrm{\cdots }\right)\left(x_t-\mu \right)$
		$=$	$\pi \left(B\right)\left(x_t-\mu \right).$

Remark. The signs in this last equation are chosen for convenience on rewriting the expansion as
3.1.7

$$x_t=\mu +{\pi }_1\left(x_{t-1}-\mu \right)+{\pi }_2\left(x_{t-2}-\mu \right)+\mathrm{\cdots }+e_t$$

- a form of the model which will be useful for prediction.

3.1.8 Definition. The invertibility condition for model 3.1.5 is that the operator $\theta \left(B\right)$ satisfies

$$\theta \left(B\right)\ne 0\text{ for }\left|B\right|\le 1.$$

Here, we are treating $B$ as a real or complex number rather than an operator.

Remark. Consider a MA(2) model where we observe $\left\{x_t\right\}$ satisfying $x_t=e_t-0.7e_{t-1}-0.6e_{t-2}$. The roots of $\theta \left(B\right)=1-0.7B-0.6B^2=0$ are given by $B=5/6$ and $B=-2$ so that $\theta \left(B\right)$ is not invertible.

The invertibility condition means that it is possible to expand

$$\theta {\left(B\right)}^{-1}=\pi \left(B\right)=1-{\pi }_{1}B-{\pi }_2{B}^2-\mathrm{\cdots }$$

with ${\pi }_j\to 0$ in order to derive 3.1.6 from 3.1.5.

One way to derive the coefficient $\left\{{\pi }_j\right\}$ of $1/\theta \left(B\right)$ is to use partial fractions which is standard in power-series expansion.

We now consider when it is appropriate to use a moving average model to represent a given series. This is based on the autocorrelation properties of the model:

3.2.1 Recall that ${\psi }_j=-{\theta }_j$ so that ${\psi }_0=1$, ${\psi }_1=-{\theta }_1$, ${\psi }_q=-{\theta }_q$. The autocovariances of a moving average model are given in terms of the coefficients ${\psi }_j$ as

$${\gamma }_{x,k}={\sigma }_e^2\left({\psi }_0{\psi }_k+\mathrm{\cdots }+{\psi }_{q-k}{\psi }_q\right)={\sigma }_e^2\left(\sum _{i=0}^{q-k}{\psi }_i{\psi }_{k+i}\right)\text{ for }k=0\mathrm{\dots }q.$$

The Characteristic property of the model is that
3.2.2

$${\rho }_k=0\text{ for }k>q.$$

Continuing the previous example, taking ${\sigma }_e^2=10$, we get

so that ${\rho }_1=-2.8/18.5=-0.15$ and ${\rho }_2=-6.0/18.5=-0.32$.

Notice that the above formula holds even if $q=\mathrm{\infty }$. In other words, if

$$x_t=\sum _{j=0}^{\mathrm{\infty }}{\psi }_j{e}_{t-j},$$

where $\left\{e_t\right\}$ is white noise, then

$${\gamma }_{x,k}={\sigma }_e^2\left(\sum _{j=0}^{\mathrm{\infty }}{\psi }_j{\psi }_{j+k}\right),{\rho }_{x,k}=\frac{{\sum }_{j=0}^{\mathrm{\infty }}{\psi }_j{\psi }_{j+k}}{{\sum }_{j=0}^{\mathrm{\infty }}{\psi }_j^2}.$$

There is a converse result, that any stationary time series $x_t$ having the property 3.2.2 may be represented by a moving average model 3.1.2. satisfying the invertibility condition.

The characteristic property is recognised by inspection of the sample acf of the data series. The sampling properties of this acf are therefore required, and these will be given shortly. But first have a look at two time series for which Moving Average models appear to be appropriate. In fact it is the successive changes, or first differences, of these series which appear to have the characteristics of Moving Average models.

Consider the ‘CABLE’ data, available on Moodle, which is clearly non-stationary. A good model for this data might be of the form $y_t=x_t{e}^{\left(a+bt\right)}$ where $\left\{x_t\right\}$ is stationary. We therefore take the logarithm and then the first differences of the series of weekly business transactions as shown in Figure 6 to make it stationary.

Figure 7 shows the sample time series properties of the differenced series, which are typical of an MA(1) model:

The model $x_t=\mu +e_t-0.713e_{t-1}$ is estimated for this series, and it has the model properties shown in the following figures.

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

The second plot shows the implied ${\psi }_k$ weights of the MA representation of 3.1.3, and the third plot shows the implied ${\pi }_k$ weights of the MA representation given in 3.1.7. The fourth plot shows the implied model spectrum which we shall define in Section 3.12.

The second example is of the annual changes in daylength, measured in milliseconds, from the dataset ‘DAYLENGT’ available in Moodle. Part of what we see is almost certainly due to the fact that the original data has been smoothed to give the series shown in the left panel of Figure 9. Each point has been replaced by its average with the previous and next point. And this procedure has then been repeated. This procedure has the property of completely annihilating any cycle of period 3.

Figure 10 shows the sample time series properties of the differenced series. Among these, the sample spectrum is very low at frequency 1/3 as a consequence of the smoothing. The sample autocorrelations are characteristic of an MA(3) model:

The MA(3) model $x_t=\mu +e_t+0.89e_{t-1}+0.87e_{t-2}+0.23e_{t-3}$ is estimated for this series, and it has the model properties shown in the following figures, which match the sample properties very well:

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

In general the sample value $r_{x,k}$ is an approximately unbiased and consistent estimate of ${\rho }_{x,k}$ for any fixed lag $k$ as the sample size $n$ tends to infinity. In the particular case of the MA($q$) model there is a useful result for the standard error of the sample acf for lags greater than the order. This helps us to weigh the statistical evidence for a hypothesised value of q.

3.2.3 For a Gaussian time series of length $n$ with acf satisfying ${\rho }_k=0$ for $k>q$

$$SE\left(r_k\right)\approx {\left\{\frac{1}{n}\left[1+2\sum _{j=1}^q{\rho }_j^2\right]\right\}}^{\frac{1}{2}}\text{ for }k>q.$$

The simplest case, $q=0$, gives the property under the hypothesis of white noise, that

$$SE\left(r_k\right)\approx \frac{1}{\sqrt{n}}$$

so it is usual to plot the sample correlations on a graph with approximate $95\%$ error limits drawn at $\pm 2/\sqrt{n}$ about the horizontal axis. If the sample acf values generally lie within these limits with the only values outside them being at unremarkable lags then white noise is a sensible conclusion. But these only apply for assessing the hypothesis of white noise. Under any other hypothesis, the limits are wider.

The sample acf may in fact lie outside these limits for many lags. We are here interested in the case when a small number of values lie clearly outside the limits at low lags. We then have good reason to suppose that a MA($q$) model is appropriate with the value of $q$ being tentatively chosen as the ‘cut-off’ lag - the highest lag with a ‘significant’ acf value.

We need then to be aware of the fact that the sampling variability of the remaining acf values, at lags greater than the assumed MA order $q$ will be increased according to 3.2.3, so that rather more values will lie outside the white noise limits. The increased $SE$ may be estimated from 3.2.3 using the sample acf values $r_j$ up to lag $q$ in place of the unknown ${\rho }_j$.

Another visual characteristic of the sample acf at these higher lags is that their statistical variablity may follow a smoother pattern than for white noise. Even when the true model is MA($q$) these patterns may be deceptively suggestive of structure in the acf, i.e. non zero values of ${\rho }_k$, which is not actually present at these higher lags.

Estimation of the parameters in a MA($q$) model may sometimes be achieved by solving the equations 3.2.1 for the unknown parameters ${\sigma }_e^2$ and ${\psi }_j=-{\theta }_j,j=1\mathrm{\dots }q$ using sample autocovariance values $C_{x,k}$ in place of ${\gamma }_k$. Unfortunately these equations may not have a real solution, though when they do it is easily found. For example in the MA($1$) case the solution for $\theta$ is found by solving

$$r_{x,1}=-\frac{\theta }{1+{\theta }^2}$$

which has no real solution if $r_{x,1}$ is outside the range $[-\frac{1}{2},\frac{1}{2}]$, and this may well happen through statistical variability in $r_{x,1}$.

Consider the model
3.3.1

$$x_t={\mu }_x+\varphi \left(x_{t-1}-{\mu }_x\right)+e_t$$

which defines successive values of the series $x_t$ in terms of previous values and a white noise series $e_t$. This is the first order autoregessive model, commonly referred to as the AR(1) model. The essential distinction with prediction models is that for the prediction equation $e_t$ is only required to be uncorrelated with $x_{t-1}$ - in general it will not be white noise. For the AR(1) model however, $e_t$ is assumed to be white noise and will, as a result of that, be independent of $x_{t-1}$.

The properties of this model may be found by successively substituting

$$x_t-{\mu }_x=e_t+\varphi \left(x_{t-1}-{\mu }_x\right)=e_t+\varphi \left(e_{t-1}+\varphi \left(x_{t-2}-{\mu }_x\right)\right)=\mathrm{\dots }$$

$$=e_t+\varphi e_{t-1}+{\varphi }^2{e}_{t-2}+\mathrm{\cdots }+{\varphi }^L{e}_{t-L}+{\varphi }^{L+1}\left(x_{t-L-1}-{\mu }_x\right)$$

Now assuming that both and $x_t$ is bounded ‘in the infinite past’ we can let $L\to \mathrm{\infty }$ and obtain the convergent sum:
3.3.2

$$x_t={\mu }_x+e_t+\varphi e_{t-1}+{\varphi }^2{e}_{t-2}+\mathrm{\cdots }={\mu }_x+\sum _{i=0}^{\mathrm{\infty }}{\varphi }^i{e}_{t-i}.$$

We can derive from this, for example, E($x_t$)=${\mu }_x$ and

$${\sigma }_x^2=\text{Var}\left(x_t\right)={\sigma }_e^2+{\varphi }^2{\sigma }_e^2+{\varphi }^4{\sigma }_e^2+\mathrm{\cdots }=\frac{{\sigma }_e^2}{1-{\varphi }^2}.$$

The representation 3.3.2 could in fact be obtained using the backward shift operator notation by writing 3.3.1 as $\left(1-\varphi B\right)\left(x_t-{\mu }_x\right)=e_t$ so that

$$x_t-{\mu }_x={\left(1-\varphi B\right)}^{-1}{e}_t=\left(1+\varphi B+{\varphi }^2{B}^2+\mathrm{\cdots }\right)e_t.$$

The variance and correlation properties can be derived directly from the model 3.3.1 without using this expansion. Because $x_{t-1}$ and $e_t$ are uncorrelated we get

$${\sigma }_x^2=\text{Var}\left(x_t\right)=\text{Var}\left(\varphi x_{t-1}+e_t\right)={\varphi }^2{\sigma }_x^2+{\sigma }_e^2\Rightarrow \left(1-{\varphi }^2\right){\sigma }_x^2={\sigma }_e^2.$$

Similarly $e_t$ is uncorrelated with $x_{t-k}$ for any $k>0$ so

	${\gamma }_{x,k}$	$=$	$\text{Cov}(x_t,x_{t-k})$
		$=$	$\text{Cov}(\varphi x_{t-1}+e_t,x_{t-k})$
		$=$	$\text{Cov}(\varphi x_{t-1},x_{t-k})$
		$=$	$\varphi \text{Cov}(x_{t-1},x_{t-k})$
		$=$	$\varphi {\gamma }_{x,k-1}.$

It follows that

$${\gamma }_{x,1}=\varphi {\sigma }_x^2,{\gamma }_{x,2}=\varphi {\gamma }_{x,1}={\varphi }^2{\sigma }_x^2,\mathrm{\dots },{\gamma }_{x,k}={\varphi }^k{\sigma }_x^2.$$

Note therefore that ${\rho }_{x,1}=\varphi$ and ${\rho }_{x,k}={\varphi }^k$ for $k>0$.

3.3.3 Definition. A time series follows an autoregressive model of order $p$ - AR(p), if it is generated by the process:

$$x_t={\mu }_x+{\varphi }_1\left(x_{t-1}-{\mu }_x\right)+{\varphi }_2\left(x_{t-2}-{\mu }_x\right)+\mathrm{\cdots }+{\varphi }_p\left(x_{t-p}-{\mu }_x\right)+e_t$$

where $e_t$ is a white noise series. Using operator notation the model may be written:
3.3.4

$$\left(1-{\varphi }_{1}B-{\varphi }_2{B}^2-\mathrm{\cdots }-{\varphi }_p{B}^p\right)\left(x_t-{\mu }_x\right)=\varphi \left(B\right)\left(x_t-{\mu }_x\right)=e_t.$$

We also require:
3.3.5 Definition. The stationarity condition for model 3.3.3 is that the operator $\varphi \left(B\right)$ satisfies

$$\varphi \left(B\right)\ne 0\text{ for }\left|B\right|\le 1.$$

Remark. The condition is the same as that required of $\theta \left(B\right)$ for a moving average model to be invertible (see 3.1.8). In this case it means that it is possible to carry out the convergent expansion following from 3.3.4:

	$\varphi \left(B\right)\left(x_t-{\mu }_x\right)$	$=$	$e_t$
	$\Rightarrow x_t-{\mu }_x$	$=$	$\varphi {\left(B\right)}^{-1}{e}_t$
		$=$	$\left(1+{\psi }_{1}B+{\psi }_2{B}^2+\mathrm{\cdots }\right)e_t$
	$\Rightarrow x_t$	$=$	${\mu }_x+{\displaystyle \sum _{i=0}^{\mathrm{\infty }}}{\psi }_i{e}_{t-i}$

where ${\psi }_j\to 0$. This is an infinite moving average representation of the model (compare with the finite represenation of 3.1.3 for MA processes). We note one consequence which may be obvious from the way the series is generated, but which we can formally verify with this expansion; the fact that $e_t$ is uncorrelated with all past values of the series:
3.3.6

$$\text{Cov}(e_t,x_{t-i})=0\text{ for }i\ge 1.$$

This is because $x_{t-i}={\mu }_x+e_{t-i}+{\psi }_1{e}_{t-i-1}+\mathrm{\cdots }$ in which all terms are uncorrelated with $e_t$.

The main contrast with moving average models is that the acf does not have a cut off, but generally decays towards higher lags. The pattern of decay is obviously reflected in the behaviour of the time series generated from this model, and includes cyclical patterns. Historically, this was a reason for the introduction of these models.

3.4.1 The Yule - Walker equations relating the coefficients and acf of the AR($p$) model.

The coefficients ${\varphi }_1\mathrm{\dots }{\varphi }_p$ determine, and are determined by, the acf values ${\rho }_{x,1}\mathrm{\dots }{\rho }_{x,p}$, by the equations

$${\rho }_{x,i}={\varphi }_1{\rho }_{x,i-1}+\mathrm{\cdots }+{\varphi }_p{\rho }_{x,i-p}\text{ for }i=1\mathrm{\dots }p.$$

Displayed as a set of linear equations for ${\varphi }_1\mathrm{\dots }{\varphi }_p$ these are

where for simplicity the subscript $x$ of ${\rho }_{x,i}$ has been omitted. We estimate the parameters ${\varphi }_1,\mathrm{\dots }{\varphi }_p$ by replacing the acf $\left\{{\rho }_k\right\}$ by sample acf $\left\{r_k\right\}$.

These equations are completed by one further useful relationship:
3.4.3

$$\frac{{\sigma }_e^2}{{\sigma }_x^2}=\left(1-{\varphi }_1{\rho }_{x,1}-{\varphi }_2{\rho }_{x,2}-\mathrm{\cdots }-{\varphi }_p{\rho }_{x,p}\right)$$

Figure 12 shows a good example of a time series that resembles a first order autoregressive process. It is a series of temperature measurements modified by wind speed and other variables for use in predicting energy demands. The sample properties, also shown in Figure 12, are typical of this process.

The model fitted to this series was $x_t={\mu }_x+0.77\left(x_{t-1}-{\mu }_x\right)+e_t$. The properties of the model, presented in the next figures, correspond very well to the sample properties of the series.

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

3.4.4 Example. The AR($2$) model:

$$x_t={\mu }_x+{\varphi }_1\left(x_{t-1}-{\mu }_x\right)+{\varphi }_2\left(x_{t-2}-{\mu }_x\right)+e_t.$$

This process has the capacity to represent a series with an irregular cycle and was put forward by Yule to model the sunspot cycle. The motivation was a discrete time version of the second order differential equation for a damped pendulum, including random fluctuations. For $p=2$ the Yule-Walker equations give:
3.4.5

If the model parameters ${\varphi }_1,{\varphi }_2$ are given, these equations supply ${\rho }_{x,1},{\rho }_{x,2}$. Then 3.4.2 generates successively ${\rho }_{x,3},{\rho }_{x,4}\mathrm{\dots }$ using ${\rho }_{x,i}={\varphi }_1{\rho }_{x,i-1}+{\varphi }_2{\rho }_{x,i-2}$.

Figure 14 shows a series which appears to follow an AR(2) model. It is an indicator of past ocean surface temperatures, obtained by measuring oxygen isotope ratios in a core from the ocean sediment. The time scale is in thousands of years. The series has an irregular cyclical appearance, reflected in its sample autocorrelations.

The coefficients of the fitted AR(2) model were found to be ${\varphi }_1=1.09$ and ${\varphi }_2=-0.42$. The next five figures show the properties of the fitted model, which are similar to those of the sample properties.

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

We began the previous section by saying that the autocorrelation properties of autoregressive models are generally quite distinct from those of moving average models. However we also have

3.5.1 The characteristic property of the AR($p$) model is that its partial autocorrelation function has a ‘cut-off’ at lag $p$, i.e.

$${\rho |}_{x,k}=0\text{ for }k>p.$$

This is because conditioning $x_t$ upon $x_{t-1},x_{t-2},\mathrm{\dots },x_{t-k+1}$ for any $k>p$ includes $x_{t-1},x_{t-2},\mathrm{\dots },x_{t-p}$. If $p$ lagged terms have been included in the predictor, the error $e_t$ is white noise and uncorrelated with all terms $x_{t-k}$ at higher lags than $p$. Including these further terms cannot improve the predictor; this implies that the partial correlation at these lags is zero.

Again there is a converse result, that any stationary time series satisfying the property 3.5.1 may be represented by the AR($p$) model 3.3.4 with $e_t$ being a white noise series, and $\varphi \left(B\right)$ satisfying the stationarity condition.

The characteristic property is recognised by inspection of the sample pacf of the data series. The sampling properties of ${r|}_{x,k}$ are that they are approximately unbiased and consistent estimates of ${\rho |}_{x,k}$ for any fixed lag as the sample size $n$ tends to infinity. When the series $x_t$ follows the AR($p$) model the following sampling property holds:

3.5.2

$$SE\left({r|}_{x,k}\right)\approx \frac{1}{\sqrt{n}}\text{ for }k>p.$$

It is therefore usual to plot the sample pacf on a graph with approximate $95\%$ error limits for zero pacf values of $\pm 2/\sqrt{n}$ about the horizontal axis. We are interested in the case when a small number of values lie clearly outside the limits at low lags. The highest lag with a significant pacf value then indicates the order $p$.

Remark. Unlike the acf which may extend to high lags before decaying effectively to zero, the pacf for most practically occurring series does not extend to very high lags. This is because

$$\text{Var}\left(x_t|x_{t-1},x_{t-2},\mathrm{\dots }{x}_{t-k}\right)=\text{Var}\left(x_t\right)\left(1-{\rho |}_{x,1}^2\right)\left(1-{\rho |}_{x,2}^2\right)\mathrm{\dots }\left(1-{\rho |}_{x,k}^2\right)$$

measures the ability to predict the series from past values, which is usually rather limited. The first few pacf values usually ‘take up’ most of this. So in practice series which do not in actual fact follow an AR model may be well approximated by one of modest order. The true pacf may not actually ‘cut-off’ at some lag, but may just become small compared with the SE of the sample value.

3.6.1 Definition. A time series follows an autoregressive moving average model of orders p and q, an ARMA($p,q$) model, if it is generated by the process

$$x_t={\mu }_x+{\varphi }_1\left(x_{t-1}-{\mu }_x\right)+{\varphi }_2\left(x_{t-2}-{\mu }_x\right)+\mathrm{\cdots }+{\varphi }_p\left(x_{t-p}-{\mu }_x\right)+e_t-{\theta }_1{e}_{t-1}-\mathrm{\cdots }-{\theta }_q{e}_{t-q}$$

where $e_t$ is a white noise series. Using operator notation the model may be written

$$\varphi \left(B\right)\left(x_t-{\mu }_x\right)=\theta \left(B\right)e_t$$

where the operators $\varphi \left(B\right)$ and $\theta \left(B\right)$ are as defined for the pure AR and MA models.

We again require conditions

• •

  for stationarity $\varphi \left(B\right)\ne 0$ for $\left|B\right|\le 1$

• •

  for invertibility $\theta \left(B\right)\ne 0$ for $\left|B\right|\le 1$.

The first ensures that we can express $x_t$ as an infinite moving average form:
3.6.2

$$x_t={\mu }_x+\varphi {\left(B\right)}^{-1}\theta \left(B\right)e_t={\mu }_x+\psi \left(B\right)e_t=e_t+{\psi }_1{e}_{t-1}+{\psi }_2{e}_{t-2}+\mathrm{\cdots }$$

and the second that we can recover $e_t$ from $x_t$ by

$$e_t=\theta {\left(B\right)}^{-1}\varphi \left(B\right)\left(x_t-{\mu }_x\right)=\pi \left(B\right)\left(x_t-{\mu }_x\right)=\left(x_t-{\mu }_x\right)-{\pi }_1\left(x_{t-1}-{\mu }_x\right)-{\pi }_2\left(x_{t-2}-{\mu }_x\right)-\mathrm{\cdots }$$

which is simply rearranged to give the infinite autoregressive form
3.6.3

$$x_t={\mu }_x+{\pi }_1\left(x_{t-1}-{\mu }_x\right)+{\pi }_2\left(x_{t-2}-{\mu }_x\right)+\mathrm{\cdots }+e_t.$$

Remark. The existence of these infinite forms means that in theory an ARMA process can be approximated with arbitrary accuracy by either a pure MA or AR model of finite order. Given a finite data series this can lead to some uncertainty as to which model to use.

3.6.4 Example. The ARMA(1,1) model as defined in 3.6.3 has

$$\psi \left(B\right)=\frac{1-\theta B}{1-\varphi B}=\left(1-\theta B\right)\left(1+\varphi B+{\varphi }^2{B}^2+\mathrm{\cdots }\right)=1+\left(\varphi -\theta \right)B+\left(\varphi -\theta \right)\varphi B^2+\mathrm{\cdots }$$

so that ${\psi }_i=\left(\varphi -\theta \right){\varphi }^{i-1}$. Similarly

$$\pi \left(B\right)=\frac{1-\varphi B}{1-\theta B}=\left(1-\varphi B\right)\left(1+\theta B+{\theta }^2{B}^2+\mathrm{\cdots }\right)=1-\left(\varphi -\theta \right)B-\left(\varphi -\theta \right)\theta B^2+\mathrm{\cdots }$$

so that ${\pi }_i=\left(\varphi -\theta \right){\theta }^{i-1}$.

Remark. Note that in this example if $\varphi =\theta$ then simply $x_t={\mu }_x+e_t$, because $\left(1-\varphi B\right)$ and $\left(1-\theta B\right)$ cancel. In general we require that $\varphi \left(B\right)$ and $\theta \left(B\right)$ have no cancelling factors, else the model would not be distinguishable from a simpler one. The problem can come with model estimation from finite data series when near-cancelling factors can occur. The signal plus noise model can give this situation if the noise almost swamps the signal.

Figure 16 in the next page shows a series of central England average temperatures, that plausibly follows such a model. The coefficients of the fitted ARMA(1,1) model were found to be $\varphi =0.849$ and $\theta =0.655$.

The properties of the fitted model, shown in the next six figures, are similar to those of the sample properties. The last plot shows a smoothed pattern of the series, estimating the underlying level which might be represented by an AR(1) process.

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

The first point to make is that in general neither of the characteristic properties of the MA or AR model holds for ARMA models. This is a negative statement but nevertheless of value in suggesting that a mixed ARMA model is needed.

A second point is that even if it appears that a pure AR or MA model is appropriate from recognition of one characteristic property or the other, a model can sometimes be improved by adding a MA term to a pure AR model or vice versa. This is because inspection of the sample acf and pacf, though providing an indication of the basic model, may ‘lose’ other indications in the sampling fluctuations.

The main concern though is that there may be no simple guide as to the choice of the model orders which are needed in order to estimate the model parameters. The acf properties of the model are considered first. These may be derived from the infinite moving average form 3.6.5 by letting $q=\mathrm{\infty }$ in 3.2.1:
3.7.1

$${\gamma }_{x,k}={\sigma }_e^2\left({\psi }_0{\psi }_k+{\psi }_1{\psi }_{k+1}+\mathrm{\cdots }\right)={\sigma }_e^2\left(\sum _{i=0}^{\mathrm{\infty }}{\psi }_i{\psi }_{k+i}\right)\text{ for }k>0.$$

These infinite sums can be avoided by solving appropriate equations similar to the Yule-Walker equations, although the infinite moving average is initially useful for deriving:
3.7.2

$$\text{Cov}(e_t,x_{t-i})=0\text{ for }i>0,\text{Cov}(x_t,e_t)={\sigma }_e^2,\text{Cov}(x_t,e_{t-i})={\psi }_i{\sigma }_e^2\text{ for }i>0.$$

3.7.3 Example. The ARMA($1,1$) model has variance and autocorrelations given by

$${\sigma }_x^2={\sigma }_e^2\frac{\left(1-2\theta \varphi +{\theta }^2\right)}{\left(1-{\varphi }^2\right)},{\rho }_{x,k}=\frac{\left(\varphi -\theta \right)\left(1-\theta \varphi \right)}{\left(1-2\theta \varphi +{\theta }^2\right)}{\varphi }^{k-1}\text{ for }k>0.$$

Note that the pattern of autocorrelations is similar to that for the AR($1$) in the sense that both exhibit geometric decay of the autocorrelation as a function of the lag $k$.

Proof. We use $\text{Cov}(x_t,e_{t-1})={\psi }_1{\sigma }_e^2=\left(\varphi -\theta \right){\sigma }_e^2$ to derive the first of the equations:

The first two can be solved for ${\sigma }_x^2$ and ${\gamma }_{x,1}$ as stated, and the third shows that ${\gamma }_{x,k}={\varphi }^{k-1}{\gamma }_{x,1}$.

For the general ARMA model the decay pattern of the acf is similar to that for the AR model with the same autoregressive part, but it is modified by the presence of the moving average part to allow greater flexibility of form. The pattern is determined by the equation relating successive acf values for the general ARMA model:
3.7.4

$${\rho }_{x,k}={\varphi }_1{\rho }_{x,k-1}+\mathrm{\cdots }+{\varphi }_p{\rho }_{k-p}\mathit{\hspace{1em}}\text{ for }k>q$$

which derives directly from the model:

This pattern of decay can help to indicate the appropriate AR order $p$. For example cyclical delay indicates that $p$ should be at least $2$. The modification at low lags to the smooth pattern of decay associated with the pure AR model can also help to indicate the MA order $q$, as the ARMA($1,1$) example illustrates. For the ARMA($1,2$) model the geometric decay of the acf starts from lag $2$.

The other statistical tools available for selecting the orders of an ARMA model are the sample pacf and sample spectrum. The sample pacf is not of great value for either the MA or ARMA model, but its pattern is generally similar to that of the ${\pi }_i$ coefficients which appear in the infinite AR representation of these models. For the ARMA($1,1$) example this pattern is approximately a geometric decay by the factor $\theta$ from the lag $1$ value of the pacf.

The remaining sections in this Chapter (3.8–3.12) are not examinable but will be needed for your project. We will discuss these sections after covering Chapters 4–5.

A suitable model for a smooth cycle is
3.8.1

$$y_t=c+R\cos 2\pi \left(ft+g\right)+e_t$$

where

• •

  The Amplitude of the cycle is $R>0$

• •

  The Frequency in cycles per unit time (that is over $t\in [0,1]$) is $f$ where .

• •

  The Phase, that is the fraction of a cycle by which it lags, is $g$, .

We restrict $R\ge 0$, $f\ge 0$ and .

We shall sometimes write $R\cos 2\pi \left(ft+g\right)=R\cos \left(\omega t+\varphi \right)$ using the angular frequency $\omega =2\pi f$ of radians per unit time and angular phase $\varphi =2\pi g$ - the notation is more convenient at times.

Rather than estimate $R$ and $\varphi$ directly we use

$$\cos \left(\omega t+\varphi \right)=\cos \omega t\cos \varphi -\sin \omega t\sin \varphi$$

to rewrite the model as
3.8.2

$$y_t=c+A\cos \omega t+B\sin \omega t+e_t$$

where $A=R\cos \varphi$ and $B=-R\sin \varphi .$ So if $\omega$ is known this is a linear regression model with coefficients $A$ and $B$.

Sometimes the fit to a smooth cycle is improved by adding further cycles known as harmonics because their frequencies are multiples of the fundamental frequency $\omega$, e.g.

$$y_t=c+A_1\cos \omega t+B_1\sin \omega t+A_2\cos 2\omega t+B_2\sin 2\omega t+A_3\cos 3\omega t+B_3\sin 3\omega t+e_t$$

Thus cycles which are asymmetric, with a slow rise and rapid fall, or sharp peaks with flat troughs, may be well fitted by including sufficient harmonics provided the pattern repeats itself periodically.