Machine Learning & Signals Learning

21 ARMA Model

This chapter extends least-squares (LS)-based modeling to systems and time-series contexts, focusing on linear prediction coefficients derived from signals. Instead of restricting to a predefined parametric form like a sinusoid, here we consider arbitrary zero-mean finite-time signals \(x[n]\) and \(y[n], n=0,\ldots , L-1\), in the presence of zero-mean noise, \(\epsilon [n]\).

21.1 Auto-Correlation Function

21.1.1 AR(1)

In system modeling and time-series analysis, a common approach is to express the current sample of a signal in terms of its past values. A simple first-order autoregressive (AR(1)) model predicts the current sample \(x[n]\) from a single past sample \(x[n-1]\), using the system model

\begin{equation} \label {eq:ar1} x[n]=h_1x[n-1] + \epsilon [n]. \end{equation}

and

\begin{equation} \hat {x}[n] = h_1x[n-1] \end{equation}

Here, \(h_1\) is the linear prediction coefficient that indicates how strongly the previous sample \(x[n-1]\) influences the current sample \(x[n]\).

To determine \(h_1\), minimization of SE loss function

\begin{equation} \loss (h_1) = \sum _n \left (x[n] - h_1x[n-1])\right )^2 \end{equation}

is used. Setting the derivative of \(\loss (h_1)\) with respect to \(h_1\) to zero leads to

\begin{equation} \label {eq-a1} \frac {d\mathcal {L}(a)}{da} =2\sum _n (x[n]-h_1x[n-1])(-x[n-1])=0 \end{equation}

Finally, the LS solution for \(h_1\) is

\begin{equation} \label {eq-acf-a1} h_1 = \frac {\sum _n x[n]x[n-1]}{\sum _n x^2[n-1]}. \end{equation}

The value \(h_1\) is termed as the (single) prediction coefficient of AR(1) model.

Matrix formulation

The coefficient may be formulated as a relation between two vectors

\begin{equation} \label {eq:ar1:matrix} \hat {\by }= \begin{bmatrix} x[1] \\ x[2] \\ \vdots \\ x[L-2]\\ x[L-1] \end {bmatrix} \quad \bx = \begin{bmatrix} x[0]\\ x[1] \\ \vdots \\ x[L-3] \\ x[L-2] \end {bmatrix}, \end{equation}

such that

\begin{equation} \hat {\by } = h_1 \bx . \end{equation}

The correspoding MMSE loss is

\begin{equation} \loss (h_1) = \norm {\by - h_1 \bx }^2 \end{equation}

and its minimum is given by normal equation, (Eq. (3.5)),

\begin{equation} \label {eq:normal_eq_vector} h_1 = {\left (\bx ^T\bx \right )^{-1}}\bx ^T\by . \end{equation}

The substitution results in the solution identical to Eq. (21.5).

Prediction with \(k\)-samples difference

To generalize beyond a single lag, we use one-coefficient prediction from the current sample \(x[n]\) and its time-shifted version \(x[n-k]\), using the \(k\)-step predictor,

\begin{equation} \hat {x}[n]=h_kx[n-k]. \end{equation}

The corresponding solution is very similar to the solution for \(k=1\),

\begin{equation} \label {eq-acf-ak} h_k = \frac {\sum _n x[n]x[n-k]}{\sum _n x^2[n-k]}. \end{equation}

The linear dependence between \(x[n]\) and \(x[n-k]\), and the prediction it enables, are illustrated in Fig. 21.1 for an AR(1) signal. The dependence weakens as the lag grows, and with it the accuracy of the prediction.

21.1.2 Auto-correlation function (ACF)

The nominator of Eq. (21.11) is termed (raw, or unscaled) auto-correlation function (ACF) at lag \(k\),

\begin{equation} R_{\bx \bx }[k] = \sum _n x[n]x[n-k],\quad k=0,\ldots ,L-1 \end{equation}

The ACF provides a measure of how linearly dependent the signal is on its shifted versions. Variants like biased, unbiased, and normalized ACF offer different normalization schemes to handle finite data length and scaling issues.

The signal energy is given by

\begin{equation} \label {eq:arma:energy} E_\bx = \sum _{n=0}^{L-1} x^2[n] = R_\bxx [0] \end{equation}

Biased auto-correlation Another useful definition is averaged sum of \(y[n]y[n-k]\),

\begin{equation} \label {eq-acf-biased} \begin{aligned} R_{\bx \bx ,biased}[k] &= \frac {1}{L}\sum _n x[n]x[n-k],\quad k=0,\ldots ,L-1\\ &= \frac {1}{L}R_{\bx \bx }[k] \end {aligned} \end{equation}

is termed biased (see Sec. 1.1.3 for bias definition) ACF. The corresponding average power is given by

\begin{equation} P_\bx = \frac {1}{L}\sum _n x^2[n] = R_{\bxx ,biased}[0] \end{equation}

Normalized auto-correlation Another version is normalized ACF of the form

\begin{equation} \label {eq-acf-norm} R_{\bx \bx ,norm}[k] = \frac {R_{\bx \bx }[k]}{R_{\bx \bx }[0]}\lesssim h_k. \end{equation}

Note, the difference between denominator above and the one in Eq. (21.11) (where the lagged sum runs over \(n=k,\ldots ,L-1\)) is

\begin{equation} \sum _{n=0}^{L-1} x^2[n] - \sum _{n=k}^{L-1}x^2[n-k] = x^2[L-k] + \cdots + x^2[L-1], \end{equation}

i.e. the last \(k\) samples of the squared signal. Under stationarity and \(k\ll L\), these boundary terms are negligible relative to \(R_{\bxx }[0]\), so \(h_k\) in Eq. (21.11) and \(R_{\bxx ,norm}[k]\) are approximately equal for a sufficiently high \(L\). Moreover, here the denominator is \(k\)-independent.

