Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC



Objectives

  1. Putting autoregressive moving average (ARMA) models into the context of linear time series models.

  2. Introduce the backshift operator, and see how it can be used to develop an algebraic approach to studying the properties of ARMA models.




4.1 Definition: Stationary causal linear process.




4.1.1 Question: When does the “stationary” here mean weak stationarity, and when does it mean strict stationary?




  • causal refers to \(\{\epsilon_n\}\) being a causal driver of \(\{Y_n\}\). The value of \(Y_n\) depends only on noise process values already determined by time \(n\). This matching a requirement for causation that causes must precede effects.

  • linear refers to linearity of \(Y_n\) as a function of \(\{\epsilon_n\}\). A linear modification of the noise process, replacing \(\{\epsilon_n\}\) by \(\{\alpha + \beta\epsilon_n\}\), results in a new linear process.

  • The autocovariance function is, \[\begin{eqnarray} \gamma_h &=& {\mathrm{Cov}}\big(Y_n,Y_{n+h}\big)\\ &=& {\mathrm{Cov}}\left(\sum_{j=0}^\infty g_j\epsilon_{n-j},\sum_{k=0}^\infty g_k\epsilon_{n+h-k}\right)\\ &=& \sum_{j=0}^\infty \sum_{k=0}^\infty g_j g_k{\mathrm{Cov}}\big(\epsilon_{n-j},\epsilon_{n+h-k}\big)\\ &=& \sum_{j=0}^\infty g_jg_{j+h} \sigma^2, \mbox{for $h\ge 0$}. \end{eqnarray}\]

  • Here, we have assumed we can move \(\sum_{j=0}^\infty \sum_{k=0}^\infty\) through \({\mathrm{Cov}}\). The interchange of expectation and infinite sums can’t be taken for granted. In this course, we do not focus on interchange issues, but we try to notice when we make assumptions.

  • In order for this autocovariance function to exist, we need \[\sum_{j=0}^\infty g_j^2 < \infty.\]

  • The interchange can be justified by requiring a stronger condition, \[\sum_{j=0}^\infty |g_j| < \infty.\]

  • The MA(q) model that we defined in equation [M3] is an example of a stationary, causal linear process.

  • The general stationary, causal linear process model [M7] can also be called the MA(\(\infty\)) model.




4.1.2 A stationary causal linear solution to the AR(1) model, and a non-causal solution

  • Recall the stochastic difference equation defining the AR(1) model,

    [M8] \({\quad\quad\quad}Y_n = {\phi}Y_{n-1}+\epsilon_n\).

  • This has a causal solution,

    [M8.1] \({\quad\quad\quad}Y_n = \sum_{j=0}^\infty {\phi}^j\epsilon_{n-j}\).

  • It also has a non-causal solution,

    [M8.1] \({\quad\quad\quad}Y_n = \sum_{j=0}^\infty {\phi}^{-j}\epsilon_{n+j}\).




4.1.3 Exercise: Work through the algebra to check that M8.1 and M8.2 both solve equation [M8].




4.1.4 Question: Assess the convergence of the infinite sums in [M8.1] and [M8.2].

  • For what values of \({\phi}\) is the causal solution [M8.1] a convergent infinite sum, meaning that it converges to a random variable with finite variance? For what values is the non-causal solution [M8.2] a convergent infinite sum?




4.1.5 Exercise: Using the MA(\(\infty\)) representation to compute the autocovariance of an ARMA model

  • The linear process representation can be a convenient way to calculate autocovariance functions. Use the linear process representation in [M8.1], together with our expression for the autocovariance of the general linear process [M7], to get an expression for the autocovariance function of the AR(1) model.




4.2 The backshift operator and the difference operator




4.3 The general ARMA model




4.3.1 Exercise: Carry out the following steps to obtain the MA(\(\infty\)) representation and the autocovariance function of the ARMA(1,1) model,

\[ Y_n = {\phi}Y_{n-1} + \epsilon_n + {\psi}\epsilon_{n-1}.\]

  1. Formally, we can write \[ (1-{\phi}B)Y_n = (1+{\psi}B)\epsilon_n,\] which algebraically is equivalent to \[Y_n = \left(\frac{1+{\psi}B}{1-{\phi}B}\right)\epsilon_n.\] We write this as \[Y_n = g(B) \epsilon_n,\] where \[ g(x) = \left(\frac{1+{\psi}x}{1-{\phi}x}\right).\] To make sense of \(g(B)\) we need to work out the Taylor series expansion, \[g(x) = g_0 + g_1 x + g_2 x^2 + g_3 x^3 + \dots\] Do this either by hand or using your favorite math software.

  2. From 1. we can get the MA(\(\infty\)) representation. Then, we can apply the general formula for the autocovariance function of an MA(\(\infty\)) process.




4.4 Causal, invertible ARMA models

