Autoregressive model

In statistics and signal processing, an autoregressive (AR) model is a type of random process which is often used to model and predict various types of natural phenomena. The autoregressive model is one of a group of linear prediction formulas that attempt to predict an output of a system based on the previous outputs.

Definition

The notation AR(p) indicates an autoregressive model of order p. The AR(p) model is defined as

$X_{t}=c+\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}\,$

where $\varphi _{1},\ldots ,\varphi _{p}$ are the parameters of the model, $c$ is a constant (often omitted for simplicity) and $\varepsilon _{t}$ is white noise.

An autoregressive model can thus be viewed as the output of an all-pole infinite impulse response filter whose input is white noise.

Some constraints are necessary on the values of the parameters of this model in order that the model remains wide-sense stationary. For example, processes in the AR(1) model with |φ₁| ≥ 1 are not stationary. More generally, for an AR(p) model to be wide-sense stationary, the roots of the polynomial $\textstyle z^{p}-\sum _{i=1}^{p}\varphi _{i}z^{p-i}$ must lie within the unit circle, i.e., each root $z_{i}$ must satisfy $|z_{i}|<1$.

Characteristic polynomial

The autocorrelation function of an AR(p) process can be expressed as

$\rho (\tau )=\sum _{k=1}^{p}a_{k}y_{k}^{-|\tau |},$

where $y_{k}$ are the roots of the polynomial

$\phi (B)=1-\sum _{k=1}^{p}\varphi _{k}B^{k}.$

The autocorrelation function of an AR(p) process is a sum of decaying exponentials.

• Each real root contributes a component to the autocorrelation function that decays exponentially.
• Similarly, each pair of complex conjugate roots contributes an exponentially damped oscillation.

Graphs of AR(p) processes

AR(0), AR(1), and AR(2) processes with white noise

The simplest AR process is AR(0), which has no dependence between the terms. Only the error/innovation/noise term contributes to the output of the process, so in the figure, AR(0) corresponds to white noise.

For an AR(1) process with a positive $\varphi$, only the previous term in the process and the noise term contribute to the output. If $\varphi$ is close to 0, then the process still looks like white noise, but as $\varphi$ approaches 1, the output gets a larger contribution from the previous term relative to the noise. This results in a "smoothing" or integration of the output, similar to a low pass filter.

For an AR(2) process, the previous two terms and the noise term contribute to the output. If both $\varphi _{1}$ and $\varphi _{2}$ are positive, the output will resemble a low pass filter, with the high frequency part of the noise decreased. If $\varphi _{1}$ is positive while $\varphi _{2}$ is negative, then the process favors changes in sign between terms of the process. The output oscillates.

Example: An AR(1)-process

An AR(1)-process is given by:

$X_{t}=c+\varphi X_{t-1}+\varepsilon _{t}\,$

where $\varepsilon _{t}$ is a white noise process with zero mean and variance $\sigma _{\varepsilon }^{2}$. (Note: The subscript on $\varphi _{1}$ has been dropped.) The process is wide-sense stationary if $|\varphi |<1$ since it is obtained as the output of a stable filter whose input is white noise. (If $\varphi =1$ then $X_{t}$ has infinite variance, and is therefore not wide sense stationary.) Consequently, assuming $|\varphi |<1$, the mean $\operatorname {E} (X_{t})$ is identical for all values of t. If the mean is denoted by $\mu$, it follows from

$\operatorname {E} (X_{t})=\operatorname {E} (c)+\varphi \operatorname {E} (X_{t-1})+\operatorname {E} (\varepsilon _{t}),$

that

$\mu =c+\varphi \mu +0,$

and hence

$\mu ={\frac {c}{1-\varphi }}.$

In particular, if $c=0$, then the mean is 0.

The variance is

${\textrm {var}}(X_{t})=\operatorname {E} (X_{t}^{2})-\mu ^{2}={\frac {\sigma _{\varepsilon }^{2}}{1-\varphi ^{2}}},$

where $\sigma _{\varepsilon }$ is the standard deviation of $\varepsilon _{t}$. This can be shown by noting that

${\textrm {var}}(X_{t})=\varphi ^{2}{\textrm {var}}(X_{t-1})+\sigma _{\varepsilon }^{2},$

and then by noticing that the quantity above is a stable fixed point of this relation.

The autocovariance is given by

$B_{n}=\operatorname {E} (X_{t+n}X_{t})-\mu ^{2}={\frac {\sigma _{\varepsilon }^{2}}{1-\varphi ^{2}}}\,\,\varphi ^{|n|}.$

It can be seen that the autocovariance function decays with a decay time (also called time constant) of $\tau =-1/\ln(\varphi )$ [to see this, write $B_{n}=K\phi ^{|n|}$ where $K$ is independent of $n$. Then note that $\phi ^{|n|}=e^{|n|\ln \phi }$ and match this to the exponential decay law $e^{-n/\tau }$].

