Time Series Modeling - ARMA

Read Post

In this post, we will continue to cover Kitagawa G. (2021) Introduction to Time Series Modeling focusing on ARMA.

ARMA Model

A model that expresses a time series yny_{n} as a linear combination of past observations yniy_{n−i} and white noise vniv_{n−i} is called an autoregressive moving average model (ARMA model):

yn=i=1maiyni+vni=1lbivni\begin{aligned} y_{n} &= \sum_{i = 1}^{m}a_{i}y_{n-i} + v_{n} - \sum_{i = 1}^{l}b_{i}v_{n-i} \end{aligned}

The AR order and the MA order (m, l) taken together are called the ARMA order.

Further we assume the error follows a normal distribution with following properties:

E[vn]=0E[vn2]=σ2\begin{aligned} E[v_{n}] &= 0\\ E[v_{n}^{2}] &= \sigma^{2} \end{aligned} E[vnvm]=0,nm\begin{aligned} E[v_{n}v_{m}] &= 0, n \neq m \end{aligned} E[vnym]=0,n>m\begin{aligned} E[v_{n}y_{m}] &= 0, n > m \end{aligned}

It should be noted that almost all analysis of stationary time series can be achieved by using AR models.

Impulse Response Function

Using the backshift operator BB, we can express the ARMA model as:

(1i=1maiBi)yn=(1i=1lbiBi)vna(B)yn=b(B)vnyn=a(B)1b(B)vn=1i=1lbiBi1i=1maiBivn\begin{aligned} \Big(1 - \sum_{i =1}^{m}a_{i}B^{i}\Big)y_{n} &= \Big(1 - \sum_{i = 1}^{l}b_{i}B^{i}\Big)v_{n}\\[5pt] a(B)y_{n} &= b(B)v_{n}\\[5pt] y_{n} &= a(B)^{-1}b(B)v_{n}\\[5pt] &= \frac{1 - \sum_{i = 1}^{l}b_{i}B^{i}}{1 - \sum_{i =1}^{m}a_{i}B^{i}}v_{n} \end{aligned}

And define g(B)g(B) as an infinite linear combination of white noise:

g(B)=a(B)1b(B)=1i=1lbiBi1i=1maiBi=i=0giBi\begin{aligned} g(B) &= a(B)^{-1}b(B)\\ &= \frac{1 - \sum_{i = 1}^{l}b_{i}B^{i}}{1 - \sum_{i =1}^{m}a_{i}B^{i}}\\ &= \sum_{i=0}^{\infty}g_{i}B^{i} \end{aligned}

Hence:

yn=i=0giBivn\begin{aligned} y_{n} &= \sum_{i=0}^{\infty}g_{i}B^{i} v_{n} \end{aligned}

The coefficients {gi;i=0,1,}\{g_{i}; i = 0, 1, \cdots\}, correspond to the influence of the noise at time n=0n = 0 to the time series at time ii, and are called the impulse response function of the ARMA model.

Here the impulse response gig_{i} is obtained by the following recursive formula:

g0=1\begin{aligned} g_{0} &= 1 \end{aligned} gi=j=1iajgijbi\begin{aligned} g_{i} &= \sum_{j=1}^{i}a_{j}g_{i-j} - b_{i} \end{aligned} i=1,2,\begin{aligned} i = 1, 2, \cdots \end{aligned} aj=0,j>m\begin{aligned} a_{j} &= 0, j > m \end{aligned} bj=0,j>l\begin{aligned} b_{j} &= 0, j > l \end{aligned}

Example: AR(1), AR(2), MA(2), and ARMA(2, 2) Models

Consider the following models:

yn=0.9yn1+vn\begin{aligned} y_{n} &= 0.9 y_{n-1} + v_{n} \end{aligned} yn=0.93yn10.81yn2+vn\begin{aligned} y_{n} &= 0.9\sqrt{3} y_{n-1} -0.81 y_{n-2} + v_{n} \end{aligned} yn=vn0.92vn1+0.81vn2\begin{aligned} y_{n} &= v_{n} - 0.9\sqrt{2} v_{n-1} + 0.81 v_{n-2} \end{aligned} yn=0.93yn10.81yn2+vn0.92vn1+0.81vn2\begin{aligned} y_{n} &= 0.9\sqrt{3} y_{n-1} -0.81 y_{n-2} + v_{n} - 0.9\sqrt{2} v_{n-1} + 0.81 v_{n-2} \end{aligned}

The impulse response function:

dat1 <- tibble(
  g = c(
    1,
    ARMAtoMA(
      ar = 0.9,  
      ma = 0,
      (n-1)
    )
  ),
  k = 1:n
)

dat2 <- tibble(
  g = c(
    1,
    ARMAtoMA(
      ar = c(0.9 * sqrt(3), -0.81),   
      ma = 0,
      (n-1)
    )
  ),
  k = 1:n
)

dat3 <- tibble(
  g = c(
    1,
    ARMAtoMA(
      ar = 0,
      ma = c(-0.9 * sqrt(2), 0.81),
      (n - 1)
    )
  ),
  k = 1:n
)

dat4 <- tibble(
  g = c(
    1,
    ARMAtoMA(
      ar = c(0.9 * sqrt(3), -0.81),    
      ma = c(-0.9 * sqrt(2), 0.81),
      (n - 1)
    )
  ),
  k = 1:n
)

The plots above show the impulse response functions obtained from the recursive relation shown previously.

The impulse response function of the MA model is non-zero only for the initial ll points. On the other hand, if the model contains an AR part, the impulse response function has non-zero values even for i>mi > m, although it gradually decays.

The Autocovariance Function

First multiply the ARMA model by ynky_{n-k} and take expectation:

yn=i=1maiyni+vni=1lbivniynynk=i=1maiyniynk+vnynki=1lbivniynkE[ynynk]=i=1maiE[yniynk]+E[vnynk]i=1lbiE[vniynk]\begin{aligned} y_{n} &= \sum_{i = 1}^{m}a_{i}y_{n-i} + v_{n} - \sum_{i = 1}^{l}b_{i}v_{n-i}\\ y_{n}y_{n-k} &= \sum_{i = 1}^{m}a_{i}y_{n-i}y_{n-k} + v_{n}y_{n-k} - \sum_{i = 1}^{l}b_{i}v_{n-i}y_{n-k}\\ E[y_{n}y_{n-k}] &= \sum_{i = 1}^{m}a_{i}E[y_{n-i}y_{n-k}] + E[v_{n}y_{n-k}] - \sum_{i = 1}^{l}b_{i}E[v_{n-i}y_{n-k}]\\ \end{aligned}

Taking result from the impulse function, we can derive the covariance between the time series ymy_{m} and white noise vnv_{n}:

