Chapter 9 Topics in Time Series Analysis

9.1.1 The AR(1) Model

The autoregressive model of order 1 is defined next. It has a closed-form expression for the population autocorrelation function and is frequently used in applications.

No subscript is necessary on the [latex]\phi[/latex] parameter because there is only one [latex]\phi[/latex] parameter in the AR(1) model. So there are two parameters that define an AR(1) model: the coefficient [latex]\phi[/latex] and the population variance of the white noise [latex]\sigma _ Z ^ {\, 2}[/latex].

The current value in the time series, X_t, is given by the parameter [latex]\phi[/latex] multiplied by the previous observed value in the time series, [latex]\phi X_{t - 1}[/latex], plus the current white noise term Z_t. This model has the form of a simple linear regression model forced through the origin in which X_t is being predicted by the previous value of the time series [latex]X_{t-1}[/latex]. The parameter [latex]\phi[/latex] plays the role of the slope of the regression line. Thinking about an AR(1) model as a simple linear regression model suggests a statistical graphic that can be helpful in determining whether it is an appropriate model for a particular time series. A plot of x_t on the vertical axis against [latex]x_{t - 1}[/latex] on the horizontal axis should be approximately linear if the AR(1) model is appropriate. The slope of the regression line on this plot corresponds to [latex]\phi[/latex], and the magnitude of the variability of the points about the regression line is determined by the population variance of the white noise [latex]\sigma _ Z ^ {\, 2}[/latex].

Some authors prefer to parameterize the AR(1) model as

$$\begin{array}{l}{X}_t={\varphi }_1{X}_{t-1}+{\varphi }_0{Z}_t,\end{array}$$

where [latex]\phi_0[/latex] and [latex]\phi_1[/latex] are real-valued parameters. We avoid this parameterization because the [latex]\phi_0[/latex] parameter is redundant in the sense that the population variance of the white noise [latex]\sigma _Z ^ {\, 2}[/latex] is absorbed into the [latex]\phi_0[/latex] parameter. Also, some authors use a – rather than a + between the terms on the right-hand side of the model.

To illustrate the thinking behind the AR(1) model in a specific context, let X_t represent the closing price of a particular stock on day t. The AR(1) model indicates that today’s closing price, denoted by X_t, equals [latex]\phi[/latex] multiplied by yesterday’s closing price ([latex]\phi X_{t - 1}[/latex]), plus today’s random white noise term Z_t.

Stationarity

One initial important question concerning the AR(1) model is whether or not the model is stationary. Consider a thought experiment that determines whether an AR(1) model is stationary for specific values of [latex]\phi[/latex]. For one particular instance, consider [latex]\phi = 0[/latex]. In this case the AR(1) time series model reduces to

$$\begin{array}{l}{X}_t=Z_t,\end{array}$$

which is a time series model consisting solely of white noise. We know from Example 7.15 that a time series model of white noise terms is stationary. Now consider another instance, [latex]\phi = 1[/latex]. In this case the AR(1) time series model reduces to

$$\begin{array}{l}{X}_t=X_{t-1}+Z_t,\end{array}$$

which indicates that each value in the time series is the previous value plus the current white noise. In this case the population variance of the process is increasing with time because the number of white noise terms accumulate over time (see Example 7.8), so the AR(1) model with [latex]\phi = 1[/latex] violates one of the stationarity conditions given in Definition 7.6. The AR(1) model with [latex]\phi = 1[/latex] can be recognized as a random walk model from Example 7.4, and it was determined to be nonstationary in Example 7.17. So we have established that the AR(1) time series model is stationary for [latex]\phi = 0[/latex] and nonstationary for [latex]\phi = 1[/latex]. We now try to determine general restrictions on [latex]\phi[/latex] associated with a stationary AR(1) time series model. We take four different approaches to establishing the values of the coefficient [latex]\phi[/latex] that lead to a stationary model. The four approaches provide a review of several concepts defined previously.

Approach 1: Causality. In the derivations concerning the AR(1) time series model that follow, it will be beneficial to write the time series value X_t as a linear combination of the current and previous white noise values. This will allow us to use the definition of causality in Definition 8.2 to determine the values of [latex]\phi[/latex] associated with a stationary AR(1) model. To begin, recall that the AR(1) model given by

$$\begin{array}{l}{X}_t=\varphi X_{t-1}+Z_t\end{array}$$

can be shifted in time and is equally valid for other t values, for example,

$$\begin{array}{ll}{X}_{t-1}& =\varphi X_{t-2}+Z_{t-1}\\ X_{t-2}& =\varphi X_{t-3}+Z_{t-2}\\ & ⋮\end{array}$$

Using successive substitutions into the AR(1) model results in

$$\begin{array}{ll}{X}_t& =\varphi X_{t-1}+Z_t\\ & =\varphi (\varphi X_{t-2}+Z_{t-1})+Z_t\\ & ={\varphi }^2{X}_{t-2}+\varphi Z_{t-1}+Z_t\\ & ={\varphi }^2(\varphi X_{t-3}+Z_{t-2})+\varphi Z_{t-1}+Z_t\\ & ={\varphi }^3{X}_{t-3}+{\varphi }^2{Z}_{t-2}+\varphi Z_{t-1}+Z_t\\ & ⋮\\ & =Z_t+\varphi Z_{t-1}+{\varphi }^2{Z}_{t-2}+{\varphi }^3{Z}_{t-3}+\cdots .\end{array}$$

This can be recognized as an MA(∞) time series model. Representing an AR(1) model as an MA(∞) model is known as duality. We now determine the constraints on the parameter [latex]\phi[/latex] which are required for stationarity. This is the form that is required for causality from Definition 8.2. The coefficients [latex]\psi_1, \, \psi_2, \, \ldots[/latex] for the AR(1) model from Definition 8.2 are

$$\begin{array}{l}{\psi }_1=\varphi ,{\psi }_2={\varphi }^2,{\psi }_3={\varphi }^3,\dots ,\end{array}$$

or in general, [latex]\psi_j = \phi ^ j[/latex], for [latex]j = 1, \, 2, \, \ldots[/latex]. Definition 8.2 requires that

$$\begin{array}{l}\sum _{j=1}^{\mathrm{\infty }}{\psi }_j^2=\sum _{j=1}^{\mathrm{\infty }}{\varphi }^{2j}={\varphi }^2+{\varphi }^4+{\varphi }^6+\cdots <\mathrm{\infty }\end{array}$$

for the time series model to be written in causal form. This summation is a geometric series that converges when [latex]|\phi| < 1[/latex], or equivalently, when [latex]-1 < \phi < 1[/latex], so these are the values of [latex]\phi[/latex] for which the AR(1) model is causal, which also implies that the model is stationary. Expressing the AR(1) model as an MA(∞) model will also be helpful in the subsequent derivation of the population autocovariance and autocorrelation functions.

Approach 2: Backshift operator. Although the purely algebraic derivation of the causal form of the AR(1) time series model using standard algebra techniques from Approach 1 works fine for establishing stationarity, there is an alternative approach which is slightly more elegant that exploits the backshift operator B. The AR(1) model

$$\begin{array}{l}{X}_t=\varphi X_{t-1}+Z_t\end{array}$$

can be rewritten as

$$\begin{array}{l}{X}_t-\varphi X_{t-1}=Z_t,\end{array}$$

which can be expressed using the backshift operator as

$$\begin{array}{l}(1-\varphi B)X_t=Z_t.\end{array}$$

The first-order polynomial [latex]\phi(B) = 1 - \phi B[/latex] generalizes to a polynomial in B of order p for an AR(p) model. Dividing both sides of this equation by [latex]1 - \phi B[/latex] gives

$$\begin{array}{l}{X}_t=\frac{Z_t}{1-\varphi B}.\end{array}$$

For values of [latex]\phi[/latex] satisfying [latex]-1 < \phi < 1[/latex], this can be recognized as a geometric series in B:

$$\begin{array}{l}{X}_t=(1+\varphi B+{\varphi }^2{B}^2+{\varphi }^3{B}^3+\cdots )Z_t.\end{array}$$

Executing the B operator converts this to the form

