{{Short description|Representation of a type of random process}} {{More footnotes needed|date=March 2011}} In statistics, an '''autoregressive''' ('''AR''') '''model''' is a modelled representation of a type of random process. It can be used to describe time-varying processes from many natural and artificial sources. The model specifies output variables that are dependent linearly on their own previous values on a stochastic basis. The model is in the form of a stochastic difference equation (or recurrence relation) which should not be confused with a differential equation. Together with the moving-average (MA) model, it is a special case and key component of the more general autoregressive–moving-average (ARMA) and autoregressive integrated moving average (ARIMA) models of time series, which have a more complicated stochastic structure; it is also a special case of the vector autoregressive model (VAR), which consists of a system of more than one interlocking stochastic difference equation in more than one evolving random variable. Another important extension is the time-varying autoregressive (TVAR) model, where the autoregressive coefficients are allowed to change over time to model evolving or non-stationary processes. TVAR models are widely applied in cases where the underlying dynamics of the system are not constant, such as in sensors time series modelling,<ref>{{cite journal |last1=Souza |first1=Douglas Baptista de |last2=Leao |first2=Bruno Paes |title=Data Augmentation of Sensor Time Series using Time-varying Autoregressive Processes |journal=Annual Conference of the PHM Society |date=26 October 2023 |volume=15 |issue=1 |doi=10.36001/phmconf.2023.v15i1.3565|doi-access=free }}</ref><ref>{{cite journal |last1=Souza |first1=Douglas Baptista de |last2=Leao |first2=Bruno Paes |title=Data Augmentation of Multivariate Sensor Time Series using Autoregressive Models and Application to Failure Prognostics |journal=Annual Conference of the PHM Society |date=5 November 2024 |volume=16 |issue=1 |doi=10.36001/phmconf.2024.v16i1.4145|arxiv=2410.16419 }}</ref> climate science,<ref>{{cite journal |last1=Diodato |first1=Nazzareno |last2=Di Salvo |first2=Cristina |last3=Bellocchi |first3=Gianni |title=Climate driven generative time-varying model for improved decadal storm power predictions in the Mediterranean |journal=Communications Earth & Environment |date=18 March 2025 |volume=6 |issue=1 |pages=212 |doi=10.1038/s43247-025-02196-2 |bibcode=2025ComEE...6..212D |doi-access=free }}</ref> economics and finance (as econometrics),<ref>{{cite journal |last1=Inayati |first1=Syarifah |last2=Iriawan |first2=Nur |title=Time-Varying Autoregressive Models for Economic Forecasting |journal=Matematika |date=31 December 2024 |pages=131–142 |doi=10.11113/matematika.v40.n3.1654|doi-access=free }}</ref><ref>{{cite journal |last1=Jia |first1=Zhixuan |last2=Li |first2=Wang |last3=Jiang |first3=Yunlong |last4=Liu |first4=Xingshen |title=The Use of Minimization Solvers for Optimizing Time-Varying Autoregressive Models and Their Applications in Finance |journal=Mathematics |date=9 July 2025 |volume=13 |issue=14 |pages=2230 |doi=10.3390/math13142230 |doi-access=free }}</ref> signal processing,<ref>{{cite journal |last1=Baptista de Souza |first1=Douglas |last2=Kuhn |first2=Eduardo Vinicius |last3=Seara |first3=Rui |title=A Time-Varying Autoregressive Model for Characterizing Nonstationary Processes |journal=IEEE Signal Processing Letters |date=January 2019 |volume=26 |issue=1 |pages=134–138 |doi=10.1109/LSP.2018.2880086 |bibcode=2019ISPL...26..134B }}</ref> telecommunications,<ref>{{cite journal |last1=Wang |first1=Shihan |last2=Chen |first2=Tao |last3=Wang |first3=Hongjian |title=IDBD-Based Beamforming Algorithm for Improving the Performance of Phased Array Radar in Nonstationary Environments |journal=Sensors |date=17 March 2023 |volume=23 |issue=6 |pages=3211 |doi=10.3390/s23063211 |pmid=36991922 |pmc=10052024 |bibcode=2023Senso..23.3211W |doi-access=free }}</ref> radar systems,<ref>{{cite journal |last1=Abramovich |first1=Yuri I. |last2=Spencer |first2=Nicholas K. |last3=Turley |first3=Michael D. E. |title=Time-Varying Autoregressive (TVAR) Models for Multiple Radar Observations |journal=IEEE Transactions on Signal Processing |date=April 2007 |volume=55 |issue=4 |pages=1298–1311 |doi=10.1109/TSP.2006.888064 |bibcode=2007ITSP...55.1298A }}</ref> and biological signals.<ref>{{cite book |last1=Gutierrez |first1=D. |last2=Salazar-Varas |first2=R. |chapter=EEG signal classification using time-varying autoregressive models and common spatial patterns |title=2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society |date=August 2011 |pages=6585–6588 |doi=10.1109/IEMBS.2011.6091624 |pmid=22255848 |isbn=978-1-4577-1589-1 }}</ref>
Unlike the moving-average (MA) model, the autoregressive model is not always stationary; non-stationarity can arise either due to the presence of a unit root or due to time-varying model parameters, as in time-varying autoregressive models.
Large language models are called autoregressive, but they are not a classical autoregressive model in this sense because they are not linear.
==Definition==
The notation <math>AR(p)</math> indicates an autoregressive model of order ''p''. The AR(''p'') model is defined as
:<math> X_t = \sum_{i=1}^p \varphi_i X_{t-i} + \varepsilon_t</math>
where <math>\varphi_1, \ldots, \varphi_p</math> are the ''parameters'' of the model, and <math>\varepsilon_t</math> is white noise.<ref>{{Cite book |last=Box |first=George E. P. |title=Time series analysis : forecasting and control |date=1994 |publisher=Prentice Hall |others=Gwilym M. Jenkins, Gregory C. Reinsel |isbn=0-13-060774-6 |edition=3rd |location=Englewood Cliffs, N.J. |pages=54 |language=en |oclc=28888762}}</ref><ref>{{Cite book |author1-link=Robert H. Shumway |last=Shumway |first=Robert H. |title=Time series analysis and its applications |date=2000 |publisher=Springer |others=David S. Stoffer |isbn=0-387-98950-1 |location=New York |pages=90–91 |language=en |oclc=42392178 }}</ref> This can be equivalently written using the backshift operator ''B'' as
:<math> X_t = \sum_{i=1}^p \varphi_i B^i X_t + \varepsilon_t </math>
so that, moving the summation term to the left side and using polynomial notation, we have
:<math>\varphi(B)X_t= \varepsilon_t</math>
An autoregressive model can thus be viewed as the output of an all-pole infinite impulse response filter whose input is white noise.
Some parameter constraints are necessary for the model to remain weak-sense stationary. For example, processes in the AR(1) model with <math>|\varphi_1 | \geq 1</math> are not stationary. More generally, for an AR(''p'') model to be weak-sense stationary, the roots of the polynomial <math>\Phi(z):=\textstyle 1 - \sum_{i=1}^p \varphi_i z^{i}</math> must lie outside the unit circle, i.e., each (complex) root <math>z_i</math> must satisfy <math>|z_i |>1</math> (see pages 89,92 <ref>{{cite book |last1=Shumway |first1=Robert H. |last2=Stoffer |first2=David |title=Time series analysis and its applications : with R examples |date=2010 |publisher=Springer |isbn=978-1441978646 |edition=3rd}}</ref>).
==Intertemporal effect of shocks==
In an AR process, a one-time shock affects values of the evolving variable infinitely far into the future. For example, consider the AR(1) model <math> X_t = \varphi_1 X_{t-1} + \varepsilon_t</math>. A non-zero value for <math>\varepsilon_t</math> at say time ''t''=1 affects <math>X_1</math> by the amount <math>\varepsilon_1</math>. Then by the AR equation for <math>X_2</math> in terms of <math>X_1</math>, this affects <math>X_2</math> by the amount <math>\varphi_1 \varepsilon_1</math>. Then by the AR equation for <math>X_3</math> in terms of <math>X_2</math>, this affects <math>X_3</math> by the amount <math>\varphi_1^2 \varepsilon_1</math>. Continuing this process shows that the effect of <math>\varepsilon_1</math> never ends, although if the process is stationary then the effect diminishes toward zero in the limit.
Because each shock affects ''X'' values infinitely far into the future from when they occur, any given value ''X''<sub>''t''</sub> is affected by shocks occurring infinitely far into the past. This can also be seen by rewriting the autoregression
:<math>\varphi (B)X_t= \varepsilon_t \,</math>
(where the constant term has been suppressed by assuming that the variable has been measured as deviations from its mean) as
:<math>X_t= \frac{1}{\varphi (B)}\varepsilon_t \, .</math>
When the polynomial division on the right side is carried out, the polynomial in the backshift operator applied to <math>\varepsilon_t</math> has an infinite order—that is, an infinite number of lagged values of <math>\varepsilon_t</math> appear on the right side of the equation.
==Characteristic polynomial== The autocorrelation function of an AR(''p'') process can be expressed as {{Citation needed|date=October 2011|reason=a_k not defined and seems wrong}} :<math>\rho(\tau) = \sum_{k=1}^p a_k y_k^{-|\tau|} ,</math>
where <math>y_k</math> are the roots of the polynomial : <math>\varphi(B) = 1- \sum_{k=1}^p \varphi_k B^k </math>
where ''B'' is the backshift operator, where <math>\varphi(\cdot)</math> is the function defining the autoregression, and where <math>\varphi_k</math> are the coefficients in the autoregression. The formula is valid only if all the roots have multiplicity 1.{{Citation needed|date=July 2022|reason=it is true but needs reference, see talk page}}
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==
thumb|right|alt="Figure has 5 plots of AR processes. AR(0) and AR(0.3) are white noise or look like white noise. AR(0.9) has some large scale oscillating structure."|AR(0); AR(1) with AR parameter 0.3; AR(1) with AR parameter 0.9; AR(2) with AR parameters 0.3 and 0.3; and AR(2) with AR parameters 0.9 and −0.8 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 <math>\varphi</math>, only the previous term in the process and the noise term contribute to the output. If <math>\varphi</math> is close to 0, then the process still looks like white noise, but as <math>\varphi</math> 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 <math>\varphi_1</math> and <math>\varphi_2</math> are positive, the output will resemble a low pass filter, with the high frequency part of the noise decreased. If <math>\varphi_1</math> is positive while <math>\varphi_2</math> is negative, then the process favors changes in sign between terms of the process. The output oscillates. This can be linked to edge detection or detection of change in direction.
==Example: An AR(1) process==
An AR(1) process is given by:<math display="block">X_t = \varphi X_{t-1}+\varepsilon_t\,</math>where <math>\varepsilon_t</math> is a white noise process with zero mean and constant variance <math>\sigma_\varepsilon^2</math>. (Note: The subscript on <math>\varphi_1</math> has been dropped.) The process is weak-sense stationary if <math>|\varphi|<1</math> since it is obtained as the output of a stable filter whose input is white noise. (If <math>\varphi=1</math> then the variance of <math>X_t</math> depends on time lag ''t'', so that the variance of the series diverges to infinity as ''t'' goes to infinity, and is therefore not weak-sense stationary.) Assuming <math>|\varphi|<1</math>, the mean <math>\operatorname{E} (X_t)</math> is identical for all values of ''t'' by definition of weak sense stationarity. If the mean is denoted by <math>\mu</math>, it follows from<math display="block">\operatorname{E} (X_t)=\varphi\operatorname{E} (X_{t-1})+\operatorname{E}(\varepsilon_t), </math>that<math display="block"> \mu=\varphi\mu+0,</math>and hence
:<math>\mu=0.</math>
The variance is
:<math>\textrm{var}(X_t)=\operatorname{E}(X_t^2)-\mu^2=\frac{\sigma_\varepsilon^2}{1-\varphi^2},</math> where <math>\sigma_\varepsilon</math> is the standard deviation of <math>\varepsilon_t</math>. This can be shown by noting that :<math>\textrm{var}(X_t) = \varphi^2\textrm{var}(X_{t-1}) + \sigma_\varepsilon^2,</math> and then by noticing that the quantity above is a stable fixed point of this relation.
The autocovariance is given by
:<math>B_n=\operatorname{E}(X_{t+n}X_t)-\mu^2=\frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\,\varphi^{|n|}.</math>
It can be seen that the autocovariance function decays with a decay time (also called time constant) of <math>\tau=1/(1-\varphi)</math>.<ref>Lai, Dihui; and Lu, Bingfeng; [https://www.soa.org/globalassets/assets/library/newsletters/predictive-analytics-and-futurism/2017/june/2017-predictive-analytics-iss15-lai-lu.pdf "Understanding Autoregressive Model for Time Series as a Deterministic Dynamic System"] {{Webarchive|url=https://web.archive.org/web/20230324041726/https://www.soa.org/globalassets/assets/library/newsletters/predictive-analytics-and-futurism/2017/june/2017-predictive-analytics-iss15-lai-lu.pdf |date=2023-03-24 }}, in ''Predictive Analytics and Futurism'', June 2017, number 15, June 2017, pages 7-9</ref>
The spectral density function is the Fourier transform of the autocovariance function. In discrete terms this will be the discrete-time Fourier transform:
:<math>\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). </math>
This expression is periodic due to the discrete nature of the <math>X_j</math>, which is manifested as the cosine term in the denominator. If we assume that the sampling time (<math>\Delta t=1</math>) is much smaller than the decay time (<math>\tau</math>), then we can use a continuum approximation to <math>B_n</math>:
:<math>B(t)\approx \frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\,\varphi^{|t|}</math>
which yields a Lorentzian profile for the spectral density:
:<math>\Phi(\omega)= \frac{1}{\sqrt{2\pi}}\,\frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\frac{\gamma}{\pi(\gamma^2+\omega^2)}</math>
where <math>\gamma=1/\tau</math> is the angular frequency associated with the decay time <math>\tau</math>.
An alternative expression for <math>X_t</math> can be derived by first substituting <math>\varphi X_{t-2}+\varepsilon_{t-1}</math> for <math>X_{t-1}</math> in the defining equation. Continuing this process ''N'' times yields
:<math>X_t=\varphi^NX_{t-N}+\sum_{k=0}^{N-1}\varphi^k\varepsilon_{t-k}.</math>
For ''N'' approaching infinity, <math>\varphi^N</math> will approach zero and:
:<math>X_t=\sum_{k=0}^\infty\varphi^k\varepsilon_{t-k}.</math>
It is seen that <math>X_t</math> is white noise convolved with the <math>\varphi^k</math> kernel plus the constant mean. If the white noise <math>\varepsilon_t</math> is a Gaussian process then <math>X_t</math> is also a Gaussian process. In other cases, the central limit theorem indicates that <math>X_t</math> will be approximately normally distributed when <math>\varphi</math> is close to one.
For <math>\varepsilon_t = 0</math>, the process <math>X_t = \varphi X_{t-1}</math> will be a geometric progression (''exponential'' growth or decay). In this case, the solution can be found analytically: <math>X_t = a \varphi^t</math> whereby <math>a</math> is an unknown constant (initial condition).
=== Explicit mean/difference form of AR(1) process ===
The AR(1) model is the discrete-time analogy of the continuous Ornstein-Uhlenbeck process. It is therefore sometimes useful to understand the properties of the AR(1) model cast in an equivalent form. In this form, the AR(1) model, with process parameter <math>\theta \in \mathbb{R}</math>, is given by
:<math>X_{t+1} = X_t + (1-\theta)(\mu - X_t) + \varepsilon_{t+1}</math>, where <math>|\theta| < 1 \,</math>, <math>\mu := E(X)</math> is the model mean, and <math>\{\varepsilon_t \}</math> is a white-noise process with zero mean and constant variance <math>\sigma</math>.
By rewriting this as <math> X_{t+1} = \theta X_t + (1 - \theta)\mu + \varepsilon_{t+1} </math> and then deriving (by induction) <math>X_{t+n} = \theta^n X_t + (1 - \theta^n) \mu + \sum_{i = 1}^n \left(\theta^{n - i} \varepsilon_{t + i} \right)</math>, one can show that
:<math> \operatorname{E}(X_{t+n} \mid X_t) = \mu\left[1-\theta^n\right] + X_t\theta^n</math> and :<math> \operatorname{Var} (X_{t+n} \mid X_t) = \sigma^2 \frac{ 1 - \theta^{2n} }{1 -\theta^2}. </math>
==Choosing the maximum lag== {{Main|Partial autocorrelation function}}
The partial autocorrelation of an AR(p) process equals zero at lags larger than ''p'', so the appropriate maximum lag ''p'' is the one after which the partial autocorrelations are all zero.
==Calculation of the AR parameters==
There are many ways to estimate the coefficients, such as the ordinary least squares procedure or method of moments (through Yule–Walker equations).
The AR(''p'') model is given by the equation
:<math> X_t = \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t.\,</math>
It is based on parameters <math>\varphi_i</math> 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=== <!-- this heading is linked from other articles -->
The Yule–Walker equations, named for Udny Yule and Gilbert Walker,<ref>Yule, G. Udny (1927) [http://visualiseur.bnf.fr/Visualiseur?Destination=Gallica&O=NUMM-56031 "On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer's Sunspot Numbers"] {{Webarchive|url=https://web.archive.org/web/20110514094546/http://visualiseur.bnf.fr/Visualiseur?Destination=Gallica&O=NUMM-56031 |date=2011-05-14 }}, ''Philosophical Transactions of the Royal Society of London'', Ser. A, Vol. 226, 267–298.]</ref><ref>Walker, Gilbert (1931) [http://visualiseur.bnf.fr/Visualiseur?Destination=Gallica&O=NUMM-56224 "On Periodicity in Series of Related Terms"] {{Webarchive|url=https://web.archive.org/web/20110607170511/http://visualiseur.bnf.fr/Visualiseur?Destination=Gallica&O=NUMM-56224 |date=2011-06-07 }}, ''Proceedings of the Royal Society of London'', Ser. A, Vol. 131, 518–532.</ref> are the following set of equations.<ref>{{cite book |last=Theodoridis |first=Sergios |title=Machine Learning: A Bayesian and Optimization Perspective |publisher=Academic Press, 2015 |chapter=Chapter 1. Probability and Stochastic Processes |pages=9–51 |isbn=978-0-12-801522-3 |date=2015-04-10 }}</ref>
:<math>\gamma_m = \sum_{k=1}^p \varphi_k \gamma_{m-k} + \sigma_\varepsilon^2\delta_{m,0},</math>
where {{nowrap|1=''m'' = 0, …, ''p''}}, yielding {{nowrap|''p'' + 1}} equations. Here <math>\gamma_m</math> is the autocovariance function of X<sub>t</sub>, <math>\sigma_\varepsilon</math> is the standard deviation of the input noise process, and <math>\delta_{m,0}</math> is the Kronecker delta function.
Because the last part of an individual equation is non-zero only if {{nowrap|1=''m'' = 0}}, the set of equations can be solved by representing the equations for {{nowrap|''m'' > 0}} in matrix form, thus getting the equation
:<math>\begin{bmatrix} \gamma_1 \\ \gamma_2 \\ \gamma_3 \\ \vdots \\ \gamma_p \\ \end{bmatrix}
=
\begin{bmatrix} \gamma_0 & \gamma_{-1} & \gamma_{-2} & \cdots \\ \gamma_1 & \gamma_0 & \gamma_{-1} & \cdots \\ \gamma_2 & \gamma_1 & \gamma_0 & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ \gamma_{p-1} & \gamma_{p-2} & \gamma_{p-3} & \cdots \\ \end{bmatrix}
\begin{bmatrix} \varphi_1 \\ \varphi_2 \\ \varphi_3 \\ \vdots \\ \varphi_p \\ \end{bmatrix}
</math>
which can be solved for all <math>\{\varphi_m; m=1,2, \dots ,p\}.</math> The remaining equation for ''m'' = 0 is
:<math>\gamma_0 = \sum_{k=1}^p \varphi_k \gamma_{-k} + \sigma_\varepsilon^2 ,</math>
which, once <math>\{\varphi_m ; m=1,2, \dots ,p \}</math> are known, can be solved for <math>\sigma_\varepsilon^2 .</math>
An alternative formulation is in terms of the autocorrelation function. The AR parameters are determined by the first ''p''+1 elements <math>\rho(\tau)</math> of the autocorrelation function. The full autocorrelation function can then be derived by recursively calculating <ref name=Storch>{{Cite book | publisher = Cambridge University Press | isbn = 0-521-01230-9 | last1 = Von Storch | first1 = Hans | first2=Francis W. |last2=Zwiers | title = Statistical analysis in climate research | year = 2001 | doi = 10.1017/CBO9780511612336 }}{{Page needed|date=March 2011}}</ref> : <math>\rho(\tau) = \sum_{k=1}^p \varphi_k \rho(k-\tau)</math> Examples for some Low-order AR(''p'') processes * ''p''=1 ** <math>\gamma_1 = \varphi_1 \gamma_0</math> ** Hence <math>\rho_1 = \gamma_1 / \gamma_0 = \varphi_1</math> * ''p''=2 ** The Yule–Walker equations for an AR(2) process are **: <math>\gamma_1 = \varphi_1 \gamma_0 + \varphi_2 \gamma_{-1}</math> **: <math>\gamma_2 = \varphi_1 \gamma_1 + \varphi_2 \gamma_0</math> *** Remember that <math>\gamma_{-k} = \gamma_k</math> *** Using the first equation yields <math>\rho_1 = \gamma_1 / \gamma_0 = \frac{\varphi_1}{1-\varphi_2}</math> *** Using the recursion formula yields <math>\rho_2 = \gamma_2 / \gamma_0 = \frac{\varphi_1^2 - \varphi_2^2 + \varphi_2}{1 - \varphi_2}</math>
===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.<ref>{{Cite web |url=http://www-stat.wharton.upenn.edu/~steele/Courses/956/Resource/YWSourceFiles/YW-Eshel.pdf |title=The Yule Walker Equations for the AR Coefficients |last=Eshel |first=Gidon |website=stat.wharton.upenn.edu |access-date=2019-01-27 |archive-date=2018-07-13 |archive-url=https://web.archive.org/web/20180713135223/http://www-stat.wharton.upenn.edu/~steele/Courses/956/Resource/YWSourceFiles/YW-Eshel.pdf |url-status=live }}</ref> 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''<sub>''t''</sub> 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: ::<math> X_t = \sum_{i=1}^p \varphi_i X_{t+i}+ \varepsilon^*_t \,.</math> :Here predicted values of ''X''<sub>''t''</sub> would be based on the ''p'' future values of the same series.{{Clarify|reason=1) the latest equation was introduced as "the backward representation of the AR model"; 2) assuming t is counted in the same direction as the progression of time, then the subscripts t-i indicate past values of X w.r.t. X_t; 3) the last sentence would indicate that - suddenly - X_t would depend on future (rather than past) values of X|date=July 2020}} This way of estimating the AR parameters is due to John Parker Burg,<ref name=Burg/> and is called the Burg method:<ref name=Brockwell/> Burg and later authors called these particular estimates "maximum entropy estimates",<ref name=Burg1/> 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.<ref name=Bos/>
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 corresponding 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 corresponding 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== thumb|right thumb|right
The power spectral density (PSD) of an AR(''p'') process with noise variance <math>\mathrm{Var}(Z_t) = \sigma_Z^2</math> is<ref name=Storch/> : <math>S(f) = \frac{\sigma_Z^2}{| 1-\sum_{k=1}^p \varphi_k e^{-i 2 \pi f k} |^2}.</math>
===AR(0)=== For white noise (AR(0)) : <math>S(f) = \sigma_Z^2.</math>
===AR(1)=== For AR(1) : <math>S(f) = \frac{\sigma_Z^2}{| 1- \varphi_1 e^{-2 \pi i f} |^2} = \frac{\sigma_Z^2}{ 1 + \varphi_1^2 - 2 \varphi_1 \cos 2 \pi f }</math> *If <math>\varphi_1 > 0</math> there is a single spectral peak at <math>f=0</math>, often referred to as red noise. As <math>\varphi_1 </math> becomes nearer 1, there is stronger power at low frequencies, i.e. larger time lags. This is then a low-pass filter, when applied to full spectrum light, everything except for the red light will be filtered. *If <math>\varphi_1 < 0</math> there is a minimum at <math>f=0</math>, often referred to as blue noise. This similarly acts as a high-pass filter, everything except for blue light will be filtered.
===AR(2)=== The behavior of an AR(2) process is determined entirely by the roots of it characteristic equation, which is expressed in terms of the lag operator as:
:<math> 1 - \varphi_1 B -\varphi_2 B^2 =0, </math>
or equivalently by the poles of its transfer function, which is defined in the Z domain by:
:<math> H_z = (1 - \varphi_1 z^{-1} -\varphi_2 z^{-2})^{-1}. </math>
It follows that the poles are values of z satisfying:
:<math> 1 - \varphi_1 z^{-1} -\varphi_2 z^{-2} = 0, </math>
which yields:
:<math> z_1,z_2 = \frac{1}{2\varphi_2}\left(\varphi_1 \pm \sqrt{\varphi_1^2 + 4\varphi_2} \, \right). </math>
<math>z_1</math> and <math>z_2</math> are the reciprocals of the characteristic roots, as well as the eigenvalues of the temporal update matrix:
:<math> \begin{bmatrix} \varphi_1 & \varphi_2 \\ 1 & 0 \end{bmatrix} </math>
AR(2) processes can be split into three groups depending on the characteristics of their roots/poles:
* When <math>\varphi_1^2 + 4\varphi_2 < 0</math>, the process has a pair of complex-conjugate poles, creating a mid-frequency peak at: :<math>f^* = \frac{1}{2\pi}\cos^{-1}\left(\frac{\varphi_1}{2\sqrt{-\varphi_2}}\right),</math>
with bandwidth about the peak inversely proportional to the moduli of the poles:
:<math>|z_1|=|z_2|=\sqrt{-\varphi_2}.</math>
The terms involving square roots are all real in the case of complex poles since they exist only when <math>\varphi_2<0</math>.
Otherwise the process has real roots, and: * When <math>\varphi_1 > 0</math> it acts as a low-pass filter on the white noise with a spectral peak at <math>f=0</math> * When <math>\varphi_1 < 0</math> it acts as a high-pass filter on the white noise with a spectral peak at <math>f=1/2</math>. The process is non-stationary when the poles are on or outside the unit circle, or equivalently when the characteristic roots are on or inside the unit circle. The process is stable when the poles are strictly within the unit circle (roots strictly outside the unit circle), or equivalently when the coefficients are in the triangle <math>-1 \le \varphi_2 \le 1 - |\varphi_1|</math>.
The full PSD function can be expressed in real form as: :<math>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)}</math>
==Implementations in statistics packages== * R – the ''stats'' package includes ''ar'' function;<ref>[http://finzi.psych.upenn.edu/R/library/stats/html/ar.html "Fit Autoregressive Models to Time Series"] {{Webarchive|url=https://web.archive.org/web/20160128234632/http://finzi.psych.upenn.edu/R/library/stats/html/ar.html |date=2016-01-28}} (in R)</ref> the ''astsa'' package includes ''sarima'' function to fit various models including AR.<ref name=":0">{{Cite web |last1=Stoffer |first1=David |last2=Poison |first2=Nicky |date=2023-01-09 |title=astsa: Applied Statistical Time Series Analysis |url=https://cran.r-project.org/web/packages/astsa/ |access-date=2023-08-20}}</ref> * MATLAB – the Econometrics Toolbox<ref>{{Cite web |url=https://www.mathworks.com/products/econometrics.html |title=Econometrics Toolbox |website=www.mathworks.com |access-date=2022-02-16 |archive-date=2023-04-16 |archive-url=https://web.archive.org/web/20230416160907/https://www.mathworks.com/products/econometrics.html |url-status=live}}</ref> and System Identification Toolbox<ref>{{Cite web |url=https://www.mathworks.com/products/sysid.html |title=System Identification Toolbox |website=www.mathworks.com |access-date=2022-02-16 |archive-date=2022-02-16 |archive-url=https://web.archive.org/web/20220216063519/https://www.mathworks.com/products/sysid.html |url-status=live}}</ref> include AR models.<ref>{{Cite web |url=https://www.mathworks.com/help/econ/autoregressive-model.html |title=Autoregressive Model - MATLAB & Simulink |website=www.mathworks.com |access-date=2022-02-16 |archive-date=2022-02-16 |archive-url=https://web.archive.org/web/20220216063648/https://www.mathworks.com/help/econ/autoregressive-model.html |url-status=live}}</ref> * MATLAB and Octave – the ''TSA'' toolbox contains several estimation functions for uni-variate, multivariate, and adaptive AR models.<ref>{{Cite web |url=http://pub.ist.ac.at/~schloegl/matlab/tsa/ |title=The Time Series Analysis (TSA) toolbox for Octave and MATLAB |website=pub.ist.ac.at |access-date=2012-04-03 |archive-date=2012-05-11 |archive-url=https://web.archive.org/web/20120511144225/http://pub.ist.ac.at/~schloegl/matlab/tsa/ |url-status=live}}</ref> * PyMC3 – the Bayesian statistics and probabilistic programming framework supports AR modes with ''p'' lags. * ''bayesloop'' – supports parameter inference and model selection for the AR-1 process with time-varying parameters.<ref>{{Cite web |url=https://github.com/christophmark/bayesloop |title=christophmark/bayesloop |date=December 7, 2021 |via=GitHub |access-date=September 4, 2018 |archive-date=September 28, 2020 |archive-url=https://web.archive.org/web/20200928085417/https://github.com/christophmark/bayesloop |url-status=live}}</ref> *Python – statsmodels.org hosts an AR model.<ref>{{Cite web |url=https://www.statsmodels.org/stable/generated/statsmodels.tsa.ar_model.AutoReg.html |title=statsmodels.tsa.ar_model.AutoReg — statsmodels 0.12.2 documentation |website=www.statsmodels.org |access-date=2021-04-29 |archive-date=2021-02-28 |archive-url=https://web.archive.org/web/20210228123354/https://www.statsmodels.org/stable/generated/statsmodels.tsa.ar_model.AutoReg.html |url-status=live}}</ref>
==Impulse response==
The impulse response of a system is the change in an evolving variable in response to a change in the value of a shock term ''k'' periods earlier, as a function of ''k''. Since the AR model is a special case of the vector autoregressive model, the computation of the impulse response in vector autoregression#impulse response applies here.
==''n''-step-ahead forecasting==
Once the parameters of the autoregression
:<math> X_t = \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t \,</math>
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 preceding values ''X''<sub>''t-i''</sub> for ''i=''1, ..., ''p'' into the autoregressive equation while setting the error term <math>\varepsilon_t</math> equal to zero (because we forecast ''X''<sub>''t''</sub> 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 preceding steps.
There are four sources of uncertainty regarding predictions obtained in this manner: (1) uncertainty as to whether the autoregressive model is the correct model; (2) uncertainty about the accuracy of the forecasted values that are used as lagged values in the right side of the autoregressive equation; (3) uncertainty about the true values of the autoregressive coefficients; and (4) uncertainty about the value of the error term <math>\varepsilon_t \,</math> for the period being predicted. Each of the last three can be quantified and combined to give a confidence interval for the ''n''-step-ahead predictions; the confidence interval will become wider as ''n'' increases because of the use of an increasing number of estimated values for the right-side variables.
==See also== * Moving average model * Linear difference equation * Predictive analytics * Linear predictive coding * Resonance * Levinson recursion * Ornstein–Uhlenbeck process * Infinite impulse response * PagedAttention / vAttention — attention algorithms for efficient serving of large language models
==Notes== {{Reflist | refs=
<ref name=Brockwell>{{cite journal |first1=Peter J. |last1=Brockwell |first2=Rainer |last2=Dahlhaus |first3=A. Alexandre |last3=Trindade |journal=Statistica Sinica |volume=15 |year=2005 |pages=197–213 |title=Modified Burg Algorithms for Multivariate Subset Autoregression |url=http://www3.stat.sinica.edu.tw/statistica/oldpdf/A15n112.pdf |url-status=dead |archive-url=https://web.archive.org/web/20121021015413/http://www3.stat.sinica.edu.tw/statistica/oldpdf/A15n112.pdf |archive-date=2012-10-21 }}</ref>
<ref name=Bos>{{Cite journal |last1=Bos |first1=Robert |last2=De Waele |first2=Stijn |last3=Broersen |first3=Piet M. T. |doi=10.1109/TIM.2002.808031 |title=Autoregressive spectral estimation by application of the Burg algorithm to irregularly sampled data |journal=IEEE Transactions on Instrumentation and Measurement |volume=51 |issue=6 |pages=1289 |year=2002 |bibcode=2002ITIM...51.1289B |url=http://resolver.tudelft.nl/uuid:870559f7-f1e9-4968-83da-edd19485eaaf |access-date=2019-12-11 |archive-date=2023-04-16 |archive-url=https://web.archive.org/web/20230416160911/https://repository.tudelft.nl/islandora/object/uuid:870559f7-f1e9-4968-83da-edd19485eaaf?collection=research |url-status=live }}</ref>
<ref name=Burg>Burg, John Parker (1968); "A new analysis technique for time series data", in ''Modern Spectrum Analysis'' (Edited by D. G. Childers), NATO Advanced Study Institute of Signal Processing with emphasis on Underwater Acoustics. IEEE Press, New York.</ref>
<ref name=Burg1>Burg, John Parker (1967) "Maximum Entropy Spectral Analysis", ''Proceedings of the 37th Meeting of the Society of Exploration Geophysicists'', Oklahoma City, Oklahoma.</ref> }}
==References== * {{cite book |last=Mills |first=Terence C. |year=1990 |title=Time Series Techniques for Economists |url=https://archive.org/details/timeseriestechni0000mill |url-access=registration |publisher=Cambridge University Press |isbn=9780521343398 }} * {{cite book |last1=Percival |first1=Donald B. |first2=Andrew T. |last2=Walden |year=1993 |title=Spectral Analysis for Physical Applications |publisher=Cambridge University Press |bibcode=1993sapa.book.....P }} * {{cite book |last1=Pandit |first1=Sudhakar M. |last2=Wu |first2=Shien-Ming |year=1983 |title=Time Series and System Analysis with Applications |publisher=John Wiley & Sons }}
==External links== * [http://paulbourke.net/miscellaneous/ar/ AutoRegression Analysis (AR)] by Paul Bourke * {{YouTube|id=b8yslhlIsA0&list=PLD15D38DC7AA3B737&index=8#t=34m25s|title=Econometrics lecture (topic: Autoregressive models)}} by Mark Thoma
{{Stochastic processes|state=collapsed}} {{Artificial intelligence navbox}}
Category:Autocorrelation Category:Signal processing