E[vnym]=E[vni=0giBivm]=i=0giE[vnvmi]={0n>mσ2gmnnm\begin{aligned} E[v_{n}y_{m}] &= E[v_{n}\sum_{i=0}^{\infty}g_{i}B^{i}v_{m}]\\ &= \sum_{i=0}^{\infty}g_{i}E[v_{n}v_{m-i}]\\ &= \begin{cases} 0 & n > m\\ \sigma^{2}g_{m-n} & n \leq m \end{cases} \end{aligned}

Since:

mi=n\begin{aligned} m-i &= n \end{aligned} i=mn\begin{aligned} i &= m - n \end{aligned}

Substituting this into the above expectation, we derive the following recurrence:

C0=i=1maiE[yniyni]+E[vnyn]i=1lbiE[vniyn]=i=1maiCi+σ2g0i=1lbign(ni)=i=1maiCi+σ2(1i=1lbigi)\begin{aligned} C_{0} &= \sum_{i = 1}^{m}a_{i}E[y_{n-i}y_{n-i}] + E[v_{n}y_{n}] - \sum_{i = 1}^{l}b_{i}E[v_{n-i}y_{n}]\\ &= \sum_{i = 1}^{m}a_{i}C_{i} + \sigma^{2}g_{0} - \sum_{i = 1}^{l}b_{i}g_{n - (n - i)}\\ &= \sum_{i = 1}^{m}a_{i}C_{i} + \sigma^{2}\Big(1 - \sum_{i=1}^{l}b_{i}g_{i}\Big) \end{aligned} Ck=i=1maiE[yniynk]+E[vnynk]i=1lbiE[vniynk]=i=1maiCki+0i=1lbignkni=i=1maiCkiσ2i=1lbigik\begin{aligned} C_{k} &= \sum_{i = 1}^{m}a_{i}E[y_{n-i}y_{n-k}] + E[v_{n}y_{n-k}] - \sum_{i = 1}^{l}b_{i}E[v_{n-i}y_{n-k}]\\ &= \sum_{i = 1}^{m}a_{i}C_{k-i} + 0 - \sum_{i = 1}^{l}b_{i}g_{n-k - n-i}\\ &= \sum_{i = 1}^{m}a_{i}C_{k-i} - \sigma^{2}\sum_{i=1}^{l}b_{i}g_{i-k} \end{aligned}

Therefore, if the orders m and l, the autoregressive and moving average coefficients aia_{i} and bib_{i}, and the innovation variance σ2\sigma^{2} of the ARMA model are given, we can obtain the autocovariance function by first computing the impulse response function g1,,glg_{1}, \cdots, g_{l} and then obtaining the autocovariance function C0,C1,C_{0},C_{1}, · · ·.

Doing this in R by using the ARMAacf() function:

n <- 20
dat1 <- tibble(
  C = ARMAacf(
    ar = 0.9,
    ma = 0,
    n
  ),
  n = 0:n
)

dat2 <- tibble(
  C = ARMAacf(
    ar = c(0.9 * sqrt(3), -0.81),
    ma = 0,
    n
  ),
  n = 0:n
)

dat3 <- tibble(
  C = ARMAacf(
    ar = 0,
    ma = c(-0.9 * sqrt(2), 0.81),
    n
  ),
  n = 0:n
)

dat4 <- tibble(
  C = ARMAacf(
    ar = c(0.9 * sqrt(3), -0.81),
    ma = c(-0.9 * sqrt(2), 0.81),    
    n
  ),
  n = 0:n
)

In particular, the following equation for the AR model obtained by putting l=0l = 0 is called the Yule-Walker equation:

C0=i=1maiCi+σ2Ck=i=1maiCki\begin{aligned} C_{0} &= \sum_{i = 1}^{m}a_{i}C_{i} + \sigma^{2}\\ C_{k} &= \sum_{i = 1}^{m}a_{i}C_{k-i} \end{aligned}

Partial Correlation

We denote the AR coefficient aia_{i} of the AR model with order m as aima_{i}^{m}. In particular, amma_{m}^{m} is called the m-th partial autocorrelation coefficient.

It can be shown the following relation holds between aim,aim1a_{i}^{m}, a_{i}^{m-1} (see Derivation of Levinson’s Algorithm):

aim=aim1ammamim1\begin{aligned} a_{i}^{m} &= a_{i}^{m-1} - a_{m}^{m}a_{m-i}^{m-1} \end{aligned} i=1,,m1\begin{aligned} i = 1, \cdots, m-1 \end{aligned}

If the m PARCOR’s, a11,,amma_{1}^{1}, \cdots, a_{m}^{m}, are given, repeated application of above equation yields the entire set of coefficients of the AR models with orders 2 through m.

On the other hand, if the coefficients a1m,,amma_{1}^{m}, \cdots, a_{m}^{m} of the AR model of the highest order are given and by solving the following equations:

aim1=aim+ammamim1(amm)2\begin{aligned} a_{i}^{m-1} &= \frac{a_{i}^{m} + a_{m}^{m}a_{m-i}^{m}}{1 - (a_{m}^{m})^{2}} \end{aligned}

The coefficients of the AR model with order m1m-1 are obtained.

To show this, given:

aim=aim1ammamim1\begin{aligned} a_{i}^{m} &= a_{i}^{m-1} - a_{m}^{m}a_{m-i}^{m-1} \end{aligned}

And substituting ii with mim-i:

amim=aim1ammam(mi)m1aim1=amim+ammaim1\begin{aligned} a_{m-i}^{m} &= a_{i}^{m-1} - a_{m}^{m}a_{m-(m-i)}^{m-1}\\ a_{i}^{m-1} &= a_{m-i}^{m} + a_{m}^{m}a_{i}^{m-1} \end{aligned}

Substituting aim1a_{i}^{m-1}:

aim=aim1amm(amim+ammaim1)=aim1ammamim(amm)2aim1=aim1(1amm)(amm)2amim\begin{aligned} a_{i}^{m} &= a_{i}^{m-1}- a_{m}^{m}(a_{m-i}^{m} + a_{m}^{m}a_{i}^{m-1})\\ &= a_{i}^{m-1}- a_{m}^{m}a_{m-i}^{m} - (a_{m}^{m})^{2}a_{i}^{m-1}\\ &= a_{i}^{m-1}(1 - a_{m}^{m}) - (a_{m}^{m})^{2}a_{m-i}^{m} \end{aligned} aim1=aim+ammamim1(amm)2\begin{aligned} a_{i}^{m-1} &= \frac{a_{i}^{m} + a_{m}^{m}a_{m-i}^{m}}{1 - (a_{m}^{m})^{2}} \end{aligned}

The PARCORs, a11,,amma_{1}^{1}, \cdots, a_{m}^{m}, can be obtained by repeating the above computation.

This reveals that estimation of the coefficients a1m,,amma_{1}^{m}, \cdots, a_{m}^{m} of the AR model of order m is equivalent to estimation of the PARCOR’s up to the order m, i.e., a11,,amma_{1}^{1}, \cdots, a_{m}^{m}.

The following R code shows the PARCOR of the four models by using the same ARMAacf() function but having the argument pacf = TRUE:

dat1 <- tibble(
  C = ARMAacf(
    ar = 0.9,
    ma = 0,
    n,
    pacf = TRUE
  ),
  n = 1:n
)

dat2 <- tibble(
  C = ARMAacf(
    ar = c(0.9 * sqrt(3), -0.81),
    ma = 0,
    n,
    pacf = TRUE
  ),
  n = 1:n
)

dat3 <- tibble(
  C = ARMAacf(
    ar = 0,
    ma = c(-0.9 * sqrt(2), 0.81),
    n,
    pacf = TRUE
  ),
  n = 1:n
)

dat4 <- tibble(
  C = ARMAacf(
    ar = c(0.9 * sqrt(3), -0.81),
    ma = c(-0.9 * sqrt(2), 0.81),    
    n,
    pacf = TRUE
  ),
  n = 1:n
)

The Power Spectrum of the ARMA Process

Recall the power spectrum is defined as:

p(f)=k=Cke2πikf\begin{aligned} p(f) &= \sum_{k=-\infty}^{\infty}C_{k}e^{-2\pi i k f} \end{aligned}

If an ARMA model of a time series is given, the power spectrum of the time series can be obtained easily:

p(f)=k=Cke2πikf=k=E[ynynk]e2πikf=k=E[j=0gjvnjp=0gpvnkp]e2πikf=k=j=0p=0gjgpE[vnjvnkp]e2πikf\begin{aligned} p(f) &= \sum_{k=-\infty}^{\infty}C_{k}e^{-2\pi i k f}\\ &= \sum_{k=-\infty}^{\infty}E[y_{n}y_{n-k}]e^{-2\pi i k f}\\ &= \sum_{k=-\infty}^{\infty}E\Big[\sum_{j=0}^{\infty}g_{j}v_{n-j}\sum_{p=0}^{\infty}g_{p}v_{n-k-p}\Big]e^{-2\pi i k f}\\ &= \sum_{k=-\infty}^{\infty}\sum_{j=0}^{\infty}\sum_{p=0}^{\infty}g_{j}g_{p}E[v_{n-j}v_{n-k-p}]e^{-2\pi i k f} \end{aligned}

Here using gp=0g_{p} = 0 for p<0p < 0 and:

nkp=njp=nk(nj)=jk\begin{aligned} n - k - p &= n - j\\ p &= n - k - (n - j)\\ &= j - k \end{aligned}

The power spectrum is expressed as:

p(f)=σ2k=j=0gjgnk(jk)e2πikf=σ2k=j=0gjgjke2πikf=σ2k=j=0gjgjke2πijf+2πijf2πikf=σ2k=j=0gjgjke2πijf2π(kj)f=σ2k=j=0gje2πijfgjke2πi(kj)f=σ2j=0p=0gje2πijfgpe2πipf=σ2j=0gje2πijf2\begin{aligned} p(f) &= \sigma^{2}\sum_{k=-\infty}^{\infty}\sum_{j=0}^{\infty}g_{j}g_{n - k - (j-k)}e^{-2\pi i k f}\\ &= \sigma^{2}\sum_{k=-\infty}^{\infty}\sum_{j=0}^{\infty}g_{j}g_{j-k}e^{-2\pi i k f}\\ &= \sigma^{2}\sum_{k=-\infty}^{\infty}\sum_{j=0}^{\infty}g_{j}g_{j-k}e^{-2\pi i j f +2 \pi ijf - 2\pi ikf}\\ &= \sigma^{2}\sum_{k=-\infty}^{\infty}\sum_{j=0}^{\infty}g_{j}g_{j-k}e^{-2\pi ijf - 2\pi(k-j)f}\\ &= \sigma^{2}\sum_{k=-\infty}^{\infty}\sum_{j=0}^{\infty}g_{j}e^{-2\pi i j f}g_{j-k}e^{-2\pi i (k-j) f}\\ &= \sigma^{2}\sum_{j=0}^{\infty}\sum_{p=0}^{\infty}g_{j}e^{-2\pi i j f}g_{p}e^{2\pi i p f}\\ &= \sigma^{2}\Big|\sum_{j=0}^{\infty}g_{j}e^{-2\pi i j f}\Big|^{2} \end{aligned}

j=0gje2πijf\sum_{j=0}^{\infty}g_{j}e^{-2\pi i j f} is the Fourier transform of the impulse response function, and is called the frequency response function.

Focusing on the frequency response function:

j=0gje2πijf=1j=1lbje2πijf1j=1maje2πijf\begin{aligned} \sum_{j=0}^{\infty}g_{j}e^{-2\pi i j f} &= \frac{1 - \sum_{j=1}^{l}b_{j}e^{-2\pi ijf}}{1 - \sum_{j=1}^{m}a_{j}e^{-2\pi i j f}} \end{aligned}

Therefore the power spectrum can be derived from the ARMA model by:

p(f)=σ21j=1lbje2πijf21j=1maje2πijf2\begin{aligned} p(f) &= \sigma^{2}\frac{\Big|1 - \sum_{j=1}^{l}b_{j}e^{-2\pi ijf}\Big|^{2}}{\Big|1 - \sum_{j=1}^{m}a_{j}e^{-2\pi i j f}\Big|^{2}} \end{aligned}

Example

Computing the log-spectrum for the 4 models:

sp1 <- astsa::arma.spec(0.9)

sp2 <- astsa::arma.spec(
  ar = c(0.9 * sqrt(3), -0.81)
)

sp3 <- astsa::arma.spec(
  ma = c(-0.9 * sqrt(2), 0.81)
)

sp4 <- astsa::arma.spec(
  ar = c(0.9 * sqrt(3), -0.81),
  ma = c(-0.9 * sqrt(2), 0.81)
)

tibble(
  freq = sp1$freq,
  spec = log10(sp1$spec)
) |>
  ggplot(aes(x = freq, y = spec)) +    
  geom_line() +
  ggtitle("AR(1)") +
  theme_tq()

The power spectrum of the AR model with order 1 does not have any peak or trough. A peak is seen in the spectrum of the 2nd order AR model and one trough is seen in the plot of the 2nd order MA model. On the other hand, the spectrum of the ARMA model of order (2,2) shown in the plot has both one peak and one trough.

The logarithm of the spectrum is expressible as:

logp(f)=logσ22log1j=1maje2πijf+2log1j=1lbje2πijf\begin{aligned} \log p(f) &= \log \sigma^{2} - 2 \log\Big|1 - \sum_{j=1}^{m}a_{j}e^{-2\pi ijf}\Big| + 2 \log\Big|1 - \sum_{j=1}^{l}b_{j}e^{-2\pi ijf}\Big| \end{aligned}

Therefore, the peak and the trough of the spectrum appear at the local minimum of 1j=1maje2πijf\mid 1 - \sum_{j=1}^{m}a_{j}e^{-2\pi ijf}\mid and 1j=1lbje2πijf\mid 1 - \sum_{j=1}^{l}b_{j}e^{-2\pi ijf}\mid, respectively.

To express k peaks or k troughs, the AR order or the MA order must be higher than or equal to 2k, respectively. Moreover, the locations and the heights of the peaks or the troughs are related to the angles and the absolute values of the complex roots of the characteristic equation.

The Characteristic Equation

The characteristics of an ARMA model are determined by the roots of the following two polynomial equations:

a(B)=1j=1majBj=0b(B)=1j=1lbjBj=0\begin{aligned} a(B) &= 1 - \sum_{j=1}^{m}a_{j}B^{j} = 0\\ b(B) &= 1 - \sum_{j=1}^{l}b_{j}B^{j} = 0 \end{aligned}

The roots of these equations are called the characteristic roots. If the roots of the characteristic equation a(B)=0a(B) = 0 of the AR operator all lie outside the unit circle, the influence of noise turbulence at a certain time decays as time progresses, and the ARMA model becomes stationary.

On the other hand, if all roots of the characteristic equation b(B)=0b(B) = 0 of the MA operator lie outside the unit circle, the coefficient of hih_{i} of b(B)1=i=0hiBib(B)^{−1} = \sum_{i=0}^{\infty}h_{i}B^{i} converges and the ARMA model can be expressed by an AR model of infinite order as:

yn=i=1hiyni+vn\begin{aligned} y_{n} &= -\sum_{i=1}^{\infty}h_{i}y_{n-i} + v_{n} \end{aligned}

In this case, the time series is called invertible.

In particular, when the angles of the complex roots of the AR operator coincide with those of the MA operator, a line-like spectral peak or trough appears.

For example, if the AR and MA coefficients of the ARMA (2,2) model are given by:

a1=0.992a2=0.992b1=0.952b2=0.952\begin{aligned} a_{1} &= 0.99\sqrt{2}\\ a_{2} &= -0.99^{2}\\ b_{1} &= 0.95\sqrt{2}\\ b_{2} &= -0.95^{2} \end{aligned}

Both characteristic equations have roots at f=0.125f = 0.125 (0.125×460=450.125 \times 460 = 45):

> polyroot(c(1, -0.99 * sqrt(2), 0.99^2)) |>    
  Arg() * 180 / pi
[1]  45 -45

> polyroot(c(1, -0.95 * sqrt(2), 0.95^2)) |>
  Arg() * 180 / pi
[1]  45 -45

And the spectral:

sp1 <- astsa::arma.spec(
  ar = c(0.99 * sqrt(2), -0.99^2)    
)

sp2 <- astsa::arma.spec(
  ma = c(-0.95 * sqrt(2), 0.95^2)
)

sp3 <- astsa::arma.spec(
  ar = c(0.99 * sqrt(2), -0.99^2),
  ma = c(-0.95 * sqrt(2), 0.95^2)
)

The peak of the spectrum (or trough) appears at f=θ2πf = \frac{\theta}{2\pi}, if the complex root of AR (or MA) operator is expressed in the form:

z=α+iβ=reiθ\begin{aligned} z &= \alpha + i\beta\\ &= re^{i\theta} \end{aligned}

Further, the closer the root r approaches to 1, the sharper the peak or the trough of the spectrum become.

The function armachar of the package TSSS computes and draw graphs of the impulse response function, the autocovariance function, the PARCORs, power spectrum and the characteristic roots of the AR operator and the MA operator of an ARMA model.

For example the following AR, MA, and ARMA models:

yn=0.9×3yn10.81yn2+vnyn=vn0.9×2+0.81vn2\begin{aligned} y_{n} &= 0.9\times \sqrt{3}y_{n-1} - 0.81y_{n-2} + v_{n}\\ y_{n} &= v_{n} - 0.9\times \sqrt{2} + 0.81v_{n-2} \end{aligned} yn=0.9×3yn10.81yn2+vn0.9×2+0.81vn2\begin{aligned} y_{n} &= 0.9\times \sqrt{3}y_{n-1} - 0.81y_{n-2} + v_{n} - 0.9\times \sqrt{2} + 0.81v_{n-2} \end{aligned}
a <- c(0.99 * sqrt(2), -0.9801)
b <- c(0.95 * sqrt(2), -0.9025)  
fit <- armachar(
  arcoef = a,
  macoef = b,
  v = 1.0,
  lag = 20
)

> names(fit)
[1] "impuls"   "acov"     "parcor"   "spec"    
[5] "croot.ar" "croot.ma"

Note that the inverse roots are graphed instead so the points are inside rather than outside the unit circle.

The Multivariate AR Model

Similar to the case of univariate time series, for multivariate time series:

yn=[yn(1)yn(l)]\begin{aligned} \textbf{y}_{n} &= \begin{bmatrix} y_{n}(1)\\ \vdots\\ y_{n}(l) \end{bmatrix} \end{aligned}

The model that expresses a present value of the time series as a linear combination of past values yn1,,ynM\textbf{y}_{n−1}, \cdots, \textbf{y}_{n−M} and white noise is called a multivariate autoregressive model (MAR model):

yn=m=1MAmynm+v\begin{aligned} \textbf{y}_{n} &= \sum_{m=1}^{M}\textbf{A}_{m}\textbf{y}_{n-m} + \textbf{v} \end{aligned}

Where Am\textbf{A}_{m} is the autoregressive coefficient matrix whose (i, j)-th element is given by am(i,j)a_{m}(i, j), and vnv_{n} is an l-dimensional white noise that satisfies:

E[vn]=[00]\begin{aligned} E[\textbf{v}_{n}] &= \begin{bmatrix} 0\\ \vdots\\ 0 \end{bmatrix} \end{aligned} E[vnvnT]=[σ11σ1lσl1σll]=W\begin{aligned} E[\textbf{v}_{n}\textbf{v}_{n}^{T}] &= \begin{bmatrix} \sigma_{11} & \cdots & \sigma_{1l}\\ \vdots & \ddots & \vdots\\ \sigma_{l1} & \cdots & \sigma_{ll} \end{bmatrix} = \textbf{W} \end{aligned} E[vnvmT]=0n×m,nm\begin{aligned} E[\textbf{v}_{n}\textbf{v}_{m}^{T}] &= \textbf{0}_{n \times m}, n \neq m \end{aligned} E[vnymT]=0n×m,n>m\begin{aligned} E[\textbf{v}_{n}\textbf{y}_{m}^{T}] &= \textbf{0}_{n\times m}, n > m \end{aligned}

The cross-covariance function of yn(i),yn(j)\textbf{y}_{n}(i), \textbf{y}_{n}(j) is defined by:

Ck(i,j)=E[yn(i)ynk(j)]\begin{aligned} \textbf{C}_{k}(i, j)&= E[\textbf{y}_{n}(i)\textbf{y}_{n-k}(j)] \end{aligned}

For the l×ll \times l matrix:

Ck=E[ynynkT]\begin{aligned} \textbf{C}_{k} = E[\textbf{y}_{n}\textbf{y}_{n-k}^{T}] \end{aligned}

The (i, j)-th component of which is Ck(i,j)\textbf{C}_{k}(i, j), is called the cross-covariance function.

Similar to the case of the univariate time series, for the multivariate AR model, Ck\textbf{C}_{k} satisfies the Yule-Walker equation:

C0=j=1MAjCj+WCk=j=1MAjCkj\begin{aligned} \textbf{C}_{0} &= \sum_{j=1}^{M}\textbf{A}_{j}\textbf{C}_{-j} + \textbf{W}\\ \textbf{C}_{k} &= \sum_{j=1}^{M}\textbf{A}_{j}\textbf{C}_{k-j} \end{aligned} k=1,2,\begin{aligned} k = 1, 2, \cdots \end{aligned}

The cross-covariance function is not symmetric. Therefore, for multivariate time series, the Yule-Walker equations for the the backward AR model and the forward AR model are different.

The Fourier transform of the cross-covariance function Ck(s,j)\textbf{C}_{k}(s, j) is called the cross-spectral density function:

psj(f)=k=Ck(s,j)e2πikf=k=Ck(s,j)cos(2πkf)ik=Ck(s,j)sin(2πkf)\begin{aligned} p_{sj}(f) &= \sum_{k=-\infty}^{\infty}\textbf{C}_{k}(s, j)e^{-2\pi i k f}\\ &= \sum_{k=-\infty}^{\infty}\textbf{C}_{k}(s, j)\textrm{cos}(2\pi k f) - i\sum_{k=-\infty}^{\infty}\textbf{C}_{k}(s, j)\textrm{sin}(2\pi k f) \end{aligned}

Since the cross-covariance function is not an even function, the cross-spectrum given by:

P(f)=[p11(f)p1l(f)pl1(f)pll(f)]\begin{aligned} \textbf{P}(f) &= \begin{bmatrix} p_{11}(f) & \cdots & p_{1l}(f)\\ \vdots & \ddots & \vdots\\ p_{l1}(f) & \cdots & p_{ll}(f) \end{bmatrix} \end{aligned}

And has an imaginary part and is a complex number.

And the relations between P(f)\mathbf{P}(f) and the cross-covariance matrix Ck\textbf{C}_{k} are given by:

P(f)=k=Cke2πikfCk=1212P(f)e2πikf df\begin{aligned} \mathbf{P}(f) &= \sum_{k=-\infty}^{\infty}\textbf{C}_{k}e^{-2\pi i k f}\\ \textbf{C}_{k} &= \int_{-\frac{1}{2}}^{\frac{1}{2}}\mathbf{P}(f)e^{2\pi i k f}\ df \end{aligned}

For time series that follow the multivariate AR model, the cross-spectrum can be obtained by:

P(f)=A(f)1W(A(f)1)\begin{aligned} \mathbf{P}(f) &= \textbf{A}(f)^{-1}\textbf{W}(\textbf{A}(f)^{-1})^{*} \end{aligned}

Where A\textbf{A}^{*} denotes the complex conjugate of the matrix A\textbf{A}, and A(f)\textbf{A}(f) denotes the l×ll \times l matrix whose (j, k)-th component is defined by:

Ajk(f)=m=0Mam(j,k)e2πimf\begin{aligned} \textbf{A}_{jk}(f) &= \sum_{m=0}^{M}a_{m}(j, k)e^{-2\pi i m f} \end{aligned}

Here, it is assumed that:

a0(j,j)=1\begin{aligned} a_{0}(j, j) &= −1 \end{aligned} a0(j,k)=0,jk\begin{aligned} a_{0}(j, k) &= 0, j \neq k \end{aligned}

Given a frequency f, the cross-spectrum is a complex number and can be expressed as:

pjk(f)=αjk(f)eiϕjk(f)\begin{aligned} p_{jk}(f) &= \alpha_{jk}(f)e^{i\phi_{jk}(f)} \end{aligned}

Where:

αjk(f)=Re(pjk(f))2+Im(pjk(f))2ϕjk(f)=arctan(Im(pjk(f))Re(pjk(f)))\begin{aligned} \alpha_{jk}(f) &= \sqrt{\mathrm{Re}(p_{jk}(f))^{2} + \mathrm{Im}(p_{jk}(f))^{2}}\\ \phi_{jk}(f) &= \mathrm{arctan}\Bigg(\frac{\mathrm{Im}(p_{jk}(f))}{\textrm{Re}(p_{jk}(f))}\Bigg) \end{aligned}

αjk(f)\alpha_{jk}(f) is called the amplitude spectrum and ϕjk(f)\phi_{jk}(f) is the phase spectrum.

The following denotes the the square of the correlation coefficient between frequency components of time series yn(j)\textbf{y}_{n}(j) and yn(k)\textbf{y}_{n}(k) at frequency f and is called the coherency:

cohjk(f)=αjk(f)2pjj(f)pkk(f)\begin{aligned} \textrm{coh}_{jk}(f) &= \frac{\alpha_{jk}(f)^{2}}{p_{jj}(f)p_{kk}(f)} \end{aligned}

For convenience, in the following A(f)1\textbf{A}(f)^{-1} will be denoted as B(f)={bjk(f)}\textbf{B}(f) = \{b_{jk}(f)\}.

If the components of the white noise vn\textbf{v}_{n} are mutually uncorrelated and the variance covariance matrix is the diagonal matrix:

W=[σ12σl2]\begin{aligned} \textbf{W} &= \begin{bmatrix} \sigma_{1}^{2} & &\\ & \ddots &\\ &&\sigma_{l}^{2} \end{bmatrix} \end{aligned}

Then the power spectrum of the i-th component of the time series can be expressed as:

pii(f)=j=1lbij(f)σj2bij(f)=j=1lbij(f)2σj2\begin{aligned} p_{ii}(f) &= \sum_{j=1}^{l}b_{ij}(f)\sigma_{j}^{2}b_{ij}(f)^{*}\\ &= \sum_{j=1}^{l}|b_{ij}(f)|^{2}\sigma_{j}^{2} \end{aligned}

This indicates that the power of the fluctuation of the component i at frequency ff can be decomposed into the effects of ll noises.

Therefore, if we define rij(f)r_{ij}(f) by:

rij(f)=bij(f)2σj2pii(f)\begin{aligned} r_{ij}(f) &= \frac{|b_{ij}(f)|^{2}\sigma_{j}^{2}}{p_{ii}(f)} \end{aligned}

It represents the ratio of the effect of vn(j)\textbf{v}_{n}(j) in the fluctuation of yn(i)\textbf{y}_{n}(i) at frequency f. This quantity rij(f)r_{ij}(f) is called the relative power contribution.

Example:

The function marspc of the package TSSS computes the cross-spectra, the coherency and the power contribution of multivariate time series.

The following arguments are required:

  • arcoef: AR coefficient matrices.
  • v: innovation covariance matrix.

And the outputs from this function are:

  • spec: cross-spectra.
  • amp: amplitude spectra.
  • phase: phase spectra.
  • coh: simple coherency.
  • power: power contribution.
  • rpower: relative power contribution.

Using the HAKUSAN dataset:

yn=[yn(Yaw)yn(Pitch)yn(Rudder)]\begin{aligned} \textbf{y}_{n} &= \begin{bmatrix} y_{n}(\text{Yaw})\\ y_{n}(\text{Pitch})\\ y_{n}(\text{Rudder}) \end{bmatrix} \end{aligned}

We take half (N = 500) of the samples at Δt=2\Delta t = 2:

yy <- HAKUSAN |>
  spread(key = key, value = value) |>
  as_tibble() |>
  select(YawRate, Rolling, Rudder) |>
  as.matrix()

nc <- dim(yy)[1]
n <- seq(1, nc, by = 2)
y <- yy[n, ]

# Fit MAR model
fit <- marfit(y, lag = 20)

> fit

 Order     AIC
    0      7091.707
    1      6238.803
    2      5275.361
    3      5173.022
    4      5135.199
    5      5136.631
    6      5121.019
    7      5105.826
    8      5096.349
    9      5087.905
   10      5093.793
   11      5091.420
   12      5099.441
   13      5097.981
   14      5100.978
   15      5113.050
   16      5116.516
   17      5129.417
   18      5136.063
   19      5143.562
   20      5157.365

 Minimum AIC = 5087.905  attained at m = 9     

 Innovation covariance matrix
  2.225414e+00 -4.712800e-01 -3.402811e-01
 -4.712800e-01  7.129234e-01  1.063104e-01
 -3.402811e-01  1.063104e-01  2.781570e+00

> names(fit)
[1] "maice.order" "aic" "v" "arcoef"     

To retrieve the spectral, coherence and power calculations:

fit <- marspc(
  fit$arcoef,
  v = fit$v,
  plot = FALSE
)

> names(fit)
[1] "spec"   "amp"    "phase"  "coh"    "power"  "rpower"    

The power spectral for the 3 series:

n <- 200
dat <- tibble(
  yaw_spec = fit$spec[, 1, 1],
  pitch_spec = fit$spec[, 2, 2],
  rudder_spec = fit$spec[, 3, 3],
  x = seq(0, 0.5, length.out = n + 1)    
)

dat |>
  ggplot(
    aes(
      x = x,
      y = log10(Mod(yaw_spec))
    )
  ) +
  geom_line() +
  xlab("Freq") +
  ylab("Log-Power") +
  ggtitle("Yaw") +
  theme_tq()

For the power spectra of the yaw rate, the maximum peak is seen in the vicinity of f=0.25f = 0.25 (8 seconds cycle) and for the role rate and the rudder angle in the vicinity of f=0.125f = 0.125 (16 seconds cycle).

Cross Phase Spectral:

dat <- tibble(
  yaw_pitch = fit$spec[, 1, 2],
  yaw_rudder = fit$spec[, 1, 3],
  pitch_rudder = fit$spec[, 2, 3],
  x = seq(0, 0.5, length.out = n + 1)
)

dat |>
  ggplot(aes(x = x, y = log10(Mod(yaw_pitch)))) +    
  geom_line() +
  xlab("Freq") +
  ylab("Log-Power") +
  ggtitle("Cross Yaw-Pitch") +
  theme_tq()

Note that phase spectra below are displayed within the range [π,π][-\pi, \pi]:

dat |>
  ggplot(aes(x = x, y = Arg(yaw_pitch))) +    
  geom_line() +
  xlab("Freq") +
  ylab("Phase") +
  ggtitle("Cross Yaw-Pitch") +
  theme_tq()

Coherencies:

dat <- tibble(
  yaw_pitch = fit$coh[, 1, 2],
  yaw_rudder = fit$coh[, 1, 3],
  pitch_rudder = fit$coh[, 2, 3],
  x = seq(0, 0.5, length.out = n + 1)    
)

dat |>
  ggplot(
    aes(
      x = x,
      y = yaw_pitch
    )
  ) +
  geom_line() +
  xlab("Freq") +
  ylab("Coherency") +
  ggtitle("Coherency Yaw-Pitch") +
  theme_tq()

Yaw and rudder angle have two significant peaks at the same frequencies. The peak of the yaw rate has slightly smaller frequency than those of yaw rate and rolling.

Power Contributions in terms of the ll noises bij(f)2\mid b_{ij}(f)\mid^{2}:

dat <- tibble(
  yaw = fit$power[, 1, 1],
  yaw_pitch = fit$power[, 1, 2],
  yaw_rudder = fit$power[, 1, 3],
  pitch = fit$power[, 2, 2],
  pitch_yaw = fit$power[, 2, 1],
  pitch_rudder = fit$power[, 2, 3],
  rudder = fit$power[, 3, 3],
  rudder_yaw = fit$power[, 3, 1],
  rudder_pitch = fit$power[, 3, 2],
  x = seq(0, 0.5, length.out = n + 1)
)

dat |>
  ggplot(
    aes(
      x = x,
      y = yaw
    )
  ) +
  geom_line(aes(col = "Yaw")) +
  geom_line(aes(y = yaw_pitch, col = "Yaw-Pitch")) +
  geom_line(aes(y = yaw_rudder, col = "Yaw-Rudder")) +    
  xlab("Freq") +
  ylab("Power (Yaw)") +
  ggtitle("Power Contribution") +
  scale_color_manual(
    name = "",
    breaks = c(
      "Yaw",
      "Yaw-Pitch",
      "Yaw-Rudder"
    ),
    values = c(
      "Yaw" = "black",
      "Yaw-Pitch" = "blue",
      "Yaw-Rudder" = "darkgreen"
    )
  ) +
  theme_tq()

The influence of the rudder angle is clearly visible for f<0.13f < 0.13. However, the influence of the rudder angle is barely noticeable in the vicinity of f=0.14f = 0.14 and f=0.23f = 0.23, where the dominant power of the yaw rate is found. This is probably explained by noting that this data set has been observed under the control of a conventional PID autopilot system that is designed to suppress the power of variation in the frequency area f<0.13f < 0.13.

The contribution of the rudder angle is almost 50 percent in the vicinity of f=0.1f = 0.1 where the power spectrum is significant. This is thought to be a side effect of the steering to suppress the variation of the yaw rate. On the other hand, for f>0.14f > 0.14, a strong influence of the yaw rate is seen, although the actual power is very small.

It can be seen that the influence of the yaw rate is very strong in the vicinity of f=0.12f = 0.12 where the power spectrum is large. Moreover, it can also be seen that in the range of f<0.08f < 0.08, the contribution of the yaw rate becomes greater as the frequency decreases.

Fitting an AR Model

Given a time series y1,,yNy_{1}, \cdots, y_{N}, we consider the problem of fitting an autoregressive model (AR model):

yn=i=1maiyni+vn\begin{aligned} y_{n} &= \sum_{i=1}^{m}a_{i}y_{n-i} + v_{n} \end{aligned}

To identify an AR model, it is necessary to determine the order m and estimate the autoregressive coefficients a1,,ama_{1}, \cdots, a_{m} and the variance σ2\sigma^{2} from the data.

θ=[a1amσ2]\begin{aligned} \boldsymbol{\theta} &= \begin{bmatrix} a_{1}\\ \vdots\\ a_{m}\\ \sigma^{2} \end{bmatrix} \end{aligned}

σ2\sigma^{2} is often called the innovation variance.

Under the assumption that the order m is given, consider the estimation of the parameter θ\boldsymbol{\theta} by the maximum likelihood method. The joint distribution of time series:

y=[y1yN]\begin{aligned} \textbf{y} &= \begin{bmatrix} y_{1}\\ \vdots\\ y_{N} \end{bmatrix} \end{aligned}

Becomes a multivariate normal distribution.

The mean vector of the autoregressive process:

E[y]=0\begin{aligned} E[\textbf{y}] &= \textbf{0} \end{aligned}

And the variance-covariance matrix is given by:

Σ=[C0C1CN1C1C0CN2CN1CN2C0]\begin{aligned} \boldsymbol{\Sigma} &= \begin{bmatrix} C_{0} & C_{1} & \cdots & C_{N-1}\\ C_{1} & C_{0} & \cdots & C_{N-2}\\ \vdots & \vdots & \ddots & \vdots\\ C_{N-1} & C_{N-2} & \cdots & C_{0} \end{bmatrix} \end{aligned}

The likelihood of the AR model is obtained by:

L(θ)=p(y1,,yNθ)=1(2π)2Σe12yTΣ1y\begin{aligned} L(\boldsymbol{\theta}) &= p(y_{1}, \cdots, y_{N}|\boldsymbol{\theta})\\ &= \frac{1}{\sqrt{(2\pi)^{2}|\boldsymbol{\Sigma}|}}e^{-\frac{1}{2}\textbf{y}^{T}\boldsymbol{\Sigma}^{-1}\textbf{y}} \end{aligned}

However, when the number of data, N, is large, the computation of the likelihood by this method becomes difficult because it involves the inversion and computation of the determinant of the N×NN \times N matrix Σ\boldsymbol{\Sigma}. Further, to obtain the maximum likelihood estimate of θ\boldsymbol{\theta}, it is necessary to apply a numerical optimization method, since the likelihood is a complicated function of the parameter.

In general, however, the likelihood of a time series model can be efficiently calculated by expressing it as a product of conditional distributions:

L(θ)=p(y1,,yNθ)=p(yNy1,,yN1,θ)p(y1,,yN1θ)=p(yNy1,,yN1,θ)p(yN1y1,,yN2,θ)p(y1,,yN2θ)  =p(y1,θ)n=2Np(yny1,,yn1,θ)\begin{aligned} L(\boldsymbol{\theta}) &= p(y_{1}, \cdots, y_{N}|\boldsymbol{\theta})\\ &= p(y_{N}|y_{1}, \cdots, y_{N-1}, \boldsymbol{\theta})p(y_{1}, \cdots, y_{N-1}|\boldsymbol{\theta})\\ &= p(y_{N}|y_{1}, \cdots, y_{N-1}, \boldsymbol{\theta})p(y_{N-1}|y_{1}, \cdots, y_{N-2}, \boldsymbol{\theta})p(y_{1}, \cdots, y_{N-2}|\boldsymbol{\theta})\\ &\ \ \vdots\\ &= p(y_{1}, \boldsymbol{\theta})\prod_{n=2}^{N}p(y_{n}|y_{1}, \cdots, y_{n-1}, \boldsymbol{\theta}) \end{aligned}

Using a Kalman filter, each term in the right-hand side can be efficiently and exactly evaluated, which makes it possible to compute the exact likelihood of the ARMA model and other time series models.

The AIC for the model is defined by:

AIC=2logL(θ^)+2(m+1)\begin{aligned} \text{AIC} &= -2 \log L(\hat{\boldsymbol{\theta}}) + 2(m+1) \end{aligned}

Yule-Walker Method

As shown previously, the autocovariance function of the AR of order m satisfies the Yule-Walker equation:

C0=i=1maiCi+σ2Cj=i=1maiCj1\begin{aligned} C_{0} &= \sum_{i=1}^{m}a_{i}C_{i} + \sigma^{2}\\ C_{j} &= \sum_{i=1}^{m}a_{i}C_{j-1} \end{aligned}

By computing the sample autocovariance functions C^k\hat{C}_{k} and substituting them into above equations, we obtain a system of linear equations for the unknown autoregressive coefficients a1,,ama_{1}, \cdots, a_{m}:

[C^0C^1C^m1C^1C^0C^m2C^m1C^m2C^0][a1a2am]=[C^1C^2C^m]\begin{aligned} \begin{bmatrix} \hat{C}_{0} & \hat{C}_{1} & \cdots \hat{C}_{m-1}\\ \hat{C}_{1} & \hat{C}_{0} & \cdots \hat{C}_{m-2}\\ \vdots & \vdots & \ddots & \vdots\\ \hat{C}_{m-1} & \hat{C}_{m-2} & \cdots \hat{C}_{0} \end{bmatrix} \begin{bmatrix} a_{1}\\ a_{2}\\ \vdots\\ a_{m} \end{bmatrix} &= \begin{bmatrix} \hat{C}_{1}\\ \hat{C}_{2}\\ \vdots\\ \hat{C}_{m} \end{bmatrix} \end{aligned}

The estimates a^i\hat{a}_{i} of the AR coefficients are obtained by solving this equation. Then from the Yule-Walker equations, an estimate of the variance is obtained by:

σ^2=C^0i=1ma^iC^i\begin{aligned} \hat{\sigma}^{2} &= \hat{C}_{0} - \sum_{i=1}^{m}\hat{a}_{i}\hat{C}_{i} \end{aligned}

We can also obtain the following Yule-Walker equation from minimizing the prediction errors:

Cj=i=1maiCji\begin{aligned} C_{j} &= \sum_{i=1}^{m}a_{i}C_{j-i} \end{aligned} E[vn2]=E[(yni=1maiyni)2]=C02i=1maiCi+i=1mj=1maiajCij\begin{aligned} E[v_{n}^{2}] &= E[(y_{n}- \sum_{i=1}^{m}a_{i}y_{n-i})^{2}]\\ &= C_{0} - 2\sum_{i=1}^{m}a_{i}C_{i} + \sum_{i=1}^{m}\sum_{j=1}^{m}a_{i}a_{j}C_{i-j} \end{aligned} E[vn2]aj=2Cj+2i=1maiCji=0\begin{aligned} \frac{\partial E[v_{n}^{2}]}{\partial a_{j}} &= -2C_{j} + 2\sum_{i=1}^{m}a_{i}C_{j-i} &= 0 \end{aligned}

Levinson’s Algorithm

To obtain the Yule-Walker estimates for an AR model of order m, it is necessary to solve a system of linear equations with m unknowns. In addition, to select the order of the AR model by the minimum AIC method, we need to evaluate the AIC values of the models with orders up to the maximum order, M. Namely, we have to estimate the coefficients by solving systems of linear equations with 1 unknown, …, M unknowns.

However, by the Levinson’s algorithm, these solutions can be obtained quite efficiently. Hereafter, the AR coefficients and the innovation variance of the AR model of order m are denoted as ajma_{j}^{m} and σm2\sigma_{m}^{2},respectively.

Then Levinson’s algorithm is defined as follows:

Set:

σ^02=C^0\begin{aligned} \hat{\sigma}_{0}^{2} &= \hat{C}_{0} \end{aligned} AIC0=N(log2πσ^02+1)+2\begin{aligned} \text{AIC}_{0} &= N(\log 2\pi\hat{\sigma}_{0}^{2} + 1) + 2 \end{aligned}

For m=1,,Mm = 1, \cdots, M

a^mm=C^mj=1m1C^mjσ^m1\begin{aligned} \hat{a}_{m}^{m} &= \frac{\hat{C}_{m} - \sum_{j=1}^{m-1}\hat{C}_{m-j}}{\hat{\sigma}_{m-1}} \end{aligned}

Recall a^mm\hat{a}_{m}^{m} is the PARCOR.

Calculate aia_{i} from i=1,,m1i = 1, \cdots, m-1:

a^im=a^im1a^mma^mim1\begin{aligned} \hat{a}_{i}^{m} &= \hat{a}_{i}^{m-1} - \hat{a}_{m}^{m}\hat{a}_{m-i}^{m-1} \end{aligned} σ^m2=σ^m12[1(a^mm)2]\begin{aligned} \hat{\sigma}_{m}^{2} &= \hat{\sigma}_{m-1}^{2}[1 - (\hat{a}_{m}^{m})^{2}] \end{aligned} AICm=N(log2πσ^m2+1)+2(m+1)\begin{aligned} \text{AIC}_{m} &= N(\log 2\pi \hat{\sigma}_{m}^{2} + 1) + 2(m+1) \end{aligned}

Derivation of Levinson’s Algorithm

For multivariate time series, the coefficients of the backward AR model are different from those of the forward AR model, the algorithm needs to treat twice as many number of recursions as the univariate case.

Assume that the autocovariance function CkC_{k}, k=0,1,k = 0, 1, \cdots is given and that the Yule-Walker estimates of the AR model of order m − 1 have already been obtained.

In the following, since we consider AR models with different orders, the AR coefficients, the prediction error and the prediction error variance of the AR model of order m − 1 are denoted by a^1m1,,a^m1m1,vnm1,σ^m12\hat{a}_{1}^{m-1}, \cdots, \hat{a}_{m-1}^{m-1},v_{n}^{m-1},\hat{\sigma}_{m-1}^{2}.

Since the coefficients a^1m1,,a^m1m1\hat{a}_{1}^{m-1}, \cdots, \hat{a}_{m-1}^{m-1} satisfy the following Yule-Walker equation of order m-1, we have:

Ck=j=1m1a^jm1Cjk\begin{aligned} C_{k} &= \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}C_{j-k} \end{aligned} k=1,,m1\begin{aligned} k &= 1, \cdots, m-1 \end{aligned}