$$\begin{array}{l}{X}_t=Z_t+\varphi Z_{t-1}+{\varphi }^2{Z}_{t-2}+{\varphi }^3{Z}_{t-3}+\cdots ,\end{array}$$

which is the same form that we encountered using the successive substitutions in the causality approach.

Approach 3: Unit roots analysis. Theorem 8.3 indicates that all AR(1) models are invertible and they are stationary when the root of

$$\begin{array}{l}\varphi (B)=1-\varphi B=0\end{array}$$

lies outside of the unit circle in the complex plane. The solution to this equation is

$$\begin{array}{l}B=\frac{1}{\varphi }.\end{array}$$

This root falls on the real axis in the complex plane and lies outside of the unit circle when

$$\begin{array}{l}-1<\varphi <1,\end{array}$$

which is consistent with Approaches 1 and 2.

Approach 4: Definition of stationarity. We can also return to first principles to establish the values of [latex]\phi[/latex] associated with a stationary AR(1) model. This approach also results in the derivation of the population autocorrelation function. Recall from Definition 7.6 that a time series model is stationary if (a) the expected value of X_t is constant for all t, and (b) the population covariance between X_s and X_t depends only on the lag [latex]|t - s|[/latex]. Using the causal formula for the AR(1) time series model expressed as an MA(∞) time series model from Approach 1, the expected value of X_t is

$$\begin{array}{ll}E\left[X_t\right]& =E[Z_t+\varphi Z_{t-1}+{\varphi }^2{Z}_{t-2}+{\varphi }^3{Z}_{t-3}+\cdots ]\\ & =E\left[Z_t\right]+\varphi E\left[\varphi Z_{t-1}\right]+{\varphi }^{2}E\left[Z_{t-2}\right]+{\varphi }^{3}E\left[Z_{t-3}\right]+\cdots \\ & =0\end{array}$$

for all values of the parameters [latex]\phi[/latex] and [latex]\sigma _ Z ^ {\, 2}[/latex], and all values of t. Again using the causal formula for the AR(1) time series model expressed as an MA(∞) time series model,

$$\begin{array}{ll}\gamma (s,t)& =\text{Cov}(X_s,X_t)\\ & =\text{Cov}(Z_s+\varphi Z_{s-1}+{\varphi }^2{Z}_{s-2}+\cdots ,Z_t+\varphi Z_{t-1}+{\varphi }^2{Z}_{t-2}+\cdots )\\ & =\text{Cov}(Z_s,{\varphi }^{|t-s|}{Z}_s)+\text{Cov}(\varphi Z_{s-1},{\varphi }^{|t-s|+1}{Z}_{s-1})+\text{Cov}({\varphi }^2{Z}_{s-2},{\varphi }^{|t-s|+2}{Z}_{s-2})+\cdots \\ & ={\varphi }^{|t-s|}{\sigma }_Z^2+{\varphi }^{|t-s|+2}{\sigma }_Z^2+{\varphi }^{|t-s|+4}{\sigma }_Z^2+\cdots \\ & =({\varphi }^{|t-s|}+{\varphi }^{|t-s|+2}+{\varphi }^{|t-s|+4}+\cdots ){\sigma }_Z^2\\ & ={\varphi }^{|t-s|}\left(\frac{1}{1-{\varphi }^2}\right){\sigma }_Z^2|t-s|=0,1,2,\dots \end{array}$$

for [latex]-1 < \phi < 1[/latex]. Since [latex]E \left[ X_t \right] = 0[/latex] for all values of t and the population autocovariance function depends only on the lag [latex]|t - s|[/latex], we conclude that the AR(1) process is stationary when [latex]-1 < \phi < 1[/latex].  So the population autocovariance function can be expressed in terms of the lag k as

$$\begin{array}{l}\gamma (k)=\left(\frac{{\varphi }^k}{1-{\varphi }^2}\right){\sigma }_Z^{2}k=0,1,2,\dots .\end{array}$$

Dividing by the population autocovariance function by

$$\begin{array}{l}\gamma (0)=\left(\frac{1}{1-{\varphi }^2}\right){\sigma }_Z^2\end{array}$$

gives the population autocorrelation function

$$\begin{array}{l}\rho (k)=\frac{\gamma (k)}{\gamma (0)}={\varphi }^{k}k=0,1,2,\dots .\end{array}$$

Based on the four approaches, we now know beyond a shadow of doubt that an AR(1) model is stationary for values of the parameter [latex]\phi[/latex] satisfying [latex]-1 < \phi < 1[/latex]. This derivation constitutes a proof of the following result, which will be stated for just the nonnegative lags. Many authors list the lags as [latex]k = \pm 1, \, \pm 2, \, \ldots \,[/latex], but we appeal to Theorem 7.1 to cover the negative lags and only report the nonnegative lags in all of the population autocorrelation functions given in this chapter.

The derivation of [latex]\rho(k) = \phi ^ k[/latex] for [latex]k = 0, \, 1, \, 2, \, \ldots[/latex] provides still further evidence of the restriction that [latex]-1 < \phi < 1[/latex]. If [latex]\phi[/latex] were equal to a value outside of this range, say [latex]\phi = 2[/latex],  this would result in population correlation values outside of the range [latex]-1 \le \rho(k) \le 1[/latex].

For all admissible values of [latex]\phi[/latex] on the interval [latex]-1 < \phi < 1[/latex], we see from the formula [latex]\rho(k) = \phi ^ k[/latex] for [latex]k = 0, \, 1, \, 2, \, \ldots[/latex] that there will be a geometric decline in the magnitude of the values in the population autocorrelation function as the lag k increases. There are two distinct cases for [latex]\phi[/latex], however, which will result in population autocorrelation functions with distinctly different shapes. The first case is [latex]0 < \phi < 1[/latex], which gives positive population autocorrelation values at all lags. This is associated with a time series that lingers on one side of the mean. How long it lingers depends on the magnitude of [latex]\phi[/latex]. Larger values of [latex]\phi[/latex] indicate that nearby observations will tend to be more likely to be on the same side of the mean, and therefore the time series will tend to linger longer on one side of the mean. The second case is [latex]-1 < \phi < 0[/latex], which gives population autocorrelation function values which alternate in sign and is associated with a time series that is likely to jump from one side of the mean to the other for adjacent observations. These two cases are illustrated in Figure 9.1 for the first 8 lags of the population autocorrelation function for [latex]\phi = 0.8[/latex] and [latex]\phi = -0.8[/latex].

Population Partial Autocorrelation Function

We now determine the population partial autocorrelation function for an AR(1) model. By Definition 7.4, the population lag 0 partial autocorrelation value is [latex]\rho ^ * (0) = 1[/latex]. The population lag 1 partial autocorrelation value is [latex]\rho ^ * (1) = \rho (1) = \phi[/latex]. The population lag 2 partial autocorrelation is

$$\begin{array}{l}{\rho }^{\ast }(2)=\frac{\left|\begin{array}{cc}1& \rho (1)\\ \rho (1)& \rho (2)\end{array}\right|}{\left|\begin{array}{cc}1& \rho (1)\\ \rho (1)& 1\end{array}\right|}=\frac{\left|\begin{array}{cc}1& \varphi \\ \varphi & {\varphi }^2\end{array}\right|}{\left|\begin{array}{cc}1& \varphi \\ \varphi & 1\end{array}\right|}=0.\end{array}$$

This is consistent with the result from Example 7.22 from first principles. Notice that the second column of the matrix in the numerator is a multiple of the first column of the matrix in the numerator. This is why the determinant of the numerator is zero. The population lag 3 partial autocorrelation is

$$\begin{array}{l}{\rho }^{\ast }(3)=\frac{\left|\begin{array}{ccc}1& \rho (1)& \rho (1)\\ \rho (1)& 1& \rho (2)\\ \rho (2)& \rho (1)& \rho (3)\end{array}\right|}{\left|\begin{array}{ccc}1& \rho (1)& \rho (2)\\ \rho (1)& 1& \rho (1)\\ \rho (2)& \rho (1)& 1\end{array}\right|}=\frac{\left|\begin{array}{ccc}1& \varphi & \varphi \\ \varphi & 1& {\varphi }^2\\ {\varphi }^2& \varphi & {\varphi }^3\end{array}\right|}{\left|\begin{array}{ccc}1& \varphi & {\varphi }^2\\ \varphi & 1& \varphi \\ {\varphi }^2& \varphi & 1\end{array}\right|}=0.\end{array}$$

