5 Forecasting (prediction, extrapolation) using ARIMA models.

We have already met forecasting in one form or another. Firstly we extrapolated the deterministic components of regression models into the future, for trend, seasonal and cyclical models. Therefore, in the model $y_t=x_t^{{}^{\prime }}\beta +{ϵ}_t$ where $\left\{{ϵ}_t\right\}$ are uncorrelated with zero mean, the forecast for $y_{t+k}=x_{t+k}^{{}^{\prime }}\beta +{ϵ}_t$ based on $y_1,\mathrm{\dots },y_t$ is $x_{t+k}^{{}^{\prime }}\widehat{\beta }$ where $\widehat{\beta }$ is the estimate of $\beta$ based on $y_1,\mathrm{\dots },y_t$.

This is still appropriate if such components are used in a model for a time series which includes ARMA error structure. But now the ARMA part also needs to be extrapolated. For uncorrelated errors the extrapolations would be simply zero. But when the errors are correlated the expected future values are not all zero.

We also came across the concept of prediction when introducing autoregressive models. These again exploit regression - lagged regression on past values to predict the future.

However the approach in this chapter is model based. Forecasting is simply asking the question about where the process (that has been modelled) is going in the future, given that you have observed it up to the present. We investigate this now for models that can be expressed in the form

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

To understand prediction of future values of the process we will now investigate how future values $x_{t+1},x_{t+2},\mathrm{\dots }$ depend in part upon the past (and present) quantities $x_t,x_{t-1},\mathrm{\dots }$ and $e_t,e_{t-1},\mathrm{\dots }$ that are known at time $t$; and in part upon the future innovations $e_{t+1},e_{t+2},\mathrm{\dots }$. We first look at particular simple examples. We omit ${\mu }_x$ for simplicity. The method we use is to write the model down with $t$ replaced by $t+1$ then $t+2$ …. If the model has terms on its right hand side which include any values in $x_{t+1},x_{t+2},\mathrm{\dots }$, we successively substitute for each of these so no such terms remain; only terms in $e_{t+1},e_{t+2},\mathrm{\dots }$ remain on the right, besides those in $x_t,x_{t-1},\mathrm{\dots }$ and $e_t,e_{t-1},\mathrm{\dots }$

5.2.1 Example: the MA(2) model:

So for $k>2$, $x_{t+k}$ depends only on future innovations; not at all on the past.

5.2.2 Example: the AR(1) model:

Thus the dependence of $x_{t+k}$ on the past is given by the coefficient ${\varphi }^k$ of $x_t$, decaying geometrically.

Now recall the infinite moving average representation of the AR(1) model, which after replacing $t$ by $t+k$ becomes

	$x_{t+k}$	$=$	${\left(1-\varphi B\right)}^{-1}{e}_{t+k}=\left(1+\varphi B+{\varphi }^2{B}^2+{\varphi }^3{B}^3+\mathrm{\cdots }\right)e_{t+k}$
		$=$	$e_{t+k}+\varphi e_{t+k-1}+{\varphi }^2{e}_{t+k-2}+{\varphi }^3{e}_{t+k-3}\mathrm{\dots },$

and note that this gives the dependence of $x_{t+k}$ upon $e_{t+1},e_{t+2},\mathrm{\dots }$ that we observed in the previous equation.

5.2.3 Example: the IMA($1,1$) model:

	$x_{t+1}$	$=$	$x_t+e_{t+1}-\theta e_t$
	$x_{t+2}=x_{t+1}+e_{t+2}-\theta e_{t+1}$	$=$	$x_t+e_{t+2}+\left(1-\theta \right)e_{t+1}-\theta e_t$
	$x_{t+3}=x_{t+2}+e_{t+3}-\theta e_{t+2}$	$=$	$x_t+e_{t+3}+\left(1-\theta \right)e_{t+2}+\left(1-\theta \right)e_{t+1}-\theta e_t$
	$x_{t+4}=x_{t+3}+e_{t+4}-\theta e_{t+3}$	$=$	$x_t+e_{t+4}+\left(1-\theta \right)e_{t+3}+\left(1-\theta \right)e_{t+2}+\left(1-\theta \right)e_{t+1}-\theta e_t$