roots <- polyroot(c(1,2,2))
roots
## [1] -0.5+0.5i -0.5-0.5i

Finding the absolute value shows that we have two roots inside the unit circle, so this MA(2) model is not invertible.

abs(roots)
## [1] 0.7071068 0.7071068




4.4.1 Question: It is undesirable to use a non-invertible model for data analysis. Why?

One answer to this question involves diagnosing model misspecification.




4.5 Reducible and irreducible ARMA models

list(AR_roots=polyroot(c(1,-5/6,1/6)),MA_roots=polyroot(c(1,-1,1/4)))
## $AR_roots
## [1] 2+0i 3+0i
## 
## $MA_roots
## [1] 2+0i 2-0i




4.6 Deterministic skeletons: Using differential equations to study ARMA models

4.6.1 Example: oscillatory behavior modeled using an AR(2) process

  • In physics, a basic model for processes that oscillate (springs, pendulums, vibrating machine parts, etc) is simple harmonic motion.

  • The differential equation for a simple harmonic motion process \(x(t)\) is

    [M10] \({\quad\quad\quad}\frac{d^2}{dt^2} x(t) = -\omega^2 x(t)\).

  • This is a second order linear differential equation with constant coefficients. Such equations have a closed form solution, which is fairly straightforward to compute once you know how!

  • The solution method is very similar to the method for solving difference equations coming up elsewhere in time series analysis, so let’s see how it is done.

  1. Look for solutions of the form \(x(t)=e^{\lambda t}\). Substituting this into the differential equation [M10] we get \[ \lambda^2 e^{\lambda t} = -\omega^2 e^{\lambda t}.\] Canceling the term \(e^{\lambda t}\), we see that this has two solutions, with \[ \lambda = \pm \omega i,\] where \(i=\sqrt{-1}\).

  2. The linearity of the differential equation means that if \(y_1(t)\) and \(y_2(t)\) are two solutions, then \(Ay_1(t)+By_2(t)\) is also a solution for any \(A\) and \(B\). So, we have a general solution to [M10] given by \[ x(t) = A e^{i\omega t} + B e^{-i\omega t}.\]

  3. Using the two identities, \[\sin(\omega t) = \frac{1}{2}\big(e^{i\omega t} - e^{-i\omega t}\big), \quad\quad \cos(\omega t) = \frac{1}{2}\big(e^{i\omega t} + e^{-i\omega t}\big), \] we can rewrite the general solution as \[ x(t) = A \sin(\omega t) + B\cos(\omega t),\] which can also be written as \[ x(t) = A\sin(\omega t + \beta).\] For the solution in this form, \(\omega\) is called the frequency, \(A\) is called the amplitude of the oscillation and \(\beta\) is called the phase. The frequency of the oscillation is determined by [M10], but the amplitude and phase are unspecfied constants. Initial conditions can be used to specify \(A\) and \(\beta\).

  • A discrete time version of [M10] is a deterministic linear difference equation, replacing \(\frac{d^2}{dt^2}\) by the second difference operator, \(\Delta^2 = (1-B)^2\). This corresponds to a deterministic model equation, \[{\quad\quad\quad}\Delta^2 y_n = - \omega^2 y_n.\]

  • Adding white noise, and expanding out \(\Delta^2 = (1-B)^2\), we get a stochastic model,

    [M11] \({\quad\quad\quad}Y_n = \frac{2Y_{n-1}}{1+\omega^2} - \frac{Y_{n-2}}{1+\omega^2} + \epsilon_n\).

  • It seems reasonable to hope that model [M11] would be a good candidate to describe systems that have semi-regular but somewhat eratic fluctuations, called quasi-periodic behavior. Such behavior is evident in business cycles or wild animal populations.

  • Let’s look at a simulation from [M11] with \(\omega=0.1\) and \(\epsilon_n\sim \mbox{IID } N[0,1]\). From our exact solution to the deterministic skeleton, we expect that the period of the oscillations (i.e., the time for each completed oscillation) should be approximately \(2\pi/\omega\).

omega <- 0.1
ar_coefs <- c(2/(1+omega^2), - 1/(1+omega^2))
set.seed(8395200)
X <- arima.sim(list(ar=ar_coefs),n=500,sd=1)
par(mfrow=c(1,2))
plot(X)
plot(ARMAacf(ar=ar_coefs,lag.max=500),type="l",ylab="ACF of X")

  • Quasi-periodic fluctuations are said to be “phase locked” as long as the random peturbations are not able to knock the oscillations away from being close to their initial phase.

  • Eventually, the randomness should mean that the process is equally likely to have any phase, regardless of the initial phase.



4.6.2 Question: What is the timescale on which the simulated model shows phase locked behavior?

  • Equivalently, on what timescale does the phase of the fluctuations lose memory of its initial phase?