Again, the determinant in the numerator is zero because the third column is a multiple of the first column. This pattern continues for the lag k population partial autocorrelation function, which has a first column of the numerator matrix [latex]\left[ 1, \, \phi , \, \phi ^ 2 , \, \ldots , \, \phi ^ {k - 1} \right] ^ \prime[/latex] and last column [latex]\left[ \phi, \, \phi ^ 2 , \, \phi ^ 3 , \, \ldots , \, \phi ^ k \right] ^ \prime[/latex]. Since the last column of the numerator matrix is a multiple of the first column of the numerator matrix, the determinant of the numerator matrix is zero. This constitutes a proof of the following result.

Figure 9.2 shows the first 8 lags of the population partial autocorrelation function for [latex]\phi = 0.8[/latex] and [latex]\phi = -0.8[/latex]. These are the same parameter settings as in Figure 9.1. Unlike the population autocorrelation function which tails off in magnitude for increasing lags, the population partial autocorrelation cuts off after lag 1. When plotting the corresponding sample analogs, it is typically easier to visually assess a function cutting off rather than tailing off, particularly if there is significant random sampling variability in the observed time series.

The Shifted AR(1) Model

The population mean function for the AR(1) model is [latex]E \left[ X_t \right] = 0[/latex]. This model is not of much use in practice because most real-world time series are not centered around zero. Adding a third parameter μ overcomes this shortcoming. Since population variance and covariance are unaffected by a shift, the associated population autocorrelation and partial autocorrelation functions remain the same as those given in Theorems 9.1 and 9.2. Likewise, the condition for stationarity is unchanged.

Simulation

An AR(1) time series can be simulated by appealing to the defining formula for the AR(1) model. Iteratively applying the defining formula for an AR(1) model

$$\begin{array}{l}{X}_t=\varphi X_{t-1}+Z_t\end{array}$$

results in the simulated values [latex]X_1, \, X_2, \, \ldots, \, X_n[/latex]. The difficult aspect of this algorithm is how to generate the first value X₁ because there is no X₀ available. For simplicity, assume that the white noise terms are Gaussian white noise. Since the expected value of X_t is [latex]E \left[ X_t \right] = 0[/latex], the population variance of X_t is

$$\begin{array}{l}V\left[X_t\right]=\gamma (0)=\left(\frac{1}{1-{\varphi }^2}\right){\sigma }_Z^2,\end{array}$$

and linear combinations of mutually independent normally distributed random variables are normal, then the first simulated observation

$$\begin{array}{l}{X}_1\sim N(0,\left(\frac{1}{1-{\varphi }^2}\right){\sigma }_Z^2).\end{array}$$

The algorithm given as pseudocode below generates an initial time series observation X₁ as indicated above, and then uses an additional [latex]n - 1[/latex] Gaussian white noise terms [latex]Z_2, \, Z_3, \, \ldots, \, Z_n[/latex] to generate the remaining time series values [latex]X_2, \, X_3, \, \ldots, \, X_n[/latex] using the AR(1) defining formula from Definition 9.1. Indentation denotes nesting in the algorithm.

The three-parameter shifted AR(1) time series model which includes a population mean parameter μ can be simulated by simply adding μ to each time series observation generated by this algorithm. The next example implements this algorithm in R.

We recommend running the simulation code from the previous example several dozen times in a loop and viewing the associated plots of x_t, r_k, and [latex]r_k^*[/latex] in search of patterns. A call to the R function Sys.sleep between the displays of the trio of plots can be used to include an artificial time delay to allow you to inspect the plots. This will allow you to see how various realizations of a simulated AR(1) time series model vary from one realization to the next. So when you then view a single realization of a real-life time series, you will have a better sense of how far these plots might deviate from their expected patterns.

There is a second way to simulate observations from an AR(1) time series. This second technique starts the time series at an initial arbitrary value, and then allows the time series to “warm up” or “burn in” for several time periods before producing the first observation X₁. A reasonable initial arbitrary value for the standard AR(1) model is 0; a reasonable initial arbitrary value for the shifted AR(1) model is μ. This is the approach taken by the built-in R function named arima.sim (for autoregressive moving average simulation), which simulates a realization of a time series. Using the arima.sim function saves a few keystrokes over the approach taken in the previous example, as illustrated next.

The remaining topics associated with the AR(1) time series model are statistical in nature: parameter estimation, model assessment, model selection, and forecasting. A sample time series that will be revisited throughout these topics will be introduced next.

Parameter Estimation

There are two parameters, [latex]\phi[/latex] and [latex]\sigma _ Z ^ {\, 2}[/latex], to estimate in the standard AR(1) model

$$\begin{array}{l}{X}_t=\varphi X_{t-1}+Z_t.\end{array}$$

There are three parameters, μ, [latex]\phi[/latex], and [latex]\sigma _ Z ^ {\, 2}[/latex], to estimate in the shifted AR(1) model

$$\begin{array}{l}{X}_t-\mu =\varphi (X_{t-1}-\mu )+Z_t.\end{array}$$

The three parameter estimation techniques outlined in Section 8.2.1 are applied to the shifted AR(1) time series model next.

Approach 1: Method of moments. In the case of the shifted AR(1) model, we match the population and sample (a) first-order moments, (b) second-order moments, and (c) lag 1 autocorrelation. Placing the population moments on the left-hand side of the equation and the associated sample moments on the right-hand side of the equation results in three equations in three unknowns:

$$\begin{array}{ll}E\left[X_t\right]& =\frac{1}{n}\sum _{t=1}^n{X}_t\\ E\left[X_t^2\right]& =\frac{1}{n}\sum _{t=1}^n{X}_t^2\\ \rho (1)& =r_1\end{array}$$

or

$$\begin{array}{rrr}\hfill \mu & \hfill =& \overline{X}\hfill \\ \hfill V\left[X_t\right]+E{\left[X_t\right]}^2=\gamma (0)+{\mu }^2=\left(\frac{1}{1-{\varphi }^2}\right){\sigma }_Z^2+{\mu }^2& \hfill =& \frac{1}{n}{\displaystyle \sum _{t=1}^n{X}_t^2}\hfill \\ \hfill \varphi & \hfill =& r_1.\hfill \end{array}$$

These equations can be solved in closed form for the three unknown parameters μ, [latex]\phi[/latex], and [latex]\sigma _ Z ^ {\, 2}[/latex] yielding the method of moments estimators

$$\begin{array}{l}\hat{\mu }=\overline{X},\hat{\varphi }=r_1,{\hat{\sigma }}_Z^2=(1-{\hat{\varphi }}^2)(\frac{1}{n}\sum _{t=1}^n{X}_t^2-{\hat{\mu }}^2)=(1-r_1^2)(\frac{1}{n}\sum _{t=1}^n{X}_t^2-{\overline{X}}^2).\end{array}$$

This constitutes a proof of the following result.

These point estimators are random variables and have been written as a function of the random time series values [latex]X_1, \, X_2, \, \ldots, \, X_n[/latex]. For observed time series values [latex]x_1, \, x_2, \, \ldots, \, x_n[/latex], the lowercase versions of the formulas will be used.

Approach 2: Least squares. Consider the shifted stationary AR(1) model

$$\begin{array}{l}{X}_t-\mu =\varphi (X_{t-1}-\mu )+Z_t.\end{array}$$

For least squares estimation, we first establish the sum of squares S as a function of the parameters μ and [latex]\phi[/latex] and use calculus to find the least squares estimators of μ and [latex]\phi[/latex]. This will result in a slight difference between the usual pattern of using the sample mean [latex]\bar x[/latex] to estimate the population mean μ. Once these least squares estimators have been determined, the population variance of the white noise [latex]\sigma _ Z ^ {\, 2}[/latex] will be estimated.