This means that vnm1v_{n}^{m-1} satisfies:

E[vnm1ynk]=E[(ynj=1m1a^jm1ynj)ynk]=0\begin{aligned} E[v_{n}^{m-1}y_{n-k}] &= E[(y_{n} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}y_{n-j})y_{n-k}] = 0 \end{aligned} k=1,,m1\begin{aligned} k &= 1, \cdots, m-1 \end{aligned}

Now we consider a backward AR model of order m − 1 that represents the present value of the time series yn with future values yn+1,,yn+m1y_{n+1}, \cdots, y_{n+m−1}:

yn=i=1m1ajyn+i+wnm1\begin{aligned} y_{n} &= \sum_{i=1}^{m-1}a_{j}y_{n+i} + w_{n}^{m-1} \end{aligned}

Then, since Ck=CkC_{−k} = C_{k} holds for a univariate time series, this backward AR model also satisfies the same Yule-Walker equation. Therefore, the backward prediction error wnm1w_{n}^{m−1} satisfies:

ynm=i=1m1ajynm+i+wnmm1\begin{aligned} y_{n-m} &= \sum_{i=1}^{m-1}a_{j}y_{n-m+i} + w_{n-m}^{m-1} \end{aligned} E[wnmm1ynk]=E[(ynmj=1m1a^jm1ynm+j)ynk]=0\begin{aligned} E[w_{n-m}^{m-1}y_{n-k}] &= E[(y_{n-m} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}y_{n-m +j})y_{n-k}] = 0 \end{aligned} k=1,,m1\begin{aligned} k &= 1, \cdots, m-1 \end{aligned}