Unbiased auto-correlation Note, that the ACF includes summation only of \(n-k\) terms. Another useful normalization is

\begin{equation} \label {eq-acf-unbiased} \begin{aligned} R_{\bx \bx ,unbiased}[k] &= \frac {1}{L-k}\sum _n x[n]x[n-k]\\ &= \frac {1}{L-k}R_{\bx \bx }[k] \end {aligned} \end{equation}

This time it is assumed that

\begin{equation} \label {eq:acf-stationarity} \begin{aligned} \frac {1}{L}\sum _n x^2[n] &\approx \frac {1}{L-k}\sum _n x^2[n-k]\\ \frac {x^2[0] +x^2[1] + \cdots x^2[L-1]}{L} &\approx \frac {x^2[k] + \cdots x^2[L-1]}{L-k} \end {aligned} \end{equation}

and the resulting expression

\begin{equation} \begin{aligned} \frac {L}{L-k}\frac {R_{\bx \bx }[k]}{R_{\bx \bx }[0]}=\frac {L}{L-k}R_{\bx \bx ,norm}[k]\\ =\frac {R_{\bx \bx ,unbiased}[k]}{R_{\bx \bx ,unbiased}[0]} \approx h_k \end {aligned} \end{equation}

is assumed to be closer approximation to \(h_k\) than the normalized auto-correlation, for a relatively small values of \(k\).

An AR(1) signal and its normalized ACF are shown in Fig. 21.2. The ACF starts at \(R_{\bxx ,norm}[0]=1\) and decays with the lag, reproducing the loss of predictability seen in Fig. 21.1.

Correlation time The correlation time, \(k_c\) is the lag where \(R_{\bxx ,norm}[k_c]\) falls below a threshold,

\begin{equation} \label {eq-correlation-time} R_{\bxx ,norm}[k_c] = 0.5 \text { or } 0.1 \text { or }\exp (-1) \end{equation}

The predictability is assumed negligible for \(k>k_c\),

\begin{equation} \label {eq-correlation-time-limit} R_{\bxx ,norm}[k>k_c]\approx 0 \end{equation}

The decision threshold depends on the field of application. For example, one of the most applicable ACF for natural processes is

\begin{equation} R_{\bxx ,norm}[k] = \exp (-\frac {k}{k_c}) \end{equation}

with \(R_{\bxx ,norm}[k_c] = \exp (-1)\).

Note, correlation time is mostly used in physics and engineering models and is less applicable in social sciences.

Correlation Coefficient Interpretation

All normalizations above still depend on the amplitude of the signal: replacing \(x[n]\) by \(ax[n]\) multiplies \(R_{\bxx }[k]\) by \(a^2\). Comparing the predictability of two different signals, or of the same signal recorded at a different gain, therefore requires a quantity on a fixed scale. The correlation coefficient achieves this by dividing the raw ACF by the energies of both vectors being compared,

\begin{equation} \label {eq:corr-coeff-def} \rho _{\bxx }[k] \triangleq \frac {\sum _n x[n]x[n-k]}{\sqrt {\sum _n x^2[n]\,\sum _n x^2[n-k]}} = \frac {R_{\bxx }[k]}{\sqrt {R_{\bxx }[0]\cdot \sum _n x^2[n-k]}}. \end{equation}

The denominator combines both energies involved: \(R_{\bxx }[0]=\sum _n x^2[n]\), the energy of the signal (Eq. (21.13)), and \(\sum _n x^2[n-k]\), the energy of its lagged copy. Taking the square root of their product, instead of dividing by either energy alone, treats the two vectors symmetrically, so that exchanging them leaves \(\rho _{\bxx }[k]\) unchanged. The result is always bounded,

\begin{equation} \label {eq:rho-bounded} -1\le \rho _{\bxx }[k]\le 1, \end{equation}

where \(\abs {\rho _{\bxx }[k]}=1\) is reached only if \(x[n]\) and \(x[n-k]\) are exactly proportional.

Geometric mean and the Cauchy–Schwarz bound (*) The geometric mean of two non-negative numbers \(a\) and \(b\) is \(\sqrt {ab}\), the multiplicative counterpart of the arithmetic mean \((a+b)/2\). The denominator of Eq. (21.24) is therefore the geometric mean of the two energies, and the bound Eq. (21.25) is the Cauchy–Schwarz inequality applied to \(x[n]\) and \(x[n-k]\).

To prove it, note that for any real \(\lambda \) a sum of squares cannot be negative,

\begin{equation} \label {eq:cs-quadratic} q(\lambda ) = \sum _n \left (\lambda \,x[n] + x[n-k]\right )^2 = \lambda ^2 R_{\bxx }[0] + 2\lambda \,R_{\bxx }[k] + \sum _n x^2[n-k] \;\ge \; 0 . \end{equation}

Assume \(R_{\bxx }[0]>0\), so that \(q(\lambda )\) is a parabola opening upwards; otherwise \(x[n]\equiv 0\) and \(\rho _{\bxx }[k]\) is undefined. A parabola that is never negative crosses zero at most once, so its discriminant cannot be positive,

\begin{equation} 4R_{\bxx }^2[k] - 4R_{\bxx }[0]\sum _n x^2[n-k] \;\le \; 0 \quad \Longleftrightarrow \quad \abs {R_{\bxx }[k]} \;\le \; \sqrt {R_{\bxx }[0]\sum _n x^2[n-k]} . \end{equation}

The ACF at any lag therefore never exceeds the geometric mean of the two energies, and dividing by it gives Eq. (21.25).

Reading the value Eq. (21.24) is the sample Pearson correlation coefficient of Sec. 2.2.2, evaluated for the pair of vectors \(x[n]\) and \(x[n-k]\). The mean subtraction appearing there is unnecessary here, because the signal is assumed zero-mean. Its interpretation carries over unchanged: the sign gives the direction of the linear association and the magnitude its strength, with \(\rho _{\bxx }[k]=0\) indicating no linear dependence at that lag. Two points are specific to the lag setting:

• • Because the amplitude cancels, \(\rho _{\bxx }[k]\) is dimensionless and unaffected by the gain of the recording, so values obtained at different lags and from different records are directly comparable.

• • Its magnitude fixes the loss achievable by the one-coefficient predictor, as developed in the Model Loss part below.

Relation to the normalized ACF Under the same stationarity assumption used for the unbiased ACF (Eq. (21.19)), \(\sum _n x^2[n-k]\approx \tfrac {L-k}{L}R_{\bxx }[0]\), and

\begin{equation} \rho _{\bxx }[k] \approx \sqrt {\frac {L}{L-k}}\,R_{\bxx ,norm}[k] \approx R_{\bxx ,norm}[k] \approx h_k \qquad (k\ll L). \end{equation}

The second step drops the factor \(\sqrt {L/(L-k)}\), which tends to \(1\) only for \(k\ll L\). For a short record or a large lag this factor is not negligible, and Eq. (21.24) has to be evaluated directly.

The summary of the different ACF normalizations and their relation to the prediction coefficient \(h_k\) is given in Table 21.1. The biased and unbiased ACFs keep the scale of the signal, as their value at \(k=0\) shows, so neither is by itself an approximation of the dimensionless \(h_k\); only their ratio to that lag-\(0\) value is comparable with it.

21.1.3 ACF Properties

The signal energy (21.13) is the highest value of ACF,

\begin{equation} R_{\bxx }[0] > R_{\bxx }[k]. \end{equation}

Theoretically, \(R_{\bxx }[0] = R_{\bxx }[k]\) may happen under certain conditions but unachievable for the practical time-limited signals.

The ACF has inherent time symmetry (backward prediction),

\begin{equation} \label {eq:acf_symmetry} R_{\bx \bx }[k] =R_{\bx \bx }[-k] \end{equation}

that reduces the need to recalculation for \(k<0\).

Minimum Loss

In matrix form (by Eq. (3.21)), the loss is given by

\begin{equation} \loss _{min}(h_k) =\be ^T\be = \by ^T\by - h_k\by ^T\bx \end{equation}

The resulting loss is given by (following Eq. (3.21))

\begin{equation} \label {eq-signals-min-J} \begin{aligned} \loss _{min}(h_k) &= \sum _{n=0}^{L-1}x^2[n] - h_k\sum _{n=0}^{L-1}x[n]x[n-k]\\[3pt] &= R_{\bx \bx }[0] - h_kR_{\bx \bx }[k]. \end {aligned} \end{equation}

Since \(\sum _n x^2[n-k]\le R_{\bxx }[0]\) implies \(h_k\ge R_{\bxx ,norm}[k]\), the loss is upper-bounded by

\begin{equation} \loss _{min}(h_k) \le R_{\bx \bx }[0]\left (1-R_{\bx \bx ,norm}^2[k]\right ), \end{equation}

and, under the stationarity approximation \(R_{\bxx ,norm}[k]\approx \rho _{\bxx }[k]\) (Eq. (21.24)),

\begin{equation} \loss _{min}(h_k) \approx R_{\bx \bx }[0]\left (1-\rho _{\bxx }^2[k]\right ). \end{equation}

The correlation coefficient \(\rho _{\bxx }[k]\) is a measure of linear dependence between \(x[n]\) and \(x[n-k]\).

• • For noiseless data and a linear relation between the samples, the prediction is perfect, \(\abs {\rho _{\bx \bx }[k]}=1\) and \(\loss _{min}=0\).

• • Without any linear dependence, \(\rho _{\bx \bx }[k]=h_k=0\) and the resulting prediction is \(\hat {x}[n] = 0\).

This can be summarized as (see also (21.13) for interpretation)

\begin{equation} \label {eq-ar-loss-bounds} 0\leq \loss _{min}(h_k) \leq R_\bxx [0] \end{equation}

The illustration of the correlation coefficient principles is presented in Fig. 21.1.

MSE and RMSE The corresponding MSE and RMSE metrics are given by

\begin{align} MSE(h_k) &= \frac {1}{L}\loss _{min}(h_k)\\ RMSE(h_k) &= \sqrt {\frac {1}{L}\loss _{min}(h_k)} \end{align}

21.1.4 Auto-covariance

For simplicity, a zero-average, \(\bar {x}[n]=0\), was assumed. When the signals is non-zero mean, the subtraction of signal average (DC level) from the signal, \(x[n] = x[n]-\bar {x}[n]\) before auto-correlation calculation is termed as auto-covariance.

Note, both auto-correlation and auto-covariance have similar abbreviation, ACF. Typically, ACF is used for auto-correlation, since majority of the signals are zero-mean.

Tip:

• • ACF calculation may be significantly speed-up with appropriate algorithms and the bounded maximum value of \(k\le k_{max}\). The value of \(k_{max}\) may be decided by Eq. (21.22).

21.2 Power Spectral Density (PSD)

21.2.1 Relation between ACF and PSD

From the signal processing point of view, the interpretation of a DFT of some general random signal is non-trivial. For example, phases \(\theta _k\) of some random origin will be quite different for each time-segment (or realization) of the signal. The Wiener–Khinchin theorem states that the PSD (pages (page for section 20.4.2) and (page for section 20.6.3)) of stationary random signal is the Fourier transform of the ACF. For the raw ACF \(R_{\bxx }[m]\), with \(m\) denoting the lag,

\begin{equation} \label {eq:wiener-khinchin} S_\bxx [k] = \DFTtr {R_{\bxx }[m]} = \abs {\DFTtr {x[n]}}^2,\quad k=0,\ldots ,N-1 \end{equation}