The sum of squared errors is

$$\begin{array}{l}S=\sum _{t=2}^n{Z}_t^2=\sum _{t=2}^n{[X_t-\mu -\varphi (X_{t-1}-\mu )]}^2.\end{array}$$

The partial derivatives of S with respect to μ and [latex]\phi[/latex] are

$$\begin{array}{l}\frac{\partial S}{\partial \mu }=\sum _{t=2}^{n}2[X_t-\mu -\varphi (X_{t-1}-\mu )](-1+\varphi )\end{array}$$

and

$$\begin{array}{l}\frac{\partial S}{\partial \varphi }=\sum _{t=2}^n-2[X_t-\mu -\varphi (X_{t-1}-\mu )](X_{t-1}-\mu ).\end{array}$$

Equating the first of the partial derivatives to zero yields

$$\begin{array}{l}\sum _{t=2}^n[X_t-\mu -\varphi (X_{t-1}-\mu )]=0\end{array}$$

or

$$\begin{array}{l}\sum _{t=2}^n{X}_t-\varphi \sum _{t=2}^n{X}_{t-1}-(n-1)\mu (1-\varphi )=0\end{array}$$

or

$$\begin{array}{l}{\overline{X}}_2-\varphi {\overline{X}}_1-\mu (1-\varphi )=0\end{array}$$

or

$$\begin{array}{l}\hat{\mu }=\frac{{\overline{X}}_2-\hat{\varphi }{\overline{X}}_1}{1-\hat{\varphi }},\end{array}$$

where

$$\begin{array}{l}{\overline{X}}_1=\frac{1}{n-1}\sum _{t=1}^{n-1}{X}_t\text{and}{\overline{X}}_2=\frac{1}{n-1}\sum _{t=2}^n{X}_t.\end{array}$$

Equating the second of the partial derivatives to zero yields

$$\begin{array}{l}\sum _{t=2}^n(X_t-\mu )(X_{t-1}-\mu )-\varphi \sum _{t=2}^n{(X_{t-1}-\mu )}^2=0\end{array}$$

or

$$\begin{array}{l}\hat{\varphi }=\frac{\sum _{t=2}^n(X_t-\hat{\mu })(X_{t-1}-\hat{\mu })}{\sum _{t=2}^n{(X_{t-1}-\hat{\mu })}^2}.\end{array}$$

So the ordinary least squares estimators for μ and [latex]\phi[/latex] can be determined by numerically solving the simultaneous equations

$$\begin{array}{l}\hat{\mu }=\frac{{\overline{X}}_2-\hat{\varphi }{\overline{X}}_1}{1-\hat{\varphi }}\text{and}\hat{\varphi }=\frac{\sum _{t=2}^n(X_t-\hat{\mu })(X_{t-1}-\hat{\mu })}{\sum _{t=2}^n{(X_{t-1}-\hat{\mu })}^2}\end{array}$$

for [latex]\hat \mu[/latex] and [latex]\hat \phi[/latex].

The last parameter to estimate is [latex]\sigma _ Z ^ {\, 2}[/latex]. Since

$$\begin{array}{l}\gamma (0)=\left(\frac{1}{1-{\varphi }^2}\right){\sigma }_Z^2\end{array}$$

for an AR(1) time series model, the population variance of the white noise can be expressed as

$$\begin{array}{l}{\sigma }_Z^2=(1-{\varphi }^2)\gamma (0).\end{array}$$

Replacing [latex]\phi[/latex] by the estimator r₁ because [latex]\rho(1) = \phi[/latex], and replacing [latex]\gamma(0) = V \left[ X_t \right][/latex] by the estimator [latex]c_0 = \frac{1}{n} \sum_{t\,=\,1}^n \left( X_t - \bar X \right) ^ 2[/latex] gives the point estimator

$$\begin{array}{l}{\hat{\sigma }}_Z^2=(1-r_1^2)c_0,\end{array}$$

which matches the method of moments estimator from Theorem 9.4. This derivation constitutes a proof of the following result.

We now apply these techniques to the beaver temperature data set from Example 9.3.

Since

$$\begin{array}{l}{\overline{X}}_1=\frac{1}{n-1}\sum _{t=1}^{n-1}{X}_t\text{and}{\overline{X}}_2=\frac{1}{n-1}\sum _{t=2}^n{X}_t\end{array}$$

contain the[latex]n - 2[/latex] common values [latex]X_2, \, X_3, \, \ldots, \, X_{n - 1}[/latex], one approximation that can be applied to the least squares estimates is to assume that [latex]\bar X _ 1 \cong \bar X _ 2 \cong \bar X[/latex] for large values of n, which allows for closed-form approximate least squares estimators:

$$\begin{array}{l}\hat{\mu }=\overline{X}\text{and}\hat{\varphi }=\frac{\sum _{t=2}^n(X_t-\overline{X})(X_{t-1}-\overline{X})}{\sum _{t=2}^n{(X_{t-1}-\overline{X})}^2}.\end{array}$$

As a secondary additional approximation, the denominator of [latex]\hat \phi[/latex] with the first approximation in place,

$$\begin{array}{l}\sum _{t=2}^n{(X_{t-1}-\overline{X})}^2,\end{array}$$

is approximately equal to

$$\begin{array}{l}\sum _{t=1}^n{(X_t-\overline{X})}^2\end{array}$$

for large values of n. With this additional assumption, the least squares estimate for [latex]\phi[/latex] reduces to the approximate least squares estimate

$$\begin{array}{l}\hat{\varphi }=\frac{c_1}{c_0}=r_1,\end{array}$$

which is the method of moments estimator of [latex]\phi[/latex] because [latex]\rho(1) = \phi[/latex] for an AR(1) model. With both approximations in place, the least squares estimators exactly match the method of moments estimators. This is why the estimates from the two techniques are so close.

Approach 3: Maximum likelihood estimation. The likelihood function is the joint probability density function of the observed values in the time series [latex]x_1, \, x_2, \, \ldots, \, x_n[/latex] in a shifted AR(1) model is

$$\begin{array}{l}L(\mu ,\varphi ,{\sigma }_Z^2)=f(x_1,x_2,\dots ,x_n),\end{array}$$

where the [latex]x_1, \, x_2, \, \ldots, \, x_n[/latex] arguments on L and the μ, [latex]\phi[/latex], and [latex]\sigma _Z ^ {\, 2}[/latex] arguments on f have been dropped for brevity. It is not possible to simply multiply the marginal probability density functions because the values in the AR(1) time series model are correlated. In order to use maximum likelihood estimation, we make the additional assumption that the white noise terms [latex]Z_1, \, Z_2, \, \ldots, \, Z_n[/latex] are in fact Gaussian white noise terms:

$$\begin{array}{l}{f}_{Z_t}(z_t)=\frac{1}{\sqrt{2\pi {\sigma }_Z^2}}{e}^{-z_t^2/\left(2{\sigma }_Z^2\right)}-\mathrm{\infty }<z_t<\mathrm{\infty }\end{array}$$

for [latex]t = 1, \, 2, \, \ldots, \, n[/latex], which is the probability density function of a [latex]N\left( 0, \, \sigma _Z ^ {\, 2} \right)[/latex] random variable. Ignoring Z₁ temporarily, the joint probability density function of the mutually independent white noise random variables [latex]Z_2, \, Z_3, \, \ldots, \, Z_n[/latex] is

$$\begin{array}{l}{f}_{Z_2,Z_3,\dots ,Z_n}(z_2,z_3,\dots ,z_n)={\left(2\pi {\sigma }_Z^2\right)}^{-(n-1)/2}{e}^{-\sum _{t=2}^n{z}_t^2/\left(2{\sigma }_Z^2\right)}\end{array}$$

for [latex]\left( z_2, \, z_3, \, \ldots , \, z_n \right) \in {\cal R} ^ {n - 1}[/latex]. The shifted AR(1) model

$$\begin{array}{l}{X}_t-\mu =\varphi (X_{t-1}-\mu )+Z_t\end{array}$$

applies for all values of t, so