Now, we define the weighted sum of the forward prediction error vnm1v_{n}^{m−1} and the backward prediction error wnmm1w_{n-m}^{m−1} by:

zn=vnm1βwnmm1=ynj=1m1a^jm1ynjβ(ynmj=1m1a^jm1ynm+j)=ynj=1m1[a^jm1βa^mjm1]ynjβynm\begin{aligned} z_{n} &= v_{n}^{m-1} - \beta w_{n-m}^{m-1}\\ &= y_{n} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}y_{n-j} - \beta\Big(y_{n-m} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}y_{n-m+j}\Big)\\ &= y_{n} - \sum_{j=1}^{m-1}[\hat{a}_{j}^{m-1} - \beta \hat{a}_{m-j}^{m-1}]y_{n-j} - \beta y_{n-m} \end{aligned}

znz_{n} is a linear combination of yn,yn1,,ynmy_{n}, y_{n−1}, \cdots, y_{n−m} and from:

E[vnm1ynk]=0E[wnmm1ynk]=0\begin{aligned} E[v_{n}^{m-1}y_{n-k}] &= 0\\ E[w_{n-m}^{m-1}y_{n-k}] &= 0 \end{aligned}

It follows that:

E[znynk]=0\begin{aligned} E[z_{n}y_{n-k}] &= 0 \end{aligned} k=1,,m1\begin{aligned} k = 1, \cdots, m-1 \end{aligned}