Note now that the infinite moving average representation of $x_{t+k}$ for the IMA(1,1) model is formally

5.2.4

We say ‘formally’ because this sum does not converge, but it does supply the correct coefficients by which the previous equation shows how the innovations $e_{t+1},e_{t+2},\mathrm{\dots }$ affect the future value $x_{t+k}$. In this sense we may still write

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

where for an ARIMA model

$$x_t=\frac{\theta \left(B\right)}{{\nabla }^d\varphi \left(B\right)}{e}_t\mathit{\hspace{1em}\hspace{1em} }\Rightarrow \mathit{\hspace{1em}\hspace{1em} }\psi \left(B\right)=\frac{\theta \left(B\right)}{{\nabla }^d\varphi \left(B\right)}.$$

Suppose we have observations of $x_t$ up to and including time $t=n$ and, through the process of estimating the model, values of $e_t$ also up to and including time $t=n$. We want to forecast future values $x_{n+k}$ for $k=1,2,3,\mathrm{\dots }$, which we call ${\widehat{x}}_{n,k}$.

The general form of the infinite MA for $x_{n+k}$ may be written

This shows explicitly that the part of $x_{n+k}$ which depends on future unknown innovations, is

$$e_{n,k}=e_{n+k}+{\psi }_1{e}_{n+k-1}+\mathrm{\cdots }+{\psi }_{k-1}{e}_{n+1}$$

and the remaining part of $x_{n+k}$ is known. In other words

$$x_{n+k}=\sum _{j=0}^{\mathrm{\infty }}{\psi }_j{e}_{n+k-j}=\sum _{j=k}^{\mathrm{\infty }}{\psi }_j{e}_{n+k-j}+\sum _{j=0}^{k-1}{\psi }_j{e}_{n+k-j}={\widehat{x}}_{n,k}+e_{n,k}.$$

The forecast error variance is therefore:

	$V_k=\text{Var}\left(e_{n,k}\right)$	$=$	$\text{Var}\left(e_{n+k}+{\psi }_1{e}_{n+k-1}+\mathrm{\cdots }+{\psi }_{k-1}{e}_{n+1}\right)$
		$=$	${\sigma }_e^2\left(1+{\psi }_1^2+{\psi }_2^2+\mathrm{\cdots }+{\psi }_{k-1}^2\right).$

When forecasts made using ARIMA models are graphed for increasing lead time $k$ from a fixed forecast origin $n$, i.e.

$${\widehat{x}}_{n,k},k=1,2,3,\mathrm{\dots }$$

the resulting graph is called a forecast function. It is usual to graph the forecast function with limits of $\pm 2\sqrt{V_k}$ to indicate the range of possible future paths of the series.

Any two such functions are usually different because they depend on the latest values of the series. However for a given model a similar pattern is usually evident in all forecast functions. To see why, consider the ARMA($p,q$) model with $t=n+k$ where $k>q$, and setting $e_{n+1}=e_{n+2}=\mathrm{\cdots }=0$ to obtain the forecasts:

$$x_{n+k}={\mu }_x+{\varphi }_1\left(x_{n+k-1}-{\mu }_x\right)+\mathrm{\cdots }+{\varphi }_p\left(x_{n+k-p}-{\mu }_x\right)+e_{n+k}-{\theta }_1{e}_{n+k-1}-\mathrm{\cdots }-{\theta }_q{e}_{n+k-q}$$

$$\Rightarrow \mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} }{\widehat{x}}_{n,k}={\mu }_x+{\varphi }_1\left({\widehat{x}}_{n,k-1}-{\mu }_x\right)+\mathrm{\cdots }+{\varphi }_p\left({\widehat{x}}_{n,k-p}-{\mu }_x\right)$$

because all terms $e_{n+k},e_{n+k-1},\mathrm{\dots }{e}_{n+k-q}$ have been set to zero. This may also be expressed:

$$\varphi \left(B\right)\left({\widehat{x}}_{n,k}-{\mu }_x\right)=0\mathit{\hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{1em}}\text{for}\mathit{\hspace{1em}\hspace{1em}}k>q.$$