$$\begin{array}{ll}{X}_2-\mu & =\varphi (X_1-\mu )+Z_2\\ X_3-\mu & =\varphi (X_2-\mu )+Z_3\\ & ⋮\\ X_n-\mu & =\varphi (X_{n-1}-\mu )+Z_n.\end{array}$$

Solving these equations for [latex]X_2, \, X_3, \, \ldots, \, X_n[/latex], consider the transformation of the [latex]Z_2, \, Z_3, \, \ldots , \, Z_n[/latex] values

$$\begin{array}{ll}{X}_2& =\mu +\varphi (X_1-\mu )+Z_2\\ X_3& =\mu +\varphi (X_2-\mu )+Z_3\\ & ⋮\\ X_n& =\mu +\varphi (X_{n-1}-\mu )+Z_n\end{array}$$

conditioned on [latex]X_1 = x_1[/latex], which is a one-to-one transformation from [latex]{\cal R}^{n - 1}[/latex] to [latex]{\cal R}^{n - 1}[/latex] with inverse transformation

$$\begin{array}{ll}{Z}_2& =X_2-\mu -\varphi (X_1-\mu )\\ Z_3& =X_3-\mu -\varphi (X_2-\mu )\\ & ⋮\\ Z_n& =X_n-\mu -\varphi (X_{n-1}-\mu )\end{array}$$

and Jacobian

$$\begin{array}{l}J=\left|\begin{array}{cccccc}1& -\varphi & 0& \cdots & 0& 0\\ 0& 1& -\varphi & \cdots & 0& 0\\ 0& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & 1& -\varphi \\ 0& 0& 0& \cdots & 0& 1\end{array}\right|=1.\end{array}$$

By the transformation technique, the joint probability density function of [latex]X_2, \, X_3, \, \ldots, \, X_n[/latex] conditioned on [latex]X_1 = x_1[/latex] is

$$\begin{array}{l}{f}_{X_2,X_3,\dots ,X_n|X_1}(x_2,x_3,\dots ,x_n|X_1=x_1)={\left(2\pi {\sigma }_Z^2\right)}^{-(n-1)/2}{e}^{-\sum _{t=2}^n{[x_t-\mu -\varphi (x_{t-1}-\mu )]}^2/\left(2{\sigma }_Z^2\right)}\end{array}$$

for [latex]\left( x_2, \, x_3, \, \ldots , \, x_n \right) \in {\cal R} ^ {n - 1}[/latex] and [latex]x_1 \in {\cal R}[/latex]. The final step in the derivation of the likelihood function involves determining the marginal distribution of X₁. Since

$$\begin{array}{l}{X}_1\sim N(\mu ,\frac{{\sigma }_Z^2}{1-{\varphi }^2}),\end{array}$$

the probability density function of X₁ is

$$\begin{array}{l}{f}_{X_1}(x_1)=\sqrt{\frac{1-{\varphi }^2}{2\pi {\sigma }_Z^2}}{e}^{-(1-{\varphi }^2)(x_1-\mu {)}^2/\left(2{\sigma }_Z^2\right)}-\mathrm{\infty }<x_1<\mathrm{\infty }.\end{array}$$

The joint probability density function of [latex]X_1, \, X_2, \, \ldots, \, X_n[/latex] is the product of the conditional probability density function and the marginal probability density function:

$$\begin{array}{l}{f}_{X_1,X_2,\dots ,X_n}(x_1,x_2,\dots ,x_n)=f_{X_2,X_3,\dots ,X_n|X_1}(x_2,x_3,\dots ,x_n|X_1=x_1)f_{X_1}(x_1)\end{array}$$

for [latex]\left( x_1, \, x_2, \, \ldots, \, x_n \right) \in {\cal R} ^ n[/latex]. So the likelihood function is

$$\begin{array}{l}L(\mu ,\varphi ,{\sigma }_Z^2)={\left(2\pi {\sigma }_Z^2\right)}^{-n/2}\sqrt{1-{\varphi }^2}{e}^{-S(\mu ,\varphi )/\left(2{\sigma }_Z^2\right)},\end{array}$$

where the unconditional sum of squares is

$$\begin{array}{l}S(\mu ,\varphi )=(1-{\varphi }^2){(x_1-\mu )}^2+\sum _{t=2}^n{[(x_t-\mu )-\varphi (x_{t-1}-\mu )]}^2.\end{array}$$

The associated log likelihood function is

$$\begin{array}{l}\ln L(\mu ,\varphi ,{\sigma }_Z^2)=-\frac{n}{2}\ln \left(2\pi {\sigma }_Z^2\right)+\frac{1}{2}\ln (1-{\varphi }^2)-\frac{S(\mu ,\varphi )}{2{\sigma }_Z^2}.\end{array}$$

The maximum likelihood estimators [latex]\hat \mu[/latex], [latex]\hat \phi[/latex], and [latex]\hat \sigma _ Z ^ {\, 2}[/latex] satisfy

$$\begin{array}{ll}\frac{\partial \ln L(\mu ,\varphi ,{\sigma }_Z^2)}{\partial \mu }& =\frac{(1-{\varphi }^2)(x_1-\mu )+(1-\varphi )\sum _{t=2}^n[(x_t-\mu )-\varphi (x_{t-1}-\mu )]}{{\sigma }_Z^2}=0\\ \frac{\partial \ln L(\mu ,\varphi ,{\sigma }_Z^2)}{\partial \varphi }& =-\frac{\varphi }{1-{\varphi }^2}+\frac{\varphi (x_1-\mu {)}^2+\sum _{t=2}^n[(x_t-\mu )-\varphi (x_{t-1}-\mu )](x_{t-1}-\mu )}{{\sigma }_Z^2}=0\\ \frac{\partial \ln L(\mu ,\varphi ,{\sigma }_Z^2)}{\partial {\sigma }_Z^2}& =-\frac{n}{2{\sigma }_Z^2}+\frac{S(\mu ,\varphi )}{2{\sigma }_Z^4}=0.\end{array}$$

Although the third equation satisfies

$$\begin{array}{l}{\hat{\sigma }}_Z^2=\frac{S(\hat{\mu },\hat{\varphi })}{n},\end{array}$$

numerical methods are required to solve the equations.

Maximum likelihood estimation will be illustrated in the next example.

Table 9.2 summarizes the point estimators that have been calculated in the previous three examples for the [latex]n = 62[/latex] beaver temperatures. The point estimators associated with the three methods are quite close for this particular time series.

The R function ar fits autoregressive models. The parameter estimates from the three previous examples could have been calculated with the following four R statements.

Table 9.3 contains the point estimates returned by the ar function. The tiny differences between some of the entries in Tables 9.2 and 9.3 might be due to slightly different approximations and/or roundoff in the optimization routines.

We have now derived and illustrated the three point estimation techniques, the method of moments, least squares, and maximum likelihood estimation, for the parameters in an AR(1) model from a realization of a time series consisting of n observations. Which of these techniques provides the best point estimators? This is not an easy question to answer because there are a large number of factors, such as the sample size n, the values of the parameters in the model, and the fact that there are three parameters to estimate. There will not necessarily be one universal answer to the question. We do a focused evaluation on the point estimator for [latex]\phi[/latex] because it typically differs for the three methods of point estimation. The mean square error associated with the point estimate for [latex]\phi[/latex] is

$$\begin{array}{l}E\left[{(\hat{\varphi }-\varphi )}^2\right].\end{array}$$

The following R code conducts a Monte Carlo simulation experiment which estimates the mean square error of the three point estimators for [latex]\phi[/latex] for 40,000 replications. We selected the time series model with [latex]\mu = 38[/latex], [latex]\phi = 0.8[/latex], [latex]\sigma _ Z = 0.13[/latex], and [latex]n = 62[/latex], which are parameters that are near the estimated parameters in the last three examples involving the time series of beaver temperatures.

After a call to set.seed(4) to establish the random number stream, three runs of this simulation yielded the following estimated mean squared error values:

$$\begin{array}{l}\text{Method of moments}:0.0135\text{Least squares}:0.0117\text{Maximum likelihood}:0.0113.\end{array}$$