The spectral density function is the Fourier transform of the autocovariance function. In discrete terms this will be the discrete-time Fourier transform:

$\Phi (\omega )={\frac {1}{\sqrt {2\pi }}}\,\sum _{n=-\infty }^{\infty }B_{n}e^{-i\omega n}={\frac {1}{\sqrt {2\pi }}}\,\left({\frac {\sigma _{\varepsilon }^{2}}{1+\varphi ^{2}-2\varphi \cos(\omega )}}\right).$

This expression is periodic due to the discrete nature of the $X_{j}$, which is manifested as the cosine term in the denominator. If we assume that the sampling time ($\Delta t=1$) is much smaller than the decay time ($\tau$), then we can use a continuum approximation to $B_{n}$:

$B(t)\approx {\frac {\sigma _{\varepsilon }^{2}}{1-\varphi ^{2}}}\,\,\varphi ^{|t|}$

which yields a Lorentzian profile for the spectral density:

$\Phi (\omega )={\frac {1}{\sqrt {2\pi }}}\,{\frac {\sigma _{\varepsilon }^{2}}{1-\varphi ^{2}}}\,{\frac {\gamma }{\pi (\gamma ^{2}+\omega ^{2})}}$

where $\gamma =1/\tau$ is the angular frequency associated with the decay time $\tau$.

An alternative expression for $X_{t}$ can be derived by first substituting $c+\varphi X_{t-2}+\varepsilon _{t-1}$ for $X_{t-1}$ in the defining equation. Continuing this process N times yields

$X_{t}=c\sum _{k=0}^{N-1}\varphi ^{k}+\varphi ^{N}X_{t-N}+\sum _{k=0}^{N-1}\varphi ^{k}\varepsilon _{t-k}.$

For N approaching infinity, $\varphi ^{N}$ will approach zero and:

$X_{t}={\frac {c}{1-\varphi }}+\sum _{k=0}^{\infty }\varphi ^{k}\varepsilon _{t-k}.$

It is seen that $X_{t}$ is white noise convolved with the $\varphi ^{k}$ kernel plus the constant mean. If the white noise $\varepsilon _{t}$ is a Gaussian process then $X_{t}$ is also a Gaussian process. In other cases, the central limit theorem indicates that $X_{t}$ will be approximately normally distributed when $\varphi$ is close to one.

Calculation of the AR parameters

There are many ways to estimate the coefficients, such as the ordinary least squares procedure, method of moments (through Yule Walker equations), or Markov chain Monte Carlo methods.

The AR(p) model is given by the equation

$X_{t}=\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}.\,$

It is based on parameters $\varphi _{i}$ where i = 1, ..., p. There is a direct correspondence between these parameters and the covariance function of the process, and this correspondence can be inverted to determine the parameters from the autocorrelation function (which is itself obtained from the covariances). This is done using the Yule-Walker equations.

Yule-Walker equations

The Yule-Walker equations are the following set of equations.

$\gamma _{m}=\sum _{k=1}^{p}\varphi _{k}\gamma _{m-k}+\sigma _{\varepsilon }^{2}\delta _{m,0},$

where m = 0, ..., p, yielding p + 1 equations. Here $\gamma _{m}$ is the autocovariance function of X_t, $\sigma _{\varepsilon }$ is the standard deviation of the input noise process, and $\delta _{m,0}$ is the Kronecker delta function.

Because the last part of an individual equation is non-zero only if m = 0, the set of equations can be solved by representing the equations for m > 0 in matrix form, thus getting the equation

${\begin{bmatrix}\gamma _{1}\\\gamma _{2}\\\gamma _{3}\\\vdots \\\gamma _{p}\\\end{bmatrix}}={\begin{bmatrix}\gamma _{0}&\gamma _{-1}&\gamma _{-2}&\dots \\\gamma _{1}&\gamma _{0}&\gamma _{-1}&\dots \\\gamma _{2}&\gamma _{1}&\gamma _{0}&\dots \\\vdots &\vdots &\vdots &\ddots \\\gamma _{p-1}&\gamma _{p-2}&\gamma _{p-3}&\dots \\\end{bmatrix}}{\begin{bmatrix}\varphi _{1}\\\varphi _{2}\\\varphi _{3}\\\vdots \\\varphi _{p}\\\end{bmatrix}}$

which can be solved for all $\{\varphi _{m};m=1,2,\cdots ,p\}.$ The remaining equation for m = 0 is

$\gamma _{0}=\sum _{k=1}^{p}\varphi _{k}\gamma _{-k}+\sigma _{\varepsilon }^{2},$

which, once $\{\varphi _{m};m=1,2,\cdots ,p\}$ are known, can be solved for $\sigma _{\varepsilon }^{2}.$