Transforming the biased ACF instead (Eq. (21.14)) scales the result by \(1/L\) and gives the raw periodogram \(\hat {S}_\bxx [k]=\abs {\DFTtr {x[n]}}^2/L\) used in Sec. 20.7.4; the two differ only by that constant.

Important properties are

• • \(S_\bxx [k]\in \Re \), due to the symmetry (Eq. (21.31)) of \(R_\bxx [k]\) it results \(\theta _k=0\forall k\).

• • \(S_\bxx [k]\ge 0\quad \forall k\).

This relation shows, that random signal also includes spectral interpretation.

Interpretation

Peaks The PSD shows where (in frequency) the signal has most of its energy. Peaks in the PSD correspond to frequencies at which the signal exhibits strong periodic or quasi-periodic components.

Periodicity Slowly decaying (long-memory) correlations in the time domain often translate into a PSD that has more energy at low frequencies (indicating slow variations in time), as in the center column of Fig. 21.3. Conversely, if the ACF shows periodicity, the PSD will have distinct peaks at the corresponding harmonic frequencies, as in the right column.

Time-shifts The ACF reveals how similar a signal is to itself at different time shifts. A slowly decaying ACF indicates strong long-term correlations, while a rapidly decaying ACF suggests only short-range predictability.

21.2.2 Frequency estimation via ACF peak

For a (quasi-)periodic signal, the ACF \(R_{\bxx }[k]\) inherits the period of \(x[n]\): its largest non-trivial peak appears at the lag corresponding to one signal period. This provides a time-domain alternative to the periodogram-peak estimator of Sec. 20.4.

Using the biased ACF \(R_{\bxx ,biased}[k]\) (Eq. (21.14)), the estimator is

\begin{equation} \label {eq:freq-acf-estimator} \hat {k}_0 = \arg \max _{k\ge k_{\min }} R_{\bxx ,biased}[k],\quad \hat {F}_0 = \frac {F_s}{\hat {k}_0},\quad \hat {\omega }_0 = \frac {2\pi }{\hat {k}_0}. \end{equation}

The lower bound \(k_{\min }\) excludes the trivial peak at \(k=0\) (where \(R_{\bxx }[0]=E_\bx \)) and the surrounding main lobe. A practical choice is to set \(k_{\min }\) just above the ACF main-lobe width, or to start the search at the first zero-crossing (or first local minimum) of \(R_{\bxx }[k]\).

Why the biased ACF The biased form is preferred because the unbiased ACF \(R_{\bxx ,unbiased}[k]\) divides by \(L-k\) (Eq. (21.18)) and so inflates the variance for large \(k\), producing spurious peaks near the tail. The biased ACF uses the constant denominator \(L\), trading a small downward bias for a much better-behaved peak search.

This is the missing fundamental case, and it is why the autocorrelation method is the canonical pitch-detection approach in speech and audio processing, where the fundamental period is read off as \(\hat {k}_0\). Selecting a peak in the frequency domain instead risks returning a harmonic, an error of a whole octave for the \(2F_0\) line.

Tip:

• • Sub-sample lag refinement: fit a parabola to \(R_{\bxx ,biased}[\hat {k}_0-1],R_{\bxx ,biased}[\hat {k}_0],R_{\bxx ,biased}[\hat {k}_0+1]\) and use its analytic vertex as a fractional-lag estimate, analogous to parabolic peak of the periodogram in (20.48).

21.3 AR(p) Model

The auto-regressive (AR) signal model describes a signal \(y[n]\) as a linear combination of its \(p\) previous samples plus noise. Formally, an \(AR(p)\) model is

\begin{equation} \begin{aligned} y[n] &=h_1y[n-1] + h_2y[n-2] + \cdots + h_{p}y[n-p] + \epsilon [n]\\ &=\sum _{m=1}^p h_my[n-m]+ \epsilon [n] \end {aligned} \end{equation}

and the prediction is of the form

\begin{equation} \hat {y}[n] =\sum _{m=1}^p h_my[n-m] \end{equation}

where \(h_1, \ldots ,h_p\) are the AR model coefficients and \(p\) is the model order that is chosen as hyper-parameter.

This model can be easily formulated in the matrix form by using \(L\) sample of \(y[n]\), by

\begin{equation} \label {eq:ar_prediction} \hat {\by }= \begin{bmatrix} \hat {y}[1] \\ \hat {y}[2]\\ \vdots \\ \hat {y}[L-2] \\ \hat {y}[L-1] \end {bmatrix},\; \bX =\begin{bmatrix} y[0] & 0 & \cdots & 0\\ y[1] & y[0] & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ y[L-3] & y[L-4] & \cdots & y[L-p-2]\\ y[L-2] & y[L-3] & \cdots & y[L-p-1]\\ \end {bmatrix},\; \bh =\begin{bmatrix} h_1 \\ h_2 \\ \vdots \\ h_{p} \end {bmatrix} \end{equation}

where the corresponding dimensions are \(\hat {\by }\in \Re ^{L-1},\bX \in \Re ^{L-1\times p},\bh \in \Re ^p\).

The AR coefficients \(\bh \) can be found by a LS regression that minimizes the MSE loss

\begin{equation} \mathcal {L}(h_i) = \sum _n \left (y[n]-\hat {y}[n]\right )^2 = \norm {\by -\bX \bh }^2. \end{equation}

The LS solution is straightforward,

\begin{equation} \label {eq:normal_eq_matrix} \bh = {\left (\bX ^T\bX \right )^{-1}}\bX ^T\by . \end{equation}

Note, the general approach for evaluation of \(\bh \) is presented in the following section.

21.3.1 Yule-Walker Form

The resulting \(\bX ^T\bX \) matrix is termed (Toeplitz) auto-correlation matrix and can be interpreted in terms of auto-correlation values,