Therefore, if the value of β\beta is determined so that E(znynm)=0E(z_{n}y_{n−m}) = 0 holds for k = m, the Yule-Walker equation of order m will be satisfied:

E[znynm]=E[ynj=1m1a^jm1ynjβ(ynmj=1m1ynm+j)]ynm=Cmj=1m1a^jm1Cmjβ(C0j=1m1a^jm1Cj)=0\begin{aligned} E[z_{n}y_{n-m}] &= E\Big[y_{n} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}y_{n-j} - \beta\Big(y_{n-m} - \sum_{j=1}^{m-1}y_{n-m+j}\Big)\Big]y_{n-m}\\ &= C_{m} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}C_{m-j} - \beta(C_{0} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}C_{j})\\ &= 0 \end{aligned}

β\beta is then determined as:

β=Cmj=1m1a^jm1CmjC0j=1m1a^jm1Cmj=Cmj=1m1a^jm1Cmjσm12\begin{aligned} \beta &= \frac{C_{m} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}C_{m-j}}{C_{0} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}C_{m-j}}\\[5pt] &= \frac{C_{m} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}C_{m-j}}{\sigma_{m-1}^{2}} \end{aligned}

Recall:

zn=ynj=1m1[a^jm1βa^m1m1]ynjβynm\begin{aligned} z_{n} &= y_{n} - \sum_{j=1}^{m-1}[\hat{a}_{j}^{m-1} - \beta \hat{a}_{m-1}^{m-1}]y_{n-j} - \beta y_{n-m} \end{aligned}

If we define β\beta as a^mm\hat{a}_{m}^{m} and a^jm\hat{a}_{j}^{m} as:

βa^mm\begin{aligned} \beta &\equiv \hat{a}_{m}^{m} \end{aligned} a^jma^jm1a^mma^mjm1\begin{aligned} \hat{a}_{j}^{m} &\equiv \hat{a}_{j}^{m-1} - \hat{a}_{m}^{m}\hat{a}_{m-j}^{m-1} \end{aligned}

Then:

E[znynk]=E[(ynj=1ma^jmynj)ynk]=E[ynynk]j=1ma^jmynjynk=Ckj=1ma^jmCkj=0\begin{aligned} E[z_{n}y_{n-k}] &= E\Big[(y_{n} - \sum_{j=1}^{m}\hat{a}_{j}^{m}y_{n-j})y_{n-k}\Big]\\ &= E[y_{n}y_{n-k}] - \sum_{j=1}^{m}\hat{a}_{j}^{m}y_{n-j}y_{n-k}\\ &= C_{k} - \sum_{j=1}^{m}\hat{a}_{j}^{m}C_{k-j}\\ &= 0 \end{aligned} zn=ynj=1m1[a^jm1βa^m1m1]ynjβynm=ynj=1ma^jmynj\begin{aligned} z_{n} &= y_{n} - \sum_{j=1}^{m-1}[\hat{a}_{j}^{m-1} - \beta \hat{a}_{m-1}^{m-1}]y_{n-j} - \beta y_{n-m}\\ &= y_{n} - \sum_{j=1}^{m}\hat{a}_{j}^{m}y_{n-j} \end{aligned}

Hence, we get the Yule-Walker equation of order m:

Ck=j=1ma^jmCkj\begin{aligned} C_{k} &= \sum_{j=1}^{m}\hat{a}_{j}^{m}C_{k-j}\\ \end{aligned} k=1,,m\begin{aligned} k = 1, \cdots, m \end{aligned}

From the Yule-Walker equation, we can express the innovation variance of the AR model as:

σ^m2=C0j=1ma^jmCj=C0j=1m(a^jm1a^mma^mjm1)Cja^mmCm=C0j=1ma^jm1Cja^mm(Cmj=1m1a^mjm1Cj)=σ^m12a^mm(Cmj=1m1a^m1m1Cj)\begin{aligned} \hat{\sigma}_{m}^{2} &= C_{0} - \sum_{j=1}^{m}\hat{a}_{j}^{m}C_{j}\\ &= C_{0} - \sum_{j=1}^{m}(\hat{a}_{j}^{m-1} - \hat{a}_{m}^{m}\hat{a}_{m-j}^{m-1})C_{j} - \hat{a}_{m}^{m}C_{m}\\ &= C_{0} - \sum_{j=1}^{m}\hat{a}_{j}^{m-1}C_{j} - \hat{a}_{m}^{m}\Bigg(C_{m} - \sum_{j=1}^{m-1}\hat{a}_{m-j}^{m-1}C_{j}\Bigg)\\ &= \hat{\sigma}_{m-1}^{2} - \hat{a}_{m}^{m}\Bigg(C_{m} - \sum_{j=1}^{m-1}\hat{a}_{m-1}^{m-1}C_{j}\Bigg) \end{aligned}

Recall that:

β=a^mm=Cmj=1m1a^m1m1Cjσ^m12\begin{aligned} \beta &= \hat{a}_{m}^{m}\\ &= \frac{C_{m} - \sum_{j=1}^{m-1}\hat{a}_{m-1}^{m-1}C_{j}}{\hat{\sigma}_{m-1}^{2}} \end{aligned} Cmj=1m1a^m1m1Cj=a^mmσ^m12\begin{aligned} C_{m} - \sum_{j=1}^{m-1}\hat{a}_{m-1}^{m-1}C_{j} &= \hat{a}_{m}^{m}\hat{\sigma}_{m-1}^{2} \end{aligned}

Hence:

σ^m2=σ^m12a^mm(Cmj=1m1a^m1m1Cj)=σ^m12a^mm(a^mmσ^m12)=σ^m12[1(a^mm)2]\begin{aligned} \hat{\sigma}_{m}^{2} &= \hat{\sigma}_{m-1}^{2} - \hat{a}_{m}^{m}\Big(C_{m} - \sum_{j=1}^{m-1}\hat{a}_{m-1}^{m-1}C_{j}\Big)\\ &= \hat{\sigma}_{m-1}^{2} - \hat{a}_{m}^{m}(\hat{a}_{m}^{m}\hat{\sigma}_{m-1}^{2})\\[5pt] &= \hat{\sigma}_{m-1}^{2}[1 - (\hat{a}_{m}^{m})^{2}] \end{aligned}

To show derivation for the AIC, we first note that:

p(ynyn1,,ynm)=12log2πσm212σm2(yni=1maiyni)2\begin{aligned} p(y_{n}|y_{n-1}, \cdots, y_{n-m}) &= -\frac{1}{2}\log 2\pi \sigma_{m}^{2} - \frac{1}{2\sigma_{m}^{2}}\Bigg(y_{n} - \sum_{i=1}^{m}a_{i}y_{n-i}\Bigg)^{2} \end{aligned}

Hence:

L(θ)=p(y1,θ)n=2Np(yny1,,yn1,θ)l(θ)=p(y1,θ)+n=2Np(yny1,,yn1,θ)=N2log2πσm2N2σm2σm2N=N2log2πσm212\begin{aligned} L(\boldsymbol{\theta}) &= -p(y_{1}, \boldsymbol{\theta})\prod_{n=2}^{N}p(y_{n}|y_{1}, \cdots, y_{n-1}, \boldsymbol{\theta})\\ l(\boldsymbol{\theta}) &= -p(y_{1}, \boldsymbol{\theta}) + \sum_{n=2}^{N}-p(y_{n}|y_{1}, \cdots, y_{n-1}, \boldsymbol{\theta})\\ &= -\frac{N}{2}\log 2\pi \sigma_{m}^{2} - \frac{N}{2\sigma_{m}^{2}}\frac{\sigma_{m}^{2}}{N}\\ &= -\frac{N}{2}\log 2\pi \sigma_{m}^{2} - \frac{1}{2} \end{aligned} AICm=2l(θ)+2(m+1)=N(log2πσ^m2+1)+2(m+1)\begin{aligned} \text{AIC}_{m} &= -2l(\boldsymbol{\theta}) + 2(m+1)\\ &= N(\log 2 \pi \hat{\sigma}_{m}^{2} + 1) + 2(m+1) \end{aligned}

Estimation of an AR Model by the Least Squares Method

The log-likelihood of the AR model:

l(θ)=n=1Nlogp(yny1,,yn1)\begin{aligned} l(\boldsymbol{\theta}) &= \sum_{n=1}^{N}\log -p(y_{n}|y_{1}, \cdots, y_{n-1}) \end{aligned} θ=[a1amσ2]\begin{aligned} \boldsymbol{\theta} &= \begin{bmatrix} a_{1}\\ \vdots\\ a_{m}\\ \sigma^{2} \end{bmatrix} \end{aligned}

Here, for the AR model of order m, since the distribution of yny_{n} is specified by the values of yn1,,ynmy_{n−1}, \cdots, y_{n−m} for n>mn > m:

p(yny1,,yn1)=p(ynynm,,yn1)=12πσ2e12σ2(yni=1maiyni)2\begin{aligned} -p(y_{n}|y_{1}, \cdots, y_{n-1}) &= -p(y_{n}|y_{n-m}, \cdots, y_{n-1})\\ &= -\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{1}{2\sigma^{2}}(y_{n} - \sum_{i=1}^{m}a_{i}y_{n-i})^{2}} \end{aligned}

And taking the log:

p(ynynm,,yn1)=12log2πσ212σ2(yni=1maiyni)2\begin{aligned} - p(y_{n}|y_{n-m}, \cdots, y_{n-1}) &= -\frac{1}{2}\log 2\pi\sigma^{2} - \frac{1}{2\sigma^{2}}(y_{n}- \sum_{i=1}^{m}a_{i}y_{n-i})^{2} \end{aligned}

Therefore, by ignoring the initial MM (MmM \geq m) terms, an approximate log-likelihood of the AR model is obtained as:

l(θ)=n=M+1Nlogp(ynynm,,yn1)=NM2log2πσ212σ2n=M+1N(yni=1maiyni)2\begin{aligned} l(\boldsymbol{\theta}) &= \sum_{n = M + 1}^{N}-\log p(y_{n}|y_{n-m}, \cdots, y_{n-1})\\ &= -\frac{N-M}{2}\log 2\pi\sigma^{2} - \frac{1}{2\sigma^{2}}\sum_{n=M+1}^{N}\Bigg(y_{n} - \sum_{i=1}^{m}a_{i}y_{n-i}\Bigg)^{2} \end{aligned}

Similar to the case of the regression model, for arbitrarily given autoregressive coefficients a1,,ama_{1}, \cdots, a_{m}, an approximate maximum likelihood estimate of the variance σ2\sigma^{2} maximizing satisfies:

l(θ)σ2=NM2σ2+12(σ2)2n=M+1N(yni=1maiyni)2=0\begin{aligned} \frac{\partial l(\boldsymbol{\theta})}{\partial \sigma^{2}} &= -\frac{N-M}{2\sigma^{2}} + \frac{1}{2(\sigma^{2})^{2}}\sum_{n=M+1}^{N}\Bigg(y_{n} - \sum_{i=1}^{m}a_{i}y_{n-i}\Bigg)^{2} = 0 \end{aligned} σ^2=1NMn=M+1N(yni=1maiyni)2\begin{aligned} \hat{\sigma}^{2} &= \frac{1}{N-M}\sum_{n=M+1}^{N}\Bigg(y_{n} - \sum_{i=1}^{m}a_{i}y_{n-i}\Bigg)^{2} \end{aligned}

Subsituting this back into the log-likelihood, the the approximate log-likelihood becomes a function of the autoregressive coefficients a1,,ama_{1}, \cdots, a_{m}:

l(a1,,am)=NM2log2πσ^2NM2\begin{aligned} l(a_{1}, \cdots, a_{m}) &= -\frac{N-M}{2}\log 2\pi\hat{\sigma}^{2} - \frac{N-M}{2} \end{aligned}

Here, since the logarithm is a monotone increasing function, maximization of the approximate log-likelihood can be achieved by minimizing the variance σ^2\hat{\sigma}^{2}. This means that the approximate maximum likelihood estimates of the AR model can be obtained by the least squares method.

To obtain the least squares estimates of the AR models with orders up to M by the Householder transformation, define the matrix Z\textbf{Z} and the vector y\textbf{y} by:

Z=[yMyM1y1yM+1yMy2yN1yN2yNM]\begin{aligned} \mathbf{Z} &= \begin{bmatrix} y_{M} & y_{M-1} & \cdots & y_{1}\\ y_{M + 1} & y_{M} & \cdots & y_{2}\\ \vdots & \vdots & \ddots & \vdots\\ y_{N-1} & y_{N-2} & \cdots & y_{N-M} \end{bmatrix} \end{aligned} y=[yM+1yM+2yN]\begin{aligned} \textbf{y} &= \begin{bmatrix} y_{M+1}\\ y_{M+2}\\ \vdots\\ y_{N} \end{bmatrix} \end{aligned}

Construct the (NM)×(M+1)(N - M)\times (M+1) matrix:

X=[Zy]=[yMyM1y1yM+1yM+1yMy2yM+2yN1yN2yNMyN]\begin{aligned} \textbf{X} &= [\textbf{Z}|\textbf{y}]\\ &= \begin{bmatrix} y_{M} & y_{M-1} & \cdots & y_{1} & y_{M+1}\\ y_{M + 1} & y_{M} & \cdots & y_{2} & y_{M+2}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ y_{N-1} & y_{N-2} & \cdots & y_{N-M} & y_{N} \end{bmatrix} \end{aligned}

And transform it to an upper triangular matrix by Householder transformation:

HX=[S0]=[s11s1Ms1,M+1sMMsM,M+1sM+1,M+10000]\begin{aligned} \textbf{HX} &= \begin{bmatrix} \textbf{S}\\ \textbf{0} \end{bmatrix}\\ &= \begin{bmatrix} s_{11} & \cdots & s_{1M} & s_{1, M+1}\\ & \ddots & \vdots & \vdots\\ & & s_{MM} & s_{M, M+1}\\ & & & s_{M+1, M+1}\\ 0 & 0 & 0 & 0 \end{bmatrix} \end{aligned}

Then, for 0jM0 \leq j \leq M, the innovation variance and the AIC of the AR model of order j are obtained by:

σ^j2=1NMi=j+1M+1si,M+12AICj=(NM)(log2πσ^j2+1)+2(j+1)\begin{aligned} \hat{\sigma}_{j}^{2} &= \frac{1}{N-M}\sum_{i=j+1}^{M+1}s_{i,M+1}^{2}\\[5pt] \text{AIC}_{j} &= (N - M)(\log 2\pi\hat{\sigma}_{j}^{2} + 1) + 2(j+1) \end{aligned}

Moreover, the least squares estimates of the autoregressive coefficients, which are the solutions of the linear equation:

[s11s1jsjj][a1aj]=[s1,M+1sj,M+1]\begin{aligned} \begin{bmatrix} s_{11} & \cdots & s_{1j}\\ & \ddots & \vdots\\ & & s_{jj} \end{bmatrix} \begin{bmatrix} a_{1}\\ \vdots\\ a_{j} \end{bmatrix} &= \begin{bmatrix} s_{1, M+1}\\ \vdots\\ s_{j, M+1} \end{bmatrix} \end{aligned}

Can be easily obtained by backward substitution as follows:

a^j=sj,M+1sjja^i=si,M+1si,i+1a^i+1si,jaj^sjj\begin{aligned} \hat{a}_{j} &= \frac{s_{j, M+1}}{s_{jj}}\\ \hat{a}_{i} &= \frac{s_{i,M+1}- s_{i, i+1}\hat{a}_{i+1} - \cdots - s_{i,j}\hat{a_{j}}}{s_{jj}} \end{aligned} i=j1,,1\begin{aligned} i = j-1, \cdots, 1 \end{aligned}

Estimation of an AR Model by the PARCOR Method

Assuming that the autocovariance functions C0,C1,C_{0},C_{1}, \cdots are given, Levinson’s algorithm can be executed by using the following relation between the coefficients of the AR model of order m − 1 and the coefficients of the AR model of order m:

ajm=ajm1ammamjm1\begin{aligned} a_{j}^{m} &= a_{j}^{m-1} - a_{m}^{m}a_{m-j}^{m-1} \end{aligned}

In the derivation of Levinson’s Algorithm section, we shown that:

β=a^mm=Cmj=1m1a^m1m1Cjσ^m12\begin{aligned} \beta &= \hat{a}_{m}^{m}\\ &= \frac{C_{m} - \sum_{j=1}^{m-1}\hat{a}_{m-1}^{m-1}C_{j}}{\hat{\sigma}_{m-1}^{2}} \end{aligned}

Therefore, if we can estimate the amma_{m}^{m}, the other coefficients can be automatically determined by using the above relation. In Levinson’s algorithm, we estimated amma_{m}^{m} by using the following formula that is obtained by substituting the sample autocovariance functions C^0,,C^m\hat{C}_{0}, \cdots, \hat{C}_{m}:

a^mm=C^mj=1m1a^jm1C^mjC^0j=1m1a^jm1C^j=C^mj=1m1a^jm1C^mjσ^m12\begin{aligned} \hat{a}_{m}^{m} &= \frac{\hat{C}_{m} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}\hat{C}_{m-j}}{\hat{C}_{0} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}\hat{C}_{j}}\\ &= \frac{\hat{C}_{m} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}\hat{C}_{m-j}}{\hat{\sigma}_{m-1}^{2}} \end{aligned}

In this section, we present another method of estimating PARCOR amma_{m}^{m} directly from the time series y1,,yNy_{1}, \cdots, y_{N} without using the sample autocovariance functions.

First of all, let wnm1w_{n}^{m-1} denote the prediction error of the backward AR model with order m − 1:

yn=j=1m1ajm1yn+j+wnm1\begin{aligned} y_{n} &= \sum_{j=1}^{m-1}a_{j}^{m-1}y_{n+j} + w_{n}^{m-1} \end{aligned}

Then:

ynm=j=1m1ajm1ynm+j+wnmm1\begin{aligned} y_{n -m} &= \sum_{j=1}^{m-1}a_{j}^{m-1}y_{n-m+j} + w_{n-m}^{m-1} \end{aligned}

In the case of a univariate time series, since the autocovariance function is an even function, the AR coefficients of the forward model coincide with those of the backward model.

From the Yule-Walker equation:

Cmj=1m1ajm1Cmj=E[(ynj=1m1ajm1ynj)ynm]=E[vnm1ynm]=E[vnm1(j=1m1ajm1ynm+j+wnmm1)]=j=1m1ajm1E[vnm1ynm+j]+E[vnm1wnmm1]=0+E[vnm1wnmm1]1Nmn=m+1Nvnm1wnmm1\begin{aligned} C_{m} - \sum_{j=1}^{m-1}a_{j}^{m-1}C_{m-j} &= E\Big[(y_{n} - \sum_{j=1}^{m-1}a_{j}^{m-1}y_{n-j})y_{n-m}\Big]\\ &= E[v_{n}^{m-1}y_{n-m}]\\ &= E[v_{n}^{m-1}\Bigg(\sum_{j=1}^{m-1}a_{j}^{m-1}y_{n-m+j} + w_{n-m}^{m-1}\Bigg)]\\ &= \sum_{j=1}^{m-1}a_{j}^{m-1}E[v_{n}^{m-1}y_{n-m+j}] + E[v_{n}^{m-1}w_{n-m}^{m-1}]\\ &= 0 + E[v_{n}^{m-1}w_{n-m}^{m-1}]\\ &\approx \frac{1}{N-m}\sum_{n=m+1}^{N}v_{n}^{m-1}w_{n-m}^{m-1} \end{aligned}

On the other hand we can express:

C0j=1m1ajm1Cj=E[ynmynmj=1m1ajm1ynm+jynm]=E[(ynmj=1m1ajm1ynm+j)ynm]=E[wnmm1ynm]=E[wnmm1j=1m1ajm1ynm+j+wnmm1]=j=1m1ajm1E[wnmm1ynm+j]+E[wnmm1wnmm1]=0+E[wnmm1]2=E[wnmm1]2\begin{aligned} C_{0} - \sum_{j=1}^{m-1}a_{j}^{m-1}C_{j} &= E[y_{n-m}y_{n-m} - \sum_{j=1}^{m-1}a_{j}^{m-1}y_{n-m+j}y_{n-m}]\\ &= E[(y_{n-m} - \sum_{j=1}^{m-1}a_{j}^{m-1}y_{n-m+j})y_{n-m}]\\ &= E[w_{n-m}^{m-1}y_{n-m}]\\ &= E\Big[w_{n-m}^{m-1}\sum_{j=1}^{m-1}a_{j}^{m-1}y_{n-m+j} + w_{n-m}^{m-1}\Big]\\ &= \sum_{j=1}^{m-1}a_{j}^{m-1}E[w_{n-m}^{m-1}y_{n-m+j}] + E[w_{n-m}^{m-1}w_{n-m}^{m-1}]\\ &= 0 + E[w_{n-m}^{m-1}]^{2}\\ &= E[w_{n-m}^{m-1}]^{2} \end{aligned}

Since from Yule-Walker equation:

C0j=1m1ajm1Cj=E[vnm1]2=E[wnmm1]2=σ2\begin{aligned} C_{0} - \sum_{j=1}^{m-1}a_{j}^{m-1}C_{j} &= E[v_{n}^{m-1}]^{2}\\ &= E[w_{n-m}^{m-1}]^{2}\\ &= \sigma^{2} \end{aligned}

Various estimates of E[vnm1wnmm1]E[v_{n}^{m-1}w_{n-m}^{m-1}] can be obtained:

1Nmn=m+1N(wnmm1)2\begin{aligned} \frac{1}{N-m}\sum_{n=m+1}^{N}(w_{n-m}^{m-1})^{2} \end{aligned} 1Nm[n=m+1N(wnmm1)2n=m+1N(vnm1)2]12\begin{aligned} \frac{1}{N-m}\Big[\sum_{n=m+1}^{N}(w_{n-m}^{m-1})^{2} \sum_{n=m+1}^{N}(v_{n}^{m-1})^{2} \Big]^{\frac{1}{2}} \end{aligned} 12(Nm)[n=m+1N(wnmm1)2+n=m+1N(vnm1)2]\begin{aligned} \frac{1}{2(N-m)}\Big[\sum_{n=m+1}^{N}(w_{n-m}^{m-1})^{2} + \sum_{n=m+1}^{N}(v_{n}^{m-1})^{2}\Big] \end{aligned}

Recall:

a^mm=C^mj=1m1a^jm1C^mjC^0j=1m1a^jm1C^j\begin{aligned} \hat{a}_{m}^{m} &= \frac{\hat{C}_{m} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}\hat{C}_{m-j}}{\hat{C}_{0} - \sum_{j=1}^{m-1}\hat{a}_{j}^{m-1}\hat{C}_{j}} \end{aligned} Cmj=1m1ajm1Cmj1Nmn=m+1Nvnm1wnmm1\begin{aligned} C_{m} - \sum_{j=1}^{m-1}a_{j}^{m-1}C_{m-j} &\approx \frac{1}{N-m}\sum_{n=m+1}^{N}v_{n}^{m-1}w_{n-m}^{m-1} \end{aligned} C0j=1m1ajm1Cjn=m+1Nwnmm1]2Nm\begin{aligned} C_{0} - \sum_{j=1}^{m-1}a_{j}^{m-1}C_{j} &\approx \frac{\sum_{n = m+1}^{N}w_{n-m}^{m-1}]^{2}}{N-m} \end{aligned}

Hence we can approximate the PARCOR in 3 ways.

The following estimate is the regression coefficient when the prediction error vnm1v_{n}^{m-1} of the forward model is regressed on the prediction error wnmm1w_{n-m}^{m-1} of the backward model.