Furthermore, confidence intervals for the three methods do not overlap. Since small values of the mean square error are preferred, we conclude that the maximum likelihood estimator is the preferred estimator for these parameter settings, followed by the least squares estimator, followed by the method of moments estimator in a distant third place.

The focus on estimation thus far has been on point estimation techniques. We also want to report some indication of the precision associated with these point estimators. In the previous example, the sampling distributions of [latex]\hat \mu[/latex], [latex]\hat \phi[/latex], and [latex]\hat \sigma _ Z ^ {\, 2}[/latex] in the AR(1) model are too complicated to derive analytically. As an illustration of a confidence interval for one of the parameters, we use the asymptotic normality of the maximum likelihood estimator of [latex]\phi[/latex] in the result:

$$\begin{array}{l}\hat{\varphi }\stackrel{D}{\to }N(\varphi ,\frac{1-{\varphi }^2}{n}).\end{array}$$

This result leads to an asymptotically exact two-sided [latex]{100(1 - \alpha)}\%[/latex] confidence interval for [latex]\phi[/latex].

This asymptotically exact confidence interval will now be illustrated with the time series of active beaver temperatures from the three previous examples.

Model Assessment

Now that techniques for point and interval estimates for the parameters in the AR(1) model have been established, we are interested in assessing the adequacy of the AR(1) time series model. This will involve an analysis of the residuals. Recall from Section 8.2.3 that the residuals are defined by

$$\begin{array}{l}\left[\text{residual}\right]=\left[\text{observed value}\right]-\left[\text{predicted value}\right]\end{array}$$

or

$$\begin{array}{l}{\hat{Z}}_t=X_t-{\hat{X}}_t.\end{array}$$

Since [latex]\hat{X} _ {t}[/latex] is the one-step-ahead forecast from the time origin [latex]t - 1[/latex], this is more clearly written as

$$\begin{array}{l}{\hat{Z}}_t=X_t-{\hat{X}}_{t-1}(1).\end{array}$$

Therefore, for the time series [latex]x_1, \, x_2, \, \ldots, \, x_n[/latex] and the fitted AR(1) model with parameter estimates [latex]\hat \mu[/latex] and [latex]\hat \phi[/latex], the residual at time t is

$$\begin{array}{l}{\hat{Z}}_t=x_t-[\hat{\mu }+\hat{\varphi }(x_{t-1}-\hat{\mu })]\end{array}$$

for [latex]t = 2, \, 3, \, \ldots, \, n[/latex] via Example 8.12. The next example shows the steps associated with assessing the adequacy of the AR(1) model for the active beaver temperature time series.

We have seen a number of indicators that the AR(1) time series model is an adequate model for the active beaver temperatures. But how do we know that there is not a better model with more terms lurking below the surface that might provide a better fit? The next subsection considers the process of model selection.

Model Selection

One way of eliminating the possibility of a better time series model is to overfit the tentative AR(1) time series model with ARMA(p, q) models of higher order. We have not yet surveyed the techniques for estimating the parameters in these models with additional terms, so for now we will let the arima function in R estimate their parameters and compare them via their AIC (Akaike’s Information Criterion) statistics. The AIC statistic was introduced in Section 8.2.4.

In some applications, just describing the time series model for the beaver temperatures in the active state with the fitted AR(1) model is adequate. In other applications, simulating the values in the fitted AR(1) model is the goal. But in many application areas, particularly economics, there is often an interest in forecasting future values of a time series from a realization. In our setting, we might be interested in this particular beaver’s future temperature based on the [latex]n = 62[/latex] temperature values collected. The next subsection considers forecasting for the AR(1) model.

Forecasting

We now pivot to the development of a procedure to forecast future values of a time series that is governed by an AR(1) model. To review the notation for forecasting, the observed time series values are [latex]x_1, \, x_2, \, \ldots, \, x_n[/latex]. The forecast is being made at time [latex]t = n[/latex]. The random future value of the time series that is h time units in the future is denoted by [latex]X_{n + h}[/latex]. The associated forecasted value is denoted by [latex]\hat{X}_{n + h}[/latex], and is the conditional expected value

$$\begin{array}{l}{\hat{X}}_{n+h}=E[X_{n+h}|X_1=x_1,X_2=x_2,\dots ,X_n=x_n].\end{array}$$

We would like to calculate this forecasted value and an associated prediction interval for the AR(1) model. As in Section 8.2.2, we assume that all parameters are known in the derivations that follow.

Recall from Example 8.12 that the forecasted value for one time unit into the future for a shifted AR(1) model is

$$\begin{array}{l}{\hat{X}}_{n+1}=\mu +\varphi (x_n-\mu ).\end{array}$$

We would like to generalize this so as to find the forecasted value h time units into the future. In other words, we want to find [latex]\hat{X} _ {n + h}[/latex]. The shifted AR(1) model is

$$\begin{array}{l}{X}_t-\mu =\varphi (X_{t-1}-\mu )+Z_t.\end{array}$$

Replacing t by [latex]n + h[/latex], which is the time value of interest, gives

$$\begin{array}{l}{X}_{n+h}-\mu =\varphi (X_{n+h-1}-\mu )+Z_{n+h}.\end{array}$$

Taking the conditional expected value of each side of this equation results in

$$\begin{array}{l}{\hat{X}}_{n+h}=\mu +\varphi ({\hat{X}}_{n+h-1}-\mu ).\end{array}$$

Iterating on this equation for time values that are sequentially one time unit closer to the present time [latex]t = n[/latex] yields

$$\begin{array}{ll}{\hat{X}}_{n+h}& =\mu +\varphi ({\hat{X}}_{n+h-1}-\mu )\\ & =\mu +\varphi [\mu +\varphi ({\hat{X}}_{n+h-2}-\mu )-\mu ]\\ & =\mu +{\varphi }^2({\hat{X}}_{n+h-2}-\mu )\\ & ⋮\\ & =\mu +{\varphi }^{h-1}({\hat{X}}_{n+1}-\mu )\\ & =\mu +{\varphi }^{h-1}[\mu +\varphi (x_n-\mu )-\mu ]\\ & =\mu +{\varphi }^h(x_n-\mu ).\end{array}$$

Notice that the forecasted value is a function of x_n, but not a function of [latex]x_1, \, x_2, \, \ldots, \, x_{n - 1}[/latex]. This is a sensible forecast in the sense that for a long time horizon h into the future and a stationary shifted AR(1) model with [latex]-1 < \phi < 1[/latex],

$$\begin{array}{l}\underset{h\to \mathrm{\infty }}{lim}{\hat{X}}_{n+h}=\mu .\end{array}$$

If you were asked to forecast your temperature one year from now, you would probably say 98.6° Fahrenheit (or whatever your average temperature might be), regardless of whether you are healthy or have a fever right now. Long-term forecasts for stationary time series models always tend to the population mean.

As is typically the case in statistics, we would like to pair our point estimator [latex]\hat{X}_{n + h}[/latex] with an interval estimator, which is a prediction interval in this setting. The prediction interval gives us an indication of the precision of the forecast. In order to derive an exact two-sided [latex]{100(1 - \alpha)}\%[/latex] prediction interval for [latex]X _ {n + h}[/latex], it is helpful to write the shifted AR(1) model as a shifted MA(∞) model. Using successive substitutions, each one time unit prior to the previous substitution,

$$\begin{array}{ll}{X}_t-\mu & =\varphi (X_{t-1}-\mu )+Z_t\\ & =\varphi [\varphi (X_{t-2}-\mu )+Z_{t-1}]+Z_t\\ & ={\varphi }^2(X_{t-2}-\mu )+Z_t+\varphi Z_{t-1}\\ & ={\varphi }^2[(\varphi X_{t-3}-\mu )+Z_{t-2}]+Z_t+\varphi Z_{t-1}\\ & ={\varphi }^3(X_{t-3}-\mu )+Z_t+\varphi Z_{t-1}+{\varphi }^2{Z}_{t-2}\\ & ⋮\end{array}$$

For [latex]-1 < \phi < 1[/latex] corresponding to a stationary shifted AR(1) model, the limiting expression for X_t is