\begin{equation} \begin{aligned} \bR &= \bX ^T\bX \\ &= \begin{bmatrix} R_{\by \by }[0] & R_{\by \by }[1] & \cdots & R_{\by \by }[p-1]\\ R_{\by \by }[1] & R_{\by \by }[0] & \cdots & R_{\by \by }[p-2]\\ \vdots & \vdots & \ddots & \vdots \\ R_{\by \by }[p-1] & R_{\by \by }[p-2] & \cdots & R_{\by \by }[0] \end {bmatrix}, \end {aligned} \end{equation}

where unscaled (or biased) \(R_{\by \by }[k]\) is defined above. It is common practice to ignore the changes in the vector lengths in calculating (biased) ACF for the same time-shift for \(p\ll L\), e.g. diagonal matrix values are

\begin{equation} R_\byy [0] = \sum _{n=0}^{L-1}y^2[n] \approx \sum _{n=0}^{L-\textbf {2}}y^2[n] \end{equation}

The vector \(\bX ^T\by \) is also comprised of the corresponding \(R_{\byy }[k]\) values,

\begin{equation} \br = \bX ^T\by = \begin{bmatrix} R_\byy [1]\\ R_\byy [2]\\ \vdots \\ R_\byy [p]\\ \end {bmatrix} \end{equation}

Using these formulations, the resulting coefficients are given by a solution of a set of linear equations,

\begin{equation} \bR \bh = \br \end{equation}

and the result is similar to Eq. (21.45),

\begin{equation} \bh = \bR ^{-1}\br \end{equation}

The use of a \(R_{\bx \bx ,biased}[k]\) (biased ACF, Eq. (21.14)) guaranties the numerically stable solution [15, Theorem 6.1].

Squared error The resulting loss (squared error) is given by (following Eqs. (3.21) and (21.33))

\begin{equation} \begin{aligned} \loss _{min} &= R_\bxx [0] - \sum _{k=1}^p h_kR_\bxx [k] \end {aligned} \end{equation}

Theoretically, higher value of \(p\) results in lower loss in the presence of the sufficiently long correlation time (Eq. (21.21)). Note, the accuracy drops for high values of \(p\) to due to reduced accuracy of \(R_\bxx [p]\).

Biased signal The presence of a non-zero average (bias) in the signal modifies the AR model formulation. Instead of assuming a zero-mean process, a constant bias term \(\mu \) is introduced,

\begin{equation} \label {eq:ar_biased} \begin{aligned} \hat {y}[n] &=\mu + h_1y[n-1] + h_2y[n-2] + \cdots + h_{p}y[n-p] + \epsilon [n] \end {aligned} \end{equation}

In matrix form, this requires adding a column of ones \(\bOne _M\) to the data matrix \(\bX \), similarly to how it is done in multivariate LS (Eq.(3.7)).

Tips:

• • For computational and memory efficiency, the Levinson-Durbin algorithm \(O\left (L^2\right )\) is often used to solve for AR coefficients ¹, instead of direct \(O\left (L^3\right )\).

• • The AR parameters \(\bh \) are also referred to as linear prediction coefficients (LPC), emphasizing their role in predicting the current value from past samples.

• • The auto-correlation matrix \(\bR \) is guaranteed to be positive-definite (and, therefore, invertible) for ordinary or biased definitions (normalization constants reduce each other in (21.45)).

• • Theoretically, the AR model coefficients can be regularized to prevent overfitting.

• • Although this section focuses on linear AR models, nonlinear variants also exist, allowing modeling of more complex dynamics.

• • Note, there are other methods for LPC calculation, such as Burg’s method, as it is done in librosa.lpc.

21.3.2 Moving Average Filter

A simple special case of an AR-like model is the moving average filter, where all coefficients are equal and sum to one, \(h_i=\dfrac {1}{p}\). In this case,

\begin{equation} \label {eq-moving-average} \begin{aligned} \hat {y}[n] &= \frac {1}{p}\left (y[n-1] + \cdots + y[n-p]\right )\\ &=\frac {1}{p}\sum _{m=1}^p y[n-m] \end {aligned} \end{equation}

This is not strictly an AR model but shares the idea of using past samples. It smooths the signal by averaging the most recent \(p\) values.

21.3.3 Nearest Neighbor (Naïve)

A simple baseline is to use the immediate past sample as the prediction with \(h_1=1\),

\begin{equation} \label {eq-ts-1nn} \hat {y}[n] = y[n-1] \end{equation}

This "1-nearest neighbor" approach ignores signal dynamics and history. While often not very accurate, it serves as a useful baseline to compare against more sophisticated models.

21.4 Linear Prediction of Sinusoidal Signal

Consider the noise-free sinusoidal signal

\begin{equation} y[n] = \cos (\omega _0 n) \end{equation}

For simplicity, phase is initially taken as zero. We aim to model it with

\begin{equation} \hat {y}[n+1] = h_0y[n] +h_1y[n-1] \end{equation}

The corresponding loss is

\begin{equation} \begin{aligned} \loss (h_0,h_1) &= \sum _n \left (y[n+1]-\hat {y}[n+1]\right )^2\\ &= \sum _n \left (y[n+1]-h_0y[n]-h_1y[n-1]\right )^2\\ \end {aligned} \end{equation}

The required minimum is given by a solution of a system of normal equations,

\begin{equation} \begin{cases} \dfrac {\partial }{\partial h_0} \loss (h_0,h_1) = 0\\[9pt] \dfrac {\partial }{\partial h_1} \loss (h_0,h_1) = 0\\ \end {cases} \end{equation}

The resulting equations are

\begin{equation} \begin{cases} \displaystyle 2\sum _n \left (y[n+1] - h_0y[n] - h_1y[n-1]\right )\cdot (-y[n]) = 0\\[16pt] \displaystyle 2\sum _n \left (y[n+1] - h_0y[n] - h_1y[n-1]\right )\cdot (-y[n-1]) = 0\\ \end {cases} \end{equation}