a^mm=n=m+1Nvnm1wnmm1n=m+1N(wnmm1)2\begin{aligned} \hat{a}_{m}^{m} &= \frac{\sum_{n=m+1}^{N}v_{n}^{m-1}w_{n-m}^{m-1}}{\sum_{n=m+1}^{N}(w_{n-m}^{m-1})^{2}} \end{aligned}

The following estimate corresponds to the definition of PARCOR, since it is the correlation coefficient of vnm1v_{n}^{m-1} and wnmm1w_{n-m}^{m-1}:

a^mm=n=m+1Nvnm1wnmm1n=m+1N(wnmm1)2n=m+1N(vnmm1)2\begin{aligned} \hat{a}_{m}^{m} &= \frac{\sum_{n=m+1}^{N}v_{n}^{m-1}w_{n-m}^{m-1}}{\sqrt{\sum_{n=m+1}^{N}(w_{n-m}^{m-1})^{2}\sum_{n=m+1}^{N}(v_{n-m}^{m-1})^{2}}} \end{aligned}

The following estimate minimizes the mean of the variances of the forward prediction errors and the backward prediction errors, and consequently Burg’s algorithm based on the maximum entropy method (MEM):

a^mm=2n=m+1Nvnm1wnmm1n=m+1N(wnmm1)2+n=m+1N(vnmm1)2\begin{aligned} \hat{a}_{m}^{m} &= 2\frac{\sum_{n=m+1}^{N}v_{n}^{m-1}w_{n-m}^{m-1}}{\sum_{n=m+1}^{N}(w_{n-m}^{m-1})^{2} + \sum_{n=m+1}^{N}(v_{n-m}^{m-1})^{2}} \end{aligned}

A generic procedure to estimate the AR model from the time series y1,,yNy_{1},\cdots, y_{N} using the PARCOR method are described below. Here, for simplicity, the mean value of the time series yny_{n} is assumed to be 0:

1) Set:

vn0=wn0=yn\begin{aligned} v_{n}^{0} &= w_{n}^{0} = y_{n} \end{aligned} n=1,,N\begin{aligned} n = 1, \cdots, N \end{aligned}

And compute:

σ^02=n=1Nyn2N\begin{aligned} \hat{\sigma}_{0}^{2} &= \frac{\sum_{n=1}^{N}y_{n}^{2}}{N} \end{aligned} AIC0=N(log2πσ^0+1)+2\begin{aligned} \text{AIC}_{0} &= N(\log 2 \pi \hat{\sigma}_{0} + 1) + 2 \end{aligned}

2) For m=1,,Mm = 1, \cdots, M, repeat the following steps:

(a) Estimate a^mm\hat{a}_{m}^{m} by any of the above 3 formulas.

See the “Estimation of an AR Model by the PARCOR Method” section for above derivation.

(b) Obtain AR coefficients a^1m,,a^m1m\hat{a}_{1}^{m}, \cdots, \hat{a}_{m-1}^{m}:

ajm=ajm1ammamjm1\begin{aligned} a_{j}^{m} &= a_{j}^{m-1} - a_{m}^{m}a_{m-j}^{m-1} \end{aligned}

(c) For n=m+1,,Nn = m + 1, \cdots, N, obtain the forward prediction error by:

vnm=vnm1a^mmwnmm1\begin{aligned} v_{n}^{m} &= v_{n}^{m-1} - \hat{a}_{m}^{m}w_{n-m}^{m-1} \end{aligned}

(d) For n=m+1,,Nn = m + 1, \cdots, N, obtain the backward prediction error by:

wnmm=wnmm1a^mmvnm1\begin{aligned} w_{n-m}^{m} &= w_{n-m}^{m-1} - \hat{a}_{m}^{m}v_{n}^{m-1} \end{aligned}

(e) Estimate the innovation variance of the AR model of order m by:

σ^m2=σ^m12[1(a^mm)2]\begin{aligned} \hat{\sigma}_{m}^{2} &= \hat{\sigma}_{m-1}^{2}[1 - (\hat{a}_{m}^{m})^{2}] \end{aligned}

(f) Obtain AIC by:

AICm=N(log2πσ^m2+1)+2(m+1)\begin{aligned} \text{AIC}_{m} &= N(\log 2\pi\hat{\sigma}_{m}^{2} + 1) + 2(m+1) \end{aligned}

Example (AR modeling for sunspot number data)

The function arfit of the package TSSS fits AR model based on the AIC criterion. The method denotes the type of estimation method:

  • Yule-Walker method (Default)
  • Least squares (Householder) method
  • PARCOR method (partial autoregression)
  • PARCOR method (PARCOR)
  • PARCOR method (Burg’s algorithm)

And the outputs are:

  • sigma2: innovation variances
  • maice.order: minimum AIC order
  • aic: AIC of the AR models
  • arcoef: AR coefficients of AR models
  • parcor: PARCORs
  • spec: power spectrum (in log scale) of the AIC best AR model
fit <- arfit(
  log10(Sunspot$value),
  lag = 20,
  method = 3
)

> names(fit)
[1] "sigma2"  "maice.order"  "aic"  "arcoef"  "parcor"     
[6] "spec"    "tsname"

Estimation of Multivariate AR Model by the Yule-Walker Method

kk denotes the number of variables (or dimensions) of a multivariate time series. The parameters of the multivariate AR model with order m:

yn=i=1mAimyni+vn\begin{aligned} \textbf{y}_{n} &= \sum_{i=1}^{m}\textbf{A}_{i}^{m}\textbf{y}_{n-i} + \textbf{v}_{n} \end{aligned} vnN(0,Vm)\begin{aligned} \textbf{v}_{n} &\sim N(\textbf{0}, \textbf{V}_{m}) \end{aligned}

Vm\textbf{V}_{m} is the covariance matrix of the innovation vn\textbf{v}_{n} and the AR coefficient matrices A1m,,Amm\textbf{A}_{1}^{m}, \cdots, \textbf{A}_{m}^{m}.

When a multivariate AR model is given, the cross-covariance function is obtained from:

C0=j=1MAjCj+WCk=j=1MAjCkj\begin{aligned} \textbf{C}_{0} &= \sum_{j=1}^{M}\textbf{A}_{j}\textbf{C}_{-j} + \textbf{W}\\ \textbf{C}_{k} &= \sum_{j=1}^{M}\textbf{A}_{j}\textbf{C}_{k-j} \end{aligned} W=E[vnvnT]=[σ11σ1lσl1σll]\begin{aligned} \textbf{W} &= E[\textbf{v}_{n}\textbf{v}_{n}^{T}]\\ &= \begin{bmatrix} \sigma_{11} & \cdots & \sigma_{1l}\\ \vdots & \ddots &\vdots\\ \sigma_{l1} & \cdots & \sigma_{ll}\\ \end{bmatrix} \end{aligned} k=1,2,\begin{aligned} k = 1, 2, \cdots \end{aligned}

On the other hand, using these equations, the estimates of the parameters of the multivariate AR model can be obtained through the sample cross-covariance function.

For actual computation, similarly to the univariate AR model, there is a computationally efficient algorithm. However, for a multivariate time series, the backward model is different from the forward model. In the case of a univariate time series, the forward AR model coincides with the backward AR model, because the autocovariance function is an even function. But this property is not satisfied by multivariate time series. Therefore, in order to derive an efficient algorithm similar to Levinson’s algorithm, we have to consider the backward multivariate AR model:

yn=i=1mBimyn+i+un\begin{aligned} \textbf{y}_{n} &= \sum_{i=1}^{m}\textbf{B}_{i}^{m}y_{n+i} + \textbf{u}_{n}\\ \end{aligned}

And we need to estimate the variance covariance matrix Um\textbf{U}_{m} and the coefficients Bim\textbf{B}_{i}^{m}, as well as Vm\textbf{V}_{m} and Aim\textbf{A}_{i}^{m} simultaneously.

yn=i=1mBimyn+i+un\begin{aligned} \textbf{y}_{n} &= \sum_{i=1}^{m}\textbf{B}_{i}^{m}\textbf{y}_{n+i} + \textbf{u}_{n}\\ \end{aligned} unN(0,Um)\begin{aligned} \textbf{u}_{n} &\sim N(\textbf{0}, \textbf{U}_{m}) \end{aligned}

1) Set:

V^0=U^0=C0\begin{aligned} \hat{\textbf{V}}_{0} &= \hat{\textbf{U}}_{0} = \textbf{C}_{0} \end{aligned}

And compute:

AIC0=N(klog2π+logV^0+k)+k(k+1)\begin{aligned} \text{AIC}_{0} &= N(k\log2\pi + \log |\hat{\textbf{V}}_{0}| + k) + k(k + 1) \end{aligned}

For m=1,,Mm = 1, \cdots, M:

a) Set:

Wm=Cmi=1m1Cmi\begin{aligned} \textbf{W}_{m} &= \textbf{C}_{m}- \sum_{i=1}^{m-1}\textbf{C}_{m-i} \end{aligned}

b) Obtain the PARCOR matrices of the forward and backward AR models by:

Amm=WmUm11Bmm=WmVm11\begin{aligned} \textbf{A}_{m}^{m} &= \textbf{W}_{m}\textbf{U}_{m-1}^{-1}\\ \textbf{B}_{m}^{m} &= \textbf{W}_{m}\textbf{V}_{m-1}^{-1}\\ \end{aligned}

c) Compute the AR coefficients of the forward and backward AR models by:

Amm=Aim1AmmBmim1Bmm=Bim1BmmAmim1\begin{aligned} \textbf{A}_{m}^{m} &= \textbf{A}_{i}^{m-1} - \textbf{A}_{m}^{m}\textbf{B}_{m-i}^{m-1}\\ \textbf{B}_{m}^{m} &= \textbf{B}_{i}^{m-1} - \textbf{B}_{m}^{m}\textbf{A}_{m-i}^{m-1} \end{aligned} i=1,,m1\begin{aligned} i &= 1, \cdots, m-1 \end{aligned}

d) Compute the innovation variance covariance matrices by:

Vm=C0i=1mAimCiTUm=C0i=1mBimCi\begin{aligned} \textbf{V}_{m} &= \textbf{C}_{0} - \sum_{i=1}^{m}\textbf{A}_{i}^{m}\textbf{C}_{i}^{T}\\ \textbf{U}_{m} &= \textbf{C}_{0} - \sum_{i=1}^{m}\textbf{B}_{i}^{m}\textbf{C}_{i} \end{aligned}

e) Compute the AIC of the AR model of order m by:

AICm=N(klog2π+logV^m+k)+k(k+1)+2k2m\begin{aligned} \text{AIC}_{m} &= N(k \log 2\pi + \log |\hat{\textbf{V}}_{m}| + k) + k(k + 1) + 2k^{2}m \end{aligned}

Example

yy <- HAKUSAN |>
  spread(key = key, value = value) |>
  as_tibble() |>
  select(YawRate, Rolling, Rudder) |>    
  as.matrix()

nc <- dim(yy)[1]
n <- seq(1, nc, by = 2)
y <- yy[n, ]

fit <- marfit(y, lag = 20)

> names(fit)
[1] "maice.order"  "aic"  "v" "arcoef"

Estimation of Multivariate AR Model by the Least Squares Method

This least square method has a significant advantage in that we can select a different order for each variable to enable more flexible modeling than the Yule-Walker method. In addition, it is also possible to specify a time-lag for a variable with respect to other variables or to specify a particular coefficient to be zero.

The function marlsq of package TSSS fits multivariate AR models by the least squares method based on the Householder transformation. The highest order of the AR model can be specified by the argument lag. If it is not given, the default value is 2N2\sqrt{N}.

Multivariate AR model was estimated for four-variate HAKUSAN data consisted of yaw rate, rolling, pitching and rudder angle by putting lag=10. The minimum AIC, 6575.741 was attained at m = 10. Using the estimated AR coefficient matrices arcoef and the innovation variance covariance matrix v, the cross-spectrum, coherency and the power contribution are obtained by the function marspc:

yy <- HAKUSAN |>
  spread(key = key, value = value) |>
  as_tibble() |>
  select(Pitching, YawRate, Rolling, Rudder) |>    
  as.matrix()

nc <- dim(yy)[1]
n <- seq(1, nc, by = 2)
y <- yy[n, ]

fit <- marlsq(y, lag = 10)
marspc(fit$arcoef, v = fit$v)

Locally Stationary AR Model

Records of real-world phenomena can mostly be categorized as nonstationary time series. The simplest approach to modeling nonstationary time series is to partition the time interval into several subintervals of appropriate size, on the assumption that the time series are stationary on each subinterval. Then, by fitting an AR model to each subinterval, we can obtain a series of models that approximate nonstationary time series.

It is assumed that the time series y1,,yNy_{1}, \cdots, y_{N} is nonstationary as a whole, but that we can consider it to be stationary on each subinterval of an appropriately constructed partition (N1++Nk=N)(N_{1} + \cdots + N_{k} = N). Such a time series that satisfies piecewise stationarity is called a locally stationary time series.