An alternative formulation is in terms of the autocorrelation function. The AR parameters are determined by the first p+1 elements $\rho (\tau )$ of the autocorrelation function. The full autocorrelation function can then be derived by recursively calculating

$\rho (\tau )=\sum _{k=1}^{p}\alpha _{k}\rho (k-\tau )$

Examples for some Low-order AR(p) processes

• p=1
  • $\gamma _{1}=\varphi _{1}\gamma _{0}$
  • Hence $\rho _{1}=\gamma _{1}/\gamma _{0}=\varphi _{1}$
• p=2
  • The Yule-Walker equations for an AR(2) process are

    $\gamma _{1}=\varphi _{1}\gamma _{0}+\varphi _{2}\gamma _{-1}$

    $\gamma _{2}=\varphi _{1}\gamma _{1}+\varphi _{2}\gamma _{0}$

    • Remember that $\gamma _{-k}=\gamma _{k}$
    • Using the first equation yields $\rho _{1}=\gamma _{1}/\gamma _{0}={\frac {\varphi _{1}}{1-\varphi _{2}}}$
    • Using the recursion formula yields $\rho _{2}=\gamma _{2}/\gamma _{0}={\frac {\varphi _{1}^{2}-\varphi _{2}^{2}+\varphi _{2}}{1-\varphi _{2}}}$

Estimation of AR parameters

The above equations (the Yule-Walker equations) provide several routes to estimating the parameters of an AR(p) model, by replacing the theoretical covariances with estimated values. Some of these variants can be described as follows:

• Estimation of autocovariances or autocorrelations. Here each of these terms is estimated separately, using conventional estimates. There are different ways of doing this and the choice between these affects the properties of the estimation scheme. For example, negative estimates of the variance can be produced by some choices.

• Formulation as a least squares regression problem in which an ordinary least squares prediction problem is constructed, basing prediction of values of X_t on the p previous values of the same series. This can be thought of as a forward-prediction scheme. The normal equations for this problem can be seen to correspond to an approximation of the matrix form of the Yule-Walker equations in which each appearance of an autocovariance of the same lag is replaced by a slightly different estimate.

• Formulation as an extended form of ordinary least squares prediction problem. Here two sets of prediction equations are combined into a single estimation scheme and a single set of normal equations. One set is the set of forward-prediction equations and the other is a corresponding set of backward prediction equations, relating to the backward representation of the AR model:

$X_{t}=c+\sum _{i=1}^{p}\varphi _{i}X_{t+i}+\varepsilon _{t}^{*}\,.$

Here predicted of values of X_t would be based on the p future values of the same series. This way of estimating the AR parameters is due to Burg, and call the Burg method: Burg and later authors called these particular estimates "maximum entropy estimates", but the reasoning behind this applies to the use of any set of estimated AR parameters. Compared to the estimation scheme using only the forward prediction equations, different estimates of the autocovariances are produced, and the estimates have different stability properties. Burg estimates are particularly associated with maximum entropy spectral estimation.

Other possible approaches to estimation include maximum likelihood estimation. Two distinct variants of maximum likelihood are available: in one (broadly equivalent to the forward prediction least squares scheme) the likelihood function considered is that corrresponding to the conditional distribution of later values in the series given the initial p values in the series; in the second, the likelihood function considered is that corrresponding to the unconditional joint distribution of all the values in the observed series. Substantial differences in the results of these approaches can occur if the observed series is short, or if the process is close to non-stationarity.

Spectrum

The power spectral density of an AR(p) process with noise variance $Var(Z_{t})=\sigma _{Z}^{2}$ is

$S(f)={\frac {\sigma _{Z}^{2}}{|1-\sum _{k=1}^{p}\varphi _{k}e^{-2\pi ikf}|^{2}}}.$

AR(0)

For white noise (AR(0))

$S(f)=\sigma _{Z}^{2}.$

AR(1)

For AR(1)

$S(f)={\frac {\sigma _{Z}^{2}}{|1-\varphi _{1}e^{-2\pi ikf}|^{2}}}={\frac {\sigma _{Z}^{2}}{1+\varphi _{1}^{2}-2\varphi _{1}cos{2\pi f}}}$

• If $\varphi _{1}>0$ there is a single spectral peak at f=0, often referred to as red noise. As $\varphi _{1}$ becomes nearer 1, there is stronger power at low frequencies, i.e. larger time lags.
• If $\varphi _{1}<0$ there is a minimum at f=0, often referred to as blue noise

AR(2)

AR(2) processes can be split into three groups depending on the characteristics of their roots:

$z_{1},z_{2}={\frac {1}{2}}\left(\varphi _{1}\pm {\sqrt {\varphi _{1}^{2}+4\varphi _{2}}}\right)$

• When $\varphi _{1}^{2}+4\varphi _{2}<0$, the process has a pair of complex-conjugate roots, creating a mid-frequency peak at:

$f^{*}={\frac {1}{2\pi }}cos^{-1}\left({\frac {(1-\varphi _{2})^{2}(4\varphi _{2}+\varphi _{1}^{2})}{4\varphi _{2}}}\right)$

Otherwise the process has real roots, and:

• When $\varphi _{1}>0$ it acts as a low-pass filter on the white noise with a spectral peak at $f=0$
• When $\varphi _{1}<0$ it acts as a high-pass filter on the white noise with a spectral peak at $f=1/2$.

The process is stable when the roots are within the unit circle, or equivalently when the coefficients are in the triangle $-1\leq \varphi _{2}\leq 1-|\varphi _{1}|$.

The full PSD function can be expressed in real form as:

$S(f)={\frac {\sigma _{Z}^{2}}{1+\varphi _{1}^{2}+\varphi _{2}^{2}-2\varphi _{1}(1-\varphi _{2})\cos(2\pi f)-2\varphi _{2}\cos(4\pi f)}}$

Implementations in statistics packages

• R, the stats package includes an ar function.
• Matlab and Octave: the TSA toolbox contains several estimation functions for uni-variate, multivariate and adaptive autoregressive models.

n-step-ahead forecasting

Once the parameters of the autoregression

$X_{t}=c+\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}\,$

have been estimated, the autoregression can be used to forecast an arbitrary number of periods into the future. First use t to refer to the first period for which data is not yet available; substitute the known prior values X_t-i for i=1, ..., p into the autoregressive equation while setting the error term $\varepsilon _{t}$ equal to zero (because we forecast X_t to equal its expected value, and the expected value of the unobserved error term is zero). The output of the autoregressive equation is the forecast for the first unobserved period. Next, use t to refer to the next period for which data is not yet available; again the autoregressive equation is used to make the forecast, with one difference: the value of X one period prior to the one now being forecast is not known, so its expected value—the predicted value arising from the previous forecasting step—is used instead. Then for future periods the same procedure is used, each time using one more forecast value on the right side of the predictive equation until, after p predictions, all p right-side values are predicted values from prior steps.

Evaluating the quality of forecasts

The predictive performance of the autoregressive model can be assessed as soon as estimation has been done if cross-validation is used. In this approach, some of the initially available data was used for parameter estimation purposes, and some (from available observations later in the data set) was held back for out-of-sample testing. Alternatively, after some time has passed after the parameter estimation was conducted, more data will have become available and predictive performance can be evaluated then using the new data.

In either case, there are two aspects of predictive performance that can be evaluated: one-step-ahead and n-step-ahead performance. For one-step-ahead performance, the estimated parameters are used in the autoregressive equation along with observed values of X for all periods prior to the one being predicted, and the output of the equation is the one-step-ahead forecast; this procedure is used to obtain forecasts for each of the out-of-sample observations. To evaluate the quality of n-step-ahead forecasts, the forecasting procedure in the previous section is employed to obtain the predictions.

Given a set of predicted values and a corresponding set of actual values for X for various time periods, a common evaluation technique is to use the mean squared prediction error; other measures are also available (see Forecasting#Forecasting accuracy).

The question of how to interpret the measured forecasting accuracy arises—for example, what is a "high" (bad) or a "low" (good) value for the mean squared prediction error? There are two possible points of comparison. First, the forecasting accuracy of an alternative model, estimated under different modeling assumptions or different estimation techniques, can be used for comparison purposes. Second, the out-of-sample accuracy measure can be compared to the same measure computed for the in-sample data points (that were used for parameter estimation) for which enough prior data values are available (that is, dropping the first p data points, for which p prior data points are not available). Since the model was estimated specifically to fit the in-sample points as well as possible, it will usually be the case that the out-of-sample predictive performance will be poorer than the in-sample predictive performance. But if the predictive quality deteriorates out-of-sample by "not very much" (which is not precisely definable), then the forecaster may be satisfied with the performance.

See also

• Moving average model
• Autoregressive moving average model
• Predictive analytics
• Linear predictive coding

References

• Mills, Terence C. (1990) Time Series Techniques for Economists. Cambridge University Press
• Percival, Donald B. and Andrew T. Walden. (1993) Spectral Analysis for Physical Applications. Cambridge University Press
• Pandit, Sudhakar M. and Wu, Shien-Ming. (1983) Time Series and System Analysis with Applications. John Wiley & Sons
• Yule, G. Udny (1927) "On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer's Sunspot Numbers", Philosophical Transactions of the Royal Society of London, Ser. A, Vol. 226, 267–298.]
• Walker, Gilbert (1931) "On Periodicity in Series of Related Terms", Proceedings of the Royal Society of London, Ser. A, Vol. 131, 518–532.

External links

• AutoRegression Analysis (AR) by Paul Bourke