These equations involve sums of trigonometrical functions. Using standard trigonometric identities and assuming a large number of samples \(L\), these sums simplify due to the oscillatory nature of sine and cosine functions.

Some important quantities are

\begin{equation} \begin{aligned} y[n+1]y[n] &=\cos (\omega _0[n])\cos (\omega _0[n+1])\\ &=\cos (\omega _0n)\left [\cos (\omega _0)\cos (\omega _0 n)-\sin (\omega _0)\sin (\omega _0 n)\right ]\\ &=\cos (\omega _0)\cos [2](\omega _0n)-\sin (\omega _0)\cos (\omega _0n)\sin (\omega _0n) \end {aligned} \end{equation}

For \(1/\omega _0\ll L\), the common assumption is

\begin{align*} \sum _n \cos [2](\omega _0n) &\approx \frac {L}{2}\\ \sum _n \cos (\omega _0n)\sin (\omega _0n) &\approx 0\\ \sum _n y[n+1]y[n] &\approx \frac {L}{2}\cos (\omega _0)\\ &\approx \sum _n y[n]y[n-1]\\ \sum _n y^2[n] &=\sum _n \cos [2](\omega _0n)\\ &\approx \sum _n y^2[n-1] \approx \frac {L}{2}\\ \sum _n y[n+1]y[n-1] &=\frac {1}{2}\sum _n \cos (2\omega _0n) + \frac {1}{2}\sum _n \cos (2\omega _0) \\ &\approx \frac {L}{2}\cos (2\omega _0) \end{align*} The approximations show that sum terms reduce to manageable forms. Substituting these approximations into the normal equations leads to a system

\begin{equation} \begin{aligned} h_0\frac {L}{2} &+ h_1\frac {L}{2}\cos (\omega _0) &&= \frac {L}{2}\cos (\omega _0)\\ h_0\frac {L}{2}\cos (\omega _0) &+h_1\frac {L}{2} &&= \frac {L}{2}\cos (2\omega _0) \end {aligned} \end{equation}

Solving these two linear equations for prediction coefficients yields

\begin{equation} \begin{aligned} h_0 &= 2\cos (\omega _0)\\ h_1 &= -1 \end {aligned} \end{equation}

This result is exact in the idealized scenario of an infinite, noise-free sinusoid.

If the signal had a phase \(\theta \),

\begin{equation} y[n] = \cos (\omega _0 n+\theta ) \end{equation}

the same coefficients still apply, only the initial conditions (the first two samples) differ,

\begin{equation} \begin{aligned} y[0] &=\cos (\theta )\\ y[1] &=\cos (\omega _0+\theta ). \end {aligned} \end{equation}

The AR(2) model still perfectly represents the sinusoidal sequence. The example of the resulting signal is presented in Fig. 21.8.

21.5 Partial auto-correlation function

The partial autocorrelation at \(k\) is the correlation that results after removing the effect of any correlations due to the terms at shorter lags,

\[ \underline {x[0]},\underbrace {x[1],x[2],\ldots ,x[j-1]}_{\text {partial out}},\underline {x[j]},x[j+1],\ldots \]

The PACF at lag \(k\) can be extracted by fitting an AR(k) model and observing the coefficient associated with \(y[n-k]\) in this model². For each \(k\),

\begin{equation} \hat {y}_{k}[n] = \phi _{k,1} y[n-1] + \phi _{k,2} y[n-2] + ... + \phi _{k,k} y[n-k] + \epsilon [n] \end{equation}

The \(k\)-th partial autocorrelation, \(\beta [k]\) is \(\phi _{k,k}\), the coefficient of \(y[n-k]\) in the AR(k) model. Each of these AR(k) models can be solved using standard LS methods.

The algorithm is as follows:

• • For the first value of PACF, \(\beta [1]\), fit AR(1) model

  \(\seteqnumber{0}{}{70}\)

  \begin{equation} \hat {y}_{1}[n] = \phi _{1,1} y[n-1] + \epsilon [n] \end{equation}

  and the coefficient \(\beta [1]=\hat {\phi }_{1,1}\) is given by the model solution, \(h_1\) (Eq. (21.5)).

• • For the second order PACF, it would be the \(\beta [2]=\phi _{2,2}\) coefficient of AR(2) model,

  \(\seteqnumber{0}{}{71}\)

  \begin{equation} \hat {y}_2[n] = \phi _{2,1} y[n-1] + \phi _{2,2} y[n-2] + \epsilon [n] \end{equation}

• • Continue this process up to the desired lag, each time extracting the coefficient of the highest-order lag as the PACF value.

Tip:

• • While the ACF reflects both direct and indirect correlations (where an earlier lag may influence a later lag through intermediate values), the PACF isolates the direct contribution of each individual lag once all shorter-term effects are factored out.

• • Though conceptually the PACF is obtained by fitting multiple AR models, practically more efficient algorithms like the Levinson-Durbin recursion can compute these values quickly and without having to solve a new LS problem at every step.

21.5.1 Relation between PACF and AR(p)

The order \(p\) of AR(p) model is related to the statistically significant (over a confidence bound) coefficients of PACF. Observing where the PACF cuts off helps identify the appropriate order \(p\) for AR(p) model fitting. An example of a synthetic signal analysis with \(p=2\) is presented in Fig. 21.9.

21.6 MA model

Another model is the moving average model (MA), where the output is a linear combination of the noise values at the different times [12, Example 4.3, pp. 90]

\begin{equation} y[n] = \epsilon [n] + c_1\epsilon [n-1] + \cdots + c_{q}\epsilon [n-q] \end{equation}