A locally stationary AR model is a nonstationary time series model which has the property that, on each appropriately constructed subinterval, it is stationary and can be modeled by an AR model on each of these subintervals.

Consider the i-th subinterval, [ni0,ni1][n_{i0}, n_{i1}]:

ni0=j=1i1Nj+1ni1=j=1iNj\begin{aligned} n_{i0} &= \sum_{j=1}^{i-1}N_{j} +1\\ n_{i1} &= \sum_{j=1}^{i}N_{j} \end{aligned}

For a locally stationary AR model, the time series yny_{n} follows an AR model on the j-th subinterval:

yn=i=1mjajiyni+vnj\begin{aligned} y_{n} &= \sum_{i=1}^{m_{j}}a_{ji}y_{n-i} + v_{nj} \end{aligned} E[vnj]=0E[vnj2]=σj\begin{aligned} E[v_{nj}] &= 0\\ E[v_{nj}^{2}] &= \sigma_{j}%+^{2} \end{aligned} E[vnjynm]=0\begin{aligned} E[v_{nj}y_{n-m}] &= 0 \end{aligned} m=1,2,\begin{aligned} m &= 1, 2, \cdots \end{aligned}

The best locally stationary AR model can, in principle, be obtained by the least squares method and AIC. However, practically speaking, it would require an enormous amount of computation to find the model that minimizes the AIC by fitting locally stationary AR models for all possible combinations of numbers of subintervals, k, and data lengths, N1,,NkN_{1}, \cdots, N_{k}. In terms of practice, the following procedure was developed to determine the dividing points of the locally stationary AR model.

Accordingly, only the points ni=iLn_{i} = iL are considered as candidates for dividing points, while the minimum unit L of division has been set beforehand. Then, the dividing points of the locally stationary AR model can automatically be decided by the following procedure:

1) Determine the basic span L and the highest order m of the AR models that are fitted to the subinterval of the length L. Here L is set to an appropriate length so that an AR model of order m can be fitted on an interval of length L.

2) Fit AR models of orders up to m to the time series y1,,yLy_{1}, \cdots, y_{L}, and compute AIC0(0),,AICm(0)\text{AIC}_{0}(0), \cdots, \text{AIC}_{m}(0) to find AIC0=minjAIC0(j)\text{AIC}_{0} = \mathrm{min}_{j} \text{AIC}_{0}(j). Further, set:

k=1n10=m+1n11=LN1=Lm\begin{aligned} k &= 1\\ n_{10} &= m + 1\\ n_{11} &= L\\ N_{1} &= L − m \end{aligned}

3) Fit AR models with orders up to m to the time series ynk1+1,,ynk1+Ly_{n_{k1} +1}, \cdots, y_{n_{k1} +L} and compute AIC1(0),,AIC1(m)\text{AIC}_{1}(0), \cdots, \text{AIC}_{1}(m) to set AIC1=minjAIC1(j)\text{AIC}_{1} = \textrm{min}_{j} \text{AIC}_{1}(j).

AIC1\text{AIC}_{1} is the AIC of a new model that was obtained under the assumption that the model changed at time nk1+1n_{k1} + 1. The AIC of the locally stationary AR model that divides the interval [nk0,nk1+L][n_{k0}, n_{k1} + L] into two subintervals [nk0,nk1][n_{k0}, n_{k1}] and [nk1+1,nk1+L][n_{k1} + 1, n_{k1} + L] is given by:

AICD=AIC0+AIC1\begin{aligned} \text{AIC}_{D} &= \text{AIC}_{0} + \text{AIC}_{1} \end{aligned}

This model is called a divided model.

4) Considering ynk0,,ynk1+Ly_{n_{k0}}, \cdots, y_{n_{k1}} + L to be a stationary interval, fit AR models of orders up to m to compute AICP(0),,AICP(m)\text{AIC}_{P}(0), \cdots, \text{AIC}_{P}(m), and then put AICP=minjAICP(j)\text{AIC}_{P} = \textrm{min}_{j} \text{AIC}_{P}(j). On the assumption that the time series on the entire interval [nk0,nk1+L][n_{k0}, n_{k1} + L] is stationary, the model is called a pooled model.

5) To judge the homogeneity of the two subintervals, compare the AICD\text{AIC}_{D} of the model of step 3 and the AICP\text{AIC}_{P} of the model of step 4.

If AICD<AICP\text{AIC}_{D} < \text{AIC}_{P}, it is judged that a divided model is better. In this case, nk1+1n_{k1} + 1 becomes the initial point of the current subinterval and we put:

k=k+1nk0=nk1,1+1nk1=nk1,1+L\begin{aligned} k &= k + 1\\ n_{k0} &= n_{k−1,1} + 1\\ n_{k1} &= n_{k−1,1} + L \end{aligned} Nk=L\begin{aligned} N_{k} &= L \end{aligned}

If AICDAICP\text{AIC}_{D} \geq \text{AIC}_{P}, a pooled model is adopted. In this case, the new subinterval [nk1+1,nk1+L][n_{k1} + 1, n_{k1} + L] is merged with the former subinterval, and [nk0,nk1+L][n_{k0}, n_{k1} + L] becomes the new current subinterval. Then we put:

nk1=nk1+LNk=Nk+L\begin{aligned} n_{k1} &= n_{k1} + L\\ N_{k} &= N_{k} + L \end{aligned} AIC0=AICP\begin{aligned} \text{AIC}_{0} &= \text{AIC}_{P} \end{aligned}

6) If we have at least L remaining additional data points, we go back to step (3).

Example: (Locally stationary modeling of seismic data)

The seismogram record involves microtremors as the noise and two types of seismic wave; the P-wave and the S-wave.

The function lsar of the package TSSS fits locally stationary AR model and automatically divides the time series into several subintervals on which the time series can be considered as stationary.

We fit a locally stationary AR model to the east-west component of a seismogram (N = 2600) with L = 150 and m = 10:

fit <- lsar(
  as.ts(MYE1F),
  max.arorder = 10,
  ns0 = 150
)

> names(fit)
[1]  "model"  "ns"     "span"   "nf"     "ms"     
[6]  "sds"    "aics"   "mp"     "sdp"    "aicp"    
[11] "spec"   "tsname"

The arguments to be specified are as follows:

  • max.arorder: highest order of AR models
  • ns0: basic local span.

The outputs from this function are:

  • model: = 1: pooled model is accepted, = 2: switched model is accepted.
  • ns: number of observations of local spans.
  • span: start points and end points of local spans.
  • nf: number of frequencies in computing power spectra.
  • ms: order of switched model.
  • sds: innovation variance of switched model.
  • aics: AIC of switched model.
  • mp: order of pooled model.
  • sdp: innovation variance of pooled model.
  • aicp: AIC of pooled model.
  • spec: local spectrum.

Plotting the frequency spectrums:

plts <- list()
for (i in seq_along(fit$spec$subspec)) {     
  spec <- fit$spec$subspec[[i]]
  n <- length(spec)

  plt <- tibble(
    freq = seq(0, 0.5, length.out = n),
    spec = spec
  ) |>
    ggplot(aes(x = freq, y = spec)) +
    geom_line() +
    ggtitle(
      paste0(
        fit$spec$start[i],
        " - ",
        fit$spec$end[i])
    ) +
    theme_tq()

  plts[[i]] <- plt
}

Structural changes have been detected at 8 points, n = 310, 610, 760, 910, 1060, 1360, 1660 and 2010. The change around n = 600 corresponds to a change in the spectrum and the variance caused by the arrival of the P-wave. The section n = 610−910 corresponds to the P-wave. Whereas the spectrum during n = 610 − 760 contains a single strong periodic component, various periodic components with different periods are intermingled during the latter half of the P-wave, n = 760 − 1060. The S-wave appears after n = 1060. We can see not only a decrease in power due to the reduction of the amplitude but also that the main peak shifts from the low frequency range to the high frequency range. After n = 2110, no significant change in the spectrum could be detected.

Precise Estimation of the Change Point

Here, we consider a method of detecting the precise time of a structural change by assuming that a structural change of the time series yny_{n} occurred within the time interval [n0,n1][n_{0}, n_{1}].

Assuming that the structural change occurred at time n;n0nn1n; n_{0} \leq n \leq n_{1}, a different AR model is fitted to each subinterval [1,n1][1, n − 1] and [n,N][n, N], respectively. Then the sum of the two AIC’s of the AR models fitted to these time series yields the AIC of a locally stationary AR model with a structural change at time n.

To obtain a precise estimate of the time of structural change based on the locally stationary AR models, we compute the AIC’s for all n such that n0nn1n_{0} \leq n \leq n_{1} to find the minimum value. With this method, because we have to estimate AR models for all n, a huge amount of computation is required. However, we can derive a computationally very efficient procedure for obtaining the AIC’s for all the locally stationary AR models by using the method of data augmentation.

Example: Estimation of the arrival times of P-wave and S-wave

The function lsar.chgpt of the package TSSS provides a precise estimate of the time of the structural change of the time series. The required parameters for this function are as follows:

  • max.arorder: highest order of AR model.
  • subinterval: a vector of the form c(n0,ne), which specifies start and end points of time interval used for modeling.
  • candidate: a vector of the form c(n1, n2), which gives the minimum and the maximum for change point. n0+2k < n1 < n2+k < ne, (k is max.arorder)

The outputs from this function are:

  • aic: AIC’s of the AR model fitted on [n1,n2].
  • aicmin: minimum AIC.
  • change.point: estimated change point.
  • subint: original sub-interval data.

The first fit exames the change points around n=600n = 600:

fit <- lsar.chgpt(
  as.ts(MYE1F),
  max.arorder = 10,
  subinterval = c(200,1000),    
  candidate = c(400,800)
)

> names(fit)
[1] "aic"  "aicmin"  "change.point"  "subint"    

> fit$change.point
[1] 630

Plotting the sub-interval from 400 to 800:

MYE1F |>
  slice(fit$subint$st:fit$subint$ed) |>    
  autoplot(value) +
  xlab("") +
  ylab("") +
  theme_tq()

And plotting the change point coinciding with the minimum AIC:

tsibble(
  aic = fit$aic,
  t = fit$subint$st:fit$subint$ed,
  index = t
) |>
  autoplot(aic) +
  geom_text(
    aes(
      y = fit$aicmin,
      x = fit$change.point,
      label = "AIC Min",
      col = "red"
    )
  ) +
  xlab("") +
  ylab("") +
  theme_tq() +
  theme(legend.position = "none")    

The first half before the change point is considered the micotremors and the latter half is considered the P-wave.

The next change point is around n=1000n = 1000:

fit <- lsar.chgpt(
  as.ts(MYE1F),
  max.arorder = 10,
  subinterval = c(600,1400),    
  candidate = c(800,1200)
)

> fit$change.point
[1] 1026

The first half before the second change point is considered the P-wave and latter half the S-wave.

Posterior Probability of the Change Point

It can be shown the likelihood of the model whose parameters are estimated by the maximum likelihood method:

p(yj)=e12AICj\begin{aligned} p(y|j) &= e^{-\frac{1}{2}\text{AIC}_{j}} \end{aligned}

The above is the likelihood of the LSAR model, which assumes that jj is the arrival time.

Therefore, if the prior distribution of the arrival time is given, then the posterior distribution of the arrival time can be obtained by:

p(jy)=p(yj)p(j)jp(yj)p(j)\begin{aligned} p(j|y) &= \frac{p(y|j)p(j)}{\sum_{j}p(y|j)p(j)} \end{aligned}

In the example shown in the next section, the uniform prior over the interval is used. It seems more reasonable to put more weight on the center of the interval. However, since the likelihood, p(yj)p(y\mid j), usually takes significant values only at narrow areas, only the local behaviour of the prior is influential to the posterior probability. Therefore as long as very smooth functions are used, the choice of the prior is not so important in the present problem.

fit1 <- lsar.chgpt(
  as.ts(MYE1F),
  max.arorder = 10,
  subinterval = c(400, 800),
  candidate = c(600, 700)
)
AICP <- fit$aic
post <- exp(-(AICP - min(AICP)) / 2)
post <- post / sum(post)

fit2 <- lsar.chgpt(
  as.ts(MYE1F),
  max.arorder = 10,
  subinterval = c(600, 1400),
  candidate = c(800, 1200)
)
AICP2 - fit2$aic
post2 <- exp(-(AICP2 - min(AICP2)) / 2)
post2 <- post / sum(post2)

In the case of the arrival of the P-wave, the AIC shown in the middle plot has a clear minimum at about 630. On the other hand, in the case of the arrival of the S-wave, change of the AIC is not so significant and actually has several local minima. Note the difference of the scale of vertical axes in the middle plots. As a result, posterior distribution of the P-wave shown in the left plot has a single peak. On the other hand, that of the S-wave has several peaks.

See Also

References

Kitagawa G. (2010) Introduction to Time Series Modeling

Jason

Passionate software developer with a background in CS, Math, and Statistics. Love challenges and solving hard quantitative problems with interest in the area of finance and ML.