This is the same recurrence relationship which generated successive autocorrelations ${\rho }_{x,k}$ for the ARMA model. If $p=1$ then a geometric decay to the mean of the series results. If $p=2$ and $\varphi \left(B\right)$ has complex factors the values will follow a damped cycle as they decay to ${\mu }_x$.

Now consider ARIMA models. For example the IMA($1,1$) model with drift:

$$x_{n+k}=x_{n+k-1}+\mu +e_{n+k}-\theta e_{n+k-1}$$

$$\Rightarrow \mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} }{\widehat{x}}_{n,k}={\widehat{x}}_{n,k-1}+\mu \mathit{\hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{1em}}\text{for}\mathit{\hspace{1em}\hspace{1em}}k>\mathrm{\hspace{0.17em}1}.$$

so that the forecast function follows a trend line with slope $\mu$. For this model we can show that $V_k={\sigma }_e^2\left\{1+\left(k-1\right)\left(1-{\theta }^2\right)\right\}$ so that the limits gradually widen out around the forecast function.

When new observations are made, for example of $x_{n+1},x_{n+2},\mathrm{\dots }$ the new innovations $e_{n+1},e_{n+2},\mathrm{\dots }$ can be regenerated and at each time point $t=n+1,\mathrm{\dots }$ in succession forecasts can be generated by the means described. It is however possible to develop equations to update the forecast function in a convenient way. The EWMA updating equation is a good example. The idea is to update from ${\widehat{x}}_{n,k}$ to ${\widehat{x}}_{n+1,k}$. From the expression given earlier for forecast errors we can get

$${\widehat{x}}_{n+1,k}={\widehat{x}}_{n,k+1}+{\psi }_k{e}_{n+1}.$$

By expressing $e_{n+1}=x_{n+1}-{\widehat{x}}_{n,1}$ and ${\widehat{x}}_{n,k+1}$ in terms of ${\widehat{x}}_{n,k},\mathrm{\dots }$ this can be re-arranged to give the updating formula. For the IMA($1,1$) model we start with ${\psi }_1=\left(1-\theta \right)$ so that

$${\widehat{x}}_{n+1,1}={\widehat{x}}_{n,2}+\left(1-\theta \right)e_{n+1}$$

and using ${\widehat{x}}_{n,2}={\widehat{x}}_{n,1}$ and $e_{n+1}=x_{n+1}-{\widehat{x}}_{n,1}$ gives

$${\widehat{x}}_{n+1,1}={\widehat{x}}_{n,1}+\left(1-\theta \right)\left(x_{n+1}-{\widehat{x}}_{n,1}\right)=\left(1-\theta \right)x_{n+1}+\theta {\widehat{x}}_{n,1}.$$

We conclude with some illustrations of forecast functions in Figure 26. In all of these the model is fitted to the data points before the forecasts start. The first of the four in Figure 26, shows forecasts of the weekly Financial Transaction series. The logarithms of the series are modelled by an IMA(1,1) model, but the forecasts and limits are transformed back to the original scale (and therefore you would not see constant forecast). The second shows forecasts of the daily temperatures, using an AR(1) model. The reversion of the forecast function to the mean level can be seen, as the error limits widen.

The third and fourth are forecasts of the annual sunspot series. An AR(2) model is used in the third, and an AR(8) in the fourth. The forecasts from the AR(2) model are seen to dampen down more quickly than those from the AR(8) model, which better captures the dynamics of the series. The error limits are actually only 90%, i.e. plus or minus 1.65 standard deviations. The errors for increasing lead-time are generally strongly correlated, so if the actual series values lie outside the limit at some point, it is quite likely that some of those close in time, will also do so.

The following four figures show forecasts of selected series using respectively an IMA(1,1), AR(1), AR(2) and AR(8) models.

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Forecasting is one of the most challenging applications of univariate time series analysis. There are many other applications, most of them being smoothing or extracting components of a series, for which model based approaches have many advantages. Classification of series is also important, with applications in medical diagnosis and seismic analysis. Multivariate time series opens up a new field of modelling the dependence between time series, which also has widespread applications.