Because the noise terms \(\epsilon [n]\) are not directly observable (unlike the past outputs in AR modeling), deriving a closed-form solution for MA parameters through simple linear regression is not as direct and is not presented here. ³ While the underlying theory is well-established, practical estimation of MA parameters often involves advanced numerical methods rather than a simple closed-form solution.

21.6.1 MA and AR relations

There exists a duality between AR and MA processes in terms of infinite expansions:

Presenting AR(p) as MA(\(\infty \)) An AR model can be represented as an infinite MA model when the autoregressive parameters are stable. Consider the recursive simple AR(1),

\begin{equation} \begin{aligned} y[n] &= h_1y[n-1] + \epsilon [n]\\ &= h_1(h_1y[n-2] + \epsilon [n-1]) + \epsilon [n]\\ &= h_1^2y[n-2] + h_1\epsilon [n-1] + \epsilon [n]\\ &= h_1^3y[n-3] + h_1^2\epsilon [n-2] + h_1\epsilon [n-1] + \epsilon [n] \end {aligned} \end{equation}

Continuing this process indefinitely and assuming stability, \(h_1<1, \lim \limits _{k\rightarrow \infty }h_1^k\rightarrow 0\), we have

\begin{equation} \begin{aligned} y[n] &= \epsilon [n] + h_1\epsilon [n-1] + h_1^2\epsilon [n-2] + h_1^3\epsilon [n-3] + \cdots \\ &=\sum _{i=0}^\infty h_1^i\epsilon [n-i] \end {aligned} \end{equation}

that it is MA(\(\infty \)) model, where the coefficients of the MA representation are the infinite geometric sequence powers of \(h_1\).

Presenting MA(1) as AR(\(\infty \)) An MA(1) process is given by

\begin{equation} x[n] = \epsilon [n] + c_1\,\epsilon [n-1]. \end{equation}

We can express the process as an infinite-order AR\((\infty )\) model.

First, note that we can write

\begin{align*} \epsilon [n] &= x[n] - c_1\,\epsilon [n-1].\\ \epsilon [n-1] &= x[n-1] - c_1\,\epsilon [n-2]\\ \epsilon [n-2] &= x[n-2] - c_1\,\epsilon [n-3] \end{align*} and so on. Substituting for \(\epsilon [n-i]\) recursively, we have

\begin{equation} \begin{aligned} \epsilon [n] &= x[n] - c_1\left (x[n-1] - c_1\,\epsilon [n-2]\right )\\[1mm] &= x[n] - c_1\,x[n-1] + c_1^2\,\epsilon [n-2]\\[1mm] &= x[n] - c_1\,x[n-1] + c_1^2\left (x[n-2] - c_1\,\epsilon [n-3]\right )\\[1mm] &= x[n] - c_1\,x[n-1] + c_1^2\,x[n-2] - c_1^3\,\epsilon [n-3]\\[1mm] &\,\,\,\vdots \\[1mm] &= \sum _{i=0}^\infty (-c_1)^i\,x[n-i]. \end {aligned} \end{equation}

Rearranging, the AR\((\infty )\) representation becomes

\begin{equation} x[n] = \epsilon [n] + \sum _{i=1}^\infty (-c_1)^i\,x[n-i]. \end{equation}

Note the invertibility condition, \(|c_1| < 1\).

Interpretation AR and MA are not mutually exclusive categories. A stable AR process can be seen as a special case of an infinite MA, and a stable MA can be thought of as an infinite AR.

21.6.2 The relation between MA(q) and ACF

The MA(q) model has a distinctive fingerprint in terms of its ACF; it is nonzero for up to lag \(q\) and essentially zero afterward (except for sampling and noise effects). Just as the partial autocorrelation function (PACF) helps determine the order \(p\) of an AR(p) process by pinpointing where its PACF cuts off (Sec. 21.5.1), the ACF helps identify the order \(q\) of an MA(q) process.

An example of a synthetic MA(4) signal analysis is presented in Fig. 21.10. After lag 4, the ACF values remain within the confidence bounds, essentially zero, suggesting the data arise from an MA(4) process.

21.7 ARMA

The ARMA(p,q) model is given by

\begin{equation} \begin{aligned} \hat {y}[n] &= h_1y[n-1] + \cdots + h_py[n-p] \\ &+ c_1\epsilon [n-1] + \cdots + c_q\epsilon [n-q] + \epsilon [n]\\ &=\sum _{i=1}^{p}h_iy[n-i] + \sum _{k=0}^qc_k\epsilon [n-k] \end {aligned}, \end{equation}

where \(c_0=1\) by definition.

One of the ways to describe MA part of ARMA, is that MA uses to model the difference unexplained by AR model,

\begin{equation} \hat {y}[n] - \sum _{i=1}^{p}h_iy[n-i] = \sum _{k=0}^qc_k\epsilon [n-k] \end{equation}

The MA part is not merely a modeling convenience. The following example shows a routine measurement setup in which it appears by itself, with no choice involved.

Three readings of Eq. (21.90) are worth noting.

• • The AR coefficient survives untouched. The sensor does not change the dynamics of the state, it only contributes the MA term, whose coefficient \(c_1\) carries the sign opposite to \(h_1\).

• • The limits behave as they should. For \(\sigma _v^2\to 0\) the correlation \(\rho _{\bu \bu }[1]\to 0\) and \(c_1\to 0\), recovering the AR(1) of Eq. (21.81). For \(\sigma _v^2\to \infty \) the correlation tends to \(-h_1/(1+h_1^2)\) and \(c_1\to -h_1\), so the AR and MA parts cancel and \(z[n]\) is white, as expected when the record is all sensor noise.

• • The order of the MA part counts the noisy measurements, not the states. One noisy sensor adds one MA coefficient.