$$\begin{array}{l}{X}_t=\mu +Z_t+\varphi Z_{t-1}+{\varphi }^2{Z}_{t-2}+\cdots ,\end{array}$$

which is a shifted MA(∞) model. Replacing t with [latex]n + h[/latex] results in

$$\begin{array}{l}{X}_{n+h}=\mu +Z_{n+h}+\varphi Z_{n+h-1}+{\varphi }^2{Z}_{n+h-2}+\cdots ,\end{array}$$

Taking the conditional variance of both sides of this equation yields

$$\begin{array}{ll}V[X_{n+h}|& X_1=x_1,X_2=x_2,\dots ,X_n=x_n]\\ & =V[\mu +Z_{n+h}+\varphi Z_{n+h-1}+{\varphi }^2{Z}_{n+h-2}+\cdots |X_1=x_1,X_2=x_2,\dots ,X_n=x_n]\\ & ={\sigma }_Z^2+{\varphi }^2{\sigma }_Z^2+{\varphi }^4{\sigma }_Z^2+\cdots +{\varphi }^{2h-2}{\sigma }_Z^2\\ & =(1+{\varphi }^2+{\varphi }^4+\cdots +{\varphi }^{2h-2}){\sigma }_Z^2\\ & =\frac{1-{\varphi }^{2h}}{1-{\varphi }^2}{\sigma }_Z^2\end{array}$$

because the error terms at time n and prior are observed and can therefore be treated as constants. Assuming Gaussian white noise terms, an exact two-sided [latex]{100(1 - \alpha)}\%[/latex] prediction interval for [latex]X _ {n + h}[/latex] is

$$\begin{array}{l}{\hat{X}}_{n+h}-z_{\alpha /2}\sqrt{\frac{1-{\varphi }^{2h}}{1-{\varphi }^2}}{\sigma }_Z<X_{n+h}<{\hat{X}}_{n+h}+z_{\alpha /2}\sqrt{\frac{1-{\varphi }^{2h}}{1-{\varphi }^2}}{\sigma }_Z.\end{array}$$

In most practical problems, the parameters in this prediction interval will be estimated from data, which results in the following approximate two-sided [latex]{100(1 - \alpha)}\%[/latex] prediction interval.

This subsection has introduced the AR(1) time series model. The key results for an AR(1) model are listed below.

• The standard AR(1) model can be written algebraically and with the backshift operator B as

  $$\begin{array}{l}{X}_t=\varphi X_{t-1}+Z_t\text{and}(1-\varphi B)X_t=Z_t,\end{array}$$

  where [latex]Z_t \sim WN \left( 0, \, \sigma _Z ^ {\, 2} \right)[/latex] and [latex]\sigma _ Z ^ {\, 2} > 0[/latex].

• The shifted AR(1) model can be written algebraically and with the backshift operator B as

  $$\begin{array}{l}{X}_t-\mu =\varphi (X_{t-1}-\mu )+Z_t\text{and}(1-\varphi B)(X_t-\mu )=Z_t.\end{array}$$

• The AR(1) model is always invertible; the AR(1) model is stationary for [latex]-1 < \phi < 1[/latex].
• The stationary shifted AR(1) model can be written as an MA(∞) model for [latex]-1 < \phi < 1[/latex] as

  $$\begin{array}{l}{X}_t=\mu +Z_t+\varphi Z_{t-1}+{\varphi }^2{Z}_{t-2}+\cdots .\end{array}$$

• The AR(1) population autocorrelation function is [latex]\rho(k) = \phi ^ k[/latex] for [latex]-1 < \phi < 1[/latex] and [latex]k = 1, \, 2, \, \ldots[/latex].
• The AR(1) population partial autocorrelation function at lag one is [latex]\rho ^ * (1) = \phi[/latex] for [latex]-1 < \phi < 1[/latex] and [latex]\rho ^ * (k) = 0[/latex] for [latex]k = 2, \, 3, \, \ldots[/latex].
• The three parameters in the shifted AR(1) model, μ, [latex]\phi[/latex], and [latex]\sigma _ Z ^ {\, 2}[/latex], can be estimated from a realization of a time series [latex]x_1, \, x_2, \, \ldots, \, x_n[/latex] by the method of moments, least squares, and maximum likelihood. The point estimators for μ, [latex]\phi[/latex], and [latex]\sigma _ Z ^ {\, 2}[/latex] are denoted by [latex]\hat \mu[/latex], [latex]\hat \phi[/latex], and [latex]\hat{\sigma} _ Z ^ {\, 2}[/latex], and are typically paired with asymptotically exact two-sided [latex]{100(1 - \alpha)}\%[/latex] confidence intervals.
• The forecast value [latex]\hat{X} _ {n + h} = \hat \mu + \hat \phi ^ h \left( x_n - \hat{\mu} \right)[/latex] for an AR(1) model approaches [latex]\hat \mu = \bar x[/latex] as [latex]h \rightarrow \infty[/latex]. The associated prediction intervals have widths that increase as h increases and approach a limit as the time horizon [latex]h \rightarrow \infty[/latex].

If the time series of interest is the daily high temperatures in July in Tuscaloosa, then an AR(1) model would be appropriate if tomorrow’s daily high temperature (X_t) can be modeled as a linear function of

• today’s high temperature ([latex]X_{t-1}[/latex]), and
• a random shock (Z_t).

But what if weather had more of a memory than just one day? What if tomorrow’s daily high temperature (X_t) is better modeled as a linear function of

• today’s high temperature ([latex]X_{t-1}[/latex]),
• yesterday’s high temperature ([latex]X_{t-2}[/latex]), and
• a random shock (Z_t).

This is an example of the thinking that lies behind the AR(2) model, which is introduced in the next section.

9.1.2 The AR(2) Model

The second-order autoregressive model, denoted by AR(2), can be used for modeling a stationary time series in instances in which the current value of the time series is a linear combination of the two previous values plus a random shock. The mathematics associated with the AR(2) model is somewhat more difficult than that associated with the AR(1) model.

There are three parameters that define an AR(2) model: the real-valued coefficients [latex]\phi_1[/latex] and [latex]\phi_2[/latex], and the population variance of the white noise [latex]\sigma _ Z ^ {\, 2}[/latex]. The AR(2) model can be written more compactly in terms of the backshift operator B as

$$\begin{array}{l}\varphi (B)X_t=Z_t,\end{array}$$

where [latex]\phi(B)[/latex] is the second-order polynomial

$$\begin{array}{l}\varphi (B)=1-{\varphi }_{1}B-{\varphi }_2{B}^2.\end{array}$$

The AR(2) model has the form of a multiple linear regression model with two independent variables and no intercept term. The current value X_t is modeled as a linear combination of the two previous values of the time series, [latex]X_{t-1}[/latex] and [latex]X_{t - 2}[/latex], plus a white noise term. The parameters [latex]\phi_1[/latex] and [latex]\phi_2[/latex] control the inclination of the regression plane in three-dimensional space. The parameter [latex]\sigma _ Z ^ {\, 2}[/latex] reflects the magnitude of the dispersion of the time series values from the regression plane.

To illustrate the thinking behind the AR(2) model in a specific context, let X_t represent the annual return of a particular stock market index in year t. The AR(2) model indicates that the annual return in year t equals [latex]\phi_1[/latex] multiplied by the previous year’s annual return ([latex]\phi_1 X_{t - 1}[/latex]), plus [latex]\phi_2[/latex] multiplied by the annual return two years prior ([latex]\phi_2 X_{t - 2}[/latex]), plus the year t random white noise term Z_t.

Stationarity

Theorem 8.3 indicates that all AR(2) models are invertible, but are stationary when the roots of

$$\begin{array}{l}\varphi (B)=1-{\varphi }_{1}B-{\varphi }_2{B}^2\end{array}$$

lie outside of the unit circle in the complex plane. Let B₁ and B₂ denote these two roots. Using the quadratic equation, the two roots are

$$\begin{array}{l}{B}_1=\frac{{\varphi }_1-\sqrt{{\varphi }_1^2+4{\varphi }_2}}{-2{\varphi }_2}\text{and}{B}_2=\frac{{\varphi }_1+\sqrt{{\varphi }_1^2+4{\varphi }_2}}{-2{\varphi }_2}.\end{array}$$