One realization of this setup, with \(h_1=0.9\), \(\sigma _\epsilon ^2=1\) and \(\sigma _v^2=2\), is shown in Fig. 21.11. The sensor noise adds a fast jitter on top of the slow swings of the hidden state, and only the jittered record is available to the model. Its AR(1) difference \(u[n]=z[n]-h_1z[n-1]\) measured on the same record follows Eq. (21.84): the estimated correlation coefficient is \(\hat {\rho }_{\bu \bu }[1]=-0.398\) against the predicted \(-0.390\), and beyond lag \(1\) it scatters about zero, never exceeding \(0.042\) in absolute value. This is the MA(1) fingerprint of Sec. 21.6.2, and it is what makes the MA term of Eq. (21.90) necessary.

As with AR models (Eq. (21.55)), a constant bias \(\mu \) can be included,

\begin{equation} \hat {y}[n] =\sum _{i=1}^{p}h_iy[n-i] + \sum _{k=0}^qc_k\epsilon [n-k] + \mu \end{equation}

21.7.1 Order Selection

Both \(p\) and \(q\) are hyper-parameters. They are fixed by fitting a grid of candidate pairs and comparing the resulting losses.

Tips:

• • The ACF and PACF fingerprints used for the pure models (Sec. 21.5.1 and Sec. 21.6.2) do not identify an ARMA model. Both of its tails decay without cutting off, so neither plot marks an order, and the two-dimensional search of Table 21.2 is used instead.

• • Fitting the largest \((p,q)\) affordable is never penalized by the in-sample loss. Comparing models of different order calls for a criterion that charges for parameters, or for a held-out record.

21.7.2 Prediction Horizon

Every fit so far was compared one step ahead. A prediction \(\ell \) steps ahead uses \(y[n],y[n-1],\ldots \) and nothing later, so the model is iterated forward one sample at a time,

\begin{equation} \label {eq:arma-horizon-rule} \hat {y}[m] = \sum _{i=1}^{p}h_i\hat {y}[m-i] + \sum _{k=1}^{q}c_k\hat {\epsilon }[m-k], \qquad m=n+1,\ldots ,n+\ell , \end{equation}

where the terms of the two sums are

\begin{equation} \label {eq:arma-horizon-sub} \hat {y}[j] = y[j],\quad j\le n \qquad \qquad \hat {\epsilon }[j] = \begin{cases} \epsilon [j], & j\le n\\ 0, & j>n \end {cases} \end{equation}

A sample past the origin \(n\) has not been observed, so the model uses its own forecast of it, supplied by the recursion itself. An innovation past the origin is white and independent of the record, so its best estimate is \(0\). These two rules are the whole of multi-step prediction, and reading them per model class is what separates the classes.

• • MA(\(q\)) forgets after \(q\) steps. Its sum has no \(\hat {y}\) terms at all, and once \(\ell >q\) every index \(m-k\) is past the origin, so every term is zeroed and the forecast is the process mean. Its memory is \(q\) samples, no fit can extend it, and beyond that its error is the full power \(\sigma _y^2\) of the signal.

• • AR(\(p\)) decays but never stops. The recursion feeds on its own output, so a forecast exists at every horizon. It decays with the poles of the model, geometrically as \(h_1^\ell \) for AR(1), and the error climbs toward \(\sigma _y^2\) gradually rather than reaching it at a fixed lag.

• • ARMA(\(p,q\)) uses both. The MA part corrects the first \(q\) steps and the AR recursion carries everything after that. Its price is that Eq. (21.94) needs \(\epsilon [j]\) for \(j\le n\), which the record does not hold: the innovations are reconstructed by the AR(\(\infty \)) expansion of Sec. 21.6.1, the first stage of the two-stage fit.

Tips:

• • The horizon is chosen by the application, and it decides the model class before any order is fitted: a model with no AR part cannot serve a horizon longer than its own order.

• • Two models that tie at \(\ell =1\) can still differ at \(\ell =5\), so a model meant to forecast far ahead is compared at the horizon it will serve.

21.7.3 Rational PSD (*)

Written in the \(a_m\) notation of Eq. (21.93), the two coefficient sets define a pair of frequency functions on the DFT grid \(\omega _k=2\pi k/N\),

\begin{equation} \label {eq:arma-poly} A(\omega _k) = \sum _{m=0}^p a_m\exp (-j\omega _km),\qquad C(\omega _k) = \sum _{m=0}^q c_m\exp (-j\omega _km), \end{equation}

that is, the DTFT of each coefficient sequence. The driving noise is white, so its PSD is flat (Sec. 21.2 and the left column of Fig. 21.3), and the PSD of the ARMA output is

\begin{equation} \label {eq:arma-psd} S_\byy [k] = \sigma _\epsilon ^2\frac {\abs {C(\omega _k)}^2}{\abs {A(\omega _k)}^2}. \end{equation}

The PSD is thus a ratio of two polynomials, which is what the ARMA coefficients buy in the frequency domain.

• • The roots of the denominator, contributed by the AR coefficients and termed poles, drive \(S_\byy [k]\) up near their frequency and produce resonance peaks.

• • The roots of the numerator, contributed by the MA coefficients and termed zeros, drive it down and produce notches.

• • A root pair of radius \(r\) at normalized frequency \(f_0\) corresponds to the coefficient triple \(\left (1,-2r\cos (2\pi f_0),r^2\right )\). Its angle sets the frequency, and its radius, approaching \(1\), sets how sharp the peak or how deep the notch becomes. A zero of radius exactly \(1\) gives an exact null, but sits on the boundary of the invertibility condition met above.

• • An AR(\(p\)) model is all-pole, \(C(\omega _k)=1\). It builds peaks directly and reproduces a notch only by spending several poles on it, which is the frequency-domain form of the order trade-off of Table 21.2. Symmetrically, an MA(\(q\)) model is all-zero and cannot produce a sharp resonance.

This combination is why ARMA is the natural model for signals whose spectra carry both resonances and anti-resonances, whereas an AR model of comparable order captures only the resonances.