Since [latex]\phi(B_1) = \phi(B_2) = 0[/latex], the quadratic function [latex]\phi(B)[/latex] can also be written in factored form as

$$\begin{array}{l}\varphi (B)=(1-B_1^{-1}B)(1-B_2^{-1}B).\end{array}$$

Equating the two versions of [latex]\phi(B)[/latex] above and matching coefficients results in

$$\begin{array}{l}{\varphi }_1=B_1^{-1}+B_2^{-1}\text{and}{\varphi }_2=-(B_1{B}_2{)}^{-1}.\end{array}$$

These two equations define the mapping from the complex plane, which contains the roots B₁ and B₂, to the plane that contains the AR(2) parameters [latex]\phi_1[/latex] and [latex]\phi_2[/latex]. To find the stationary region, we must find the mapping of the part of the complex plane outside of the unit circle to the [latex](\phi_1, \, \phi_2)[/latex] plane. The mapping yields a triangular-shaped stationary region, as specified in the following result.

Population Autocorrelation Function

Now that the stationary region for an AR(2) time series model has been established, we turn to the derivation of the population autocorrelation function. Assuming that the parameters [latex]\phi_1[/latex] and [latex]\phi_2[/latex] fall in the stationary region, the AR(2) model

$$\begin{array}{l}{X}_t={\varphi }_1{X}_{t-1}+{\varphi }_2{X}_{t-2}+Z_t\end{array}$$

can be multiplied by [latex]X_{t - k}[/latex] to give

$$\begin{array}{l}{X}_t{X}_{t-k}={\varphi }_1{X}_{t-1}{X}_{t-k}+{\varphi }_2{X}_{t-2}{X}_{t-k}+Z_t{X}_{t-k}.\end{array}$$

Taking the expected value of both sides of this equation results in the recursive equation

$$\begin{array}{l}\gamma (k)={\varphi }_1\gamma (k-1)+{\varphi }_2\gamma (k-2)\end{array}$$

for [latex]k = 1, \, 2, \, \ldots[/latex] because Z_t has expected value zero and is independent of [latex]X_{t - k}[/latex]. Dividing both sides of this equation by [latex]\gamma(0) = V \left[ X_t \right][/latex] gives the recursive equation

$$\begin{array}{l}\rho (k)={\varphi }_1\rho (k-1)+{\varphi }_2\rho (k-2)\end{array}$$

for [latex]k = 1, \, 2, \, \ldots[/latex]. These linear equations, whether written in terms of [latex]\gamma(k)[/latex] or [latex]\rho(k)[/latex], are known in time series analysis as the Yule–Walker equations after British statisticians George Udny Yule and Sir Gilbert Walker. Once the first two values of [latex]\gamma(k)[/latex] or [latex]\rho(k)[/latex] are known, these recursive equations can be used to calculate subsequent values. The next two paragraphs focus on determining the first two values of [latex]\gamma(k)[/latex] and [latex]\rho(k)[/latex], respectively.

For a stationary AR(2) time series model, we derive expressions for [latex]\gamma(0)[/latex] and [latex]\gamma(1)[/latex]. The AR(2) model is

$$\begin{array}{l}{X}_t={\varphi }_1{X}_{t-1}+{\varphi }_2{X}_{t-2}+Z_t.\end{array}$$

Squaring both sides of this equation and taking the expected value of both sides gives

$$\begin{array}{l}\gamma (0)={\varphi }_1^2\gamma (0)+{\varphi }_2^2\gamma (0)+{\sigma }_Z^2+2{\varphi }_1{\varphi }_2\gamma (1).\end{array}$$

Using the symmetry of the population autocovariance function, the Yule–Walker equation with [latex]k = 1[/latex] is

$$\begin{array}{l}\gamma (1)={\varphi }_1\gamma (0)+{\varphi }_2\gamma (1)\Rightarrow \gamma (1)=\frac{{\varphi }_1\gamma (0)}{1-{\varphi }_2}.\end{array}$$

Replacing this expression for [latex]\gamma(1)[/latex] in the previous equation gives

$$\begin{array}{l}\gamma (0)={\varphi }_1^2\gamma (0)+{\varphi }_2^2\gamma (0)+{\sigma }_Z^2+2{\varphi }_1{\varphi }_2\frac{{\varphi }_1\gamma (0)}{1-{\varphi }_2}.\end{array}$$

Moving all terms involving [latex]\gamma(0)[/latex] to the left-hand side of this equation gives

$$\begin{array}{l}\gamma (0)[1-{\varphi }_1^2-{\varphi }_2^2-2{\varphi }_1^2{\varphi }_2\frac{1}{1-{\varphi }_2}]={\sigma }_Z^2.\end{array}$$

Solving this equation for [latex]\gamma(0)[/latex],

$$\begin{array}{ll}\gamma (0)& =\frac{(1-{\varphi }_2){\sigma }_Z^2}{(1-{\varphi }_2)-(1-{\varphi }_2){\varphi }_1^2-(1-{\varphi }_2){\varphi }_2^2-2{\varphi }_1^2{\varphi }_2}\\ & =\frac{(1-{\varphi }_2){\sigma }_Z^2}{1-{\varphi }_2-{\varphi }_1^2+{\varphi }_1^2{\varphi }_2-{\varphi }_2^2+{\varphi }_2^3-2{\varphi }_1^2{\varphi }_2}\\ & =\frac{(1-{\varphi }_2){\sigma }_Z^2}{(1+{\varphi }_2)(1+{\varphi }_1-{\varphi }_2)(1-{\varphi }_1-{\varphi }_2)}.\end{array}$$

An expression for [latex]\gamma(1)[/latex] is

$$\begin{array}{l}\gamma (1)=\frac{{\varphi }_1\gamma (0)}{1-{\varphi }_2}=\frac{{\varphi }_1{\sigma }_Z^2}{(1+{\varphi }_2)(1+{\varphi }_1-{\varphi }_2)(1-{\varphi }_1-{\varphi }_2)}.\end{array}$$

These two values can be used as arguments in the Yule–Walker equations to obtain subsequent values for [latex]\gamma(k)[/latex].

We now turn to the problem of finding [latex]\rho(1)[/latex] and [latex]\rho(2)[/latex]. The first two Yule–Walker equations in terms of [latex]\rho(k)[/latex] are

$$\begin{array}{ll}\rho (1)& ={\varphi }_1\rho (0)+{\varphi }_2\rho (-1)\\ \rho (2)& ={\varphi }_1\rho (1)+{\varphi }_2\rho (0).\end{array}$$

Since [latex]\rho(0) = 1[/latex] and [latex]\rho(-k) = \rho(k)[/latex] via Theorem 7.1, these equations reduce to

$$\begin{array}{ll}\rho (1)& ={\varphi }_1+{\varphi }_2\rho (1)\\ \rho (2)& ={\varphi }_1\rho (1)+{\varphi }_2,\end{array}$$

which are easily solved for [latex]\rho(1)[/latex] and [latex]\rho(2)[/latex]:

$$\begin{array}{l}\rho (1)=\frac{{\varphi }_1}{1-{\varphi }_2}\text{and}\rho (2)=\frac{{\varphi }_1^2}{1-{\varphi }_2}+{\varphi }_2.\end{array}$$

A general formula for [latex]\rho(k)[/latex] exists, but it can involve complex numbers and is unwieldy. An exercise concerning its calculation is given at the end of the chapter. These results are summarized in the following theorem.

We now focus in on the values of [latex]\rho(1)[/latex] and [latex]\rho(2)[/latex]. We can solve for [latex]\phi_1[/latex] and [latex]\phi_2[/latex] in terms of [latex]\rho(1)[/latex] and [latex]\rho(2)[/latex] as

$$\begin{array}{l}{\varphi }_1=\frac{\rho (1)(1-\rho (2))}{1-\rho (1{)}^2}\text{and}{\varphi }_2=\frac{\rho (2)-\rho (1{)}^2}{1-\rho (1{)}^2}.\end{array}$$