ARIMA and Seasonal ARIMA Models

Coefficients:

Suppose that the model for your data is white noise. If this is true for every $t$, then it is true for $t-1$ as well, in other words:

$x_t = w_t$

$x_{t-1} = w_{t-1}$

Let's multiply both sides of the second equation by 0.5:

$0.5x_{t-1} = 0.5w_{t-1}$

Next, we will move both terms over to one side:

$0 = -0.5x_{t-1} + 0.5w_{t-1}$

Because the data is white noise, $x_t=w_t$, so we can add $x_t$ to the left side and $w_t$ to the right side:

$x_t = -0.5x_{t-1} + w_t + 0.5w_{t-1}$

This is an ARMA(1, 1)! The problem is that we know it is white noise because of the original equations. If we looked at the ACF, what would we see? You would see the ACF corresponding to white noise, a spike at zero, and then nothing else. This also means if we take the white noise process and you try to fit in an ARMA(1, 1), R will do it and will come up with coefficients that look something like what we have above. This is one of the reasons why we need to look at the ACF and the PACF plots, and other diagnostics. We prefer a model with the fewest parameters. This example also says that for certain parameter values, ARMA models can appear very similar to one another.

R Code for Example 3-1

Here's how we accomplished the work for the example in this lesson.

We first loaded the astsa library discussed in Lesson 1. It's a set of scripts written by Stoffer, one of the textbook's authors. If you installed the astsa package during Lesson 1, then you only need the library command.

Use the command library("astsa"). This makes the downloaded routines accessible.

The session for creating Example 1 then proceeds as follows:

xerie = scan("eriedata.dat") #reads the data
xerie = ts (xerie) # makes sure xerie is a time series object
plot (xerie, type = "b") # plots xerie
acf2 (xerie) # author's routine for graphing both the ACF and the PACF
sarima (xerie, 1, 0, 0) # this is the AR(1) model estimated with the author's routine
sarima (xerie, 0, 0, 1) # this is the incorrect MA(1) model
sarima (xerie, 1, 0, 1) # this is the over-parameterized ARMA(1,1) model 

sarima.for(xerie, 4, 1, 0, 0) # four forecasts from an AR(1) model for the erie data 

When you request a graph of the ACF values, "significance" limits are shown by R and by Minitab. In general, the limits for the autocorrelation are placed at $0 ± 2$ standard errors of $r_k$. The formula used for standard error depends upon the situation.

• Within the ACF of residuals as part of the ARIMA routine, the standard errors are determined, assuming the residuals are white noise. The approximate formula for any lag is that s.e. of $r_k=1/\sqrt{n}$.
• For the ACF of raw data (the ACF command), the standard error at a lag k is found as if the right model was an MA(k-1). This allows the possible interpretation that if all autocorrelations past a certain lag are within limits, the model might be an MA of order defined by the last significant autocorrelation.

What are standardized residuals in a time series framework? One of the things that we need to look at when we look at the diagnostics from a regression fit is a graph of the standardized residuals. Let's review what this is for regular regression where the standard deviation is $\sigma$. The standardized residual at observation i

$\dfrac{y_i - \beta_0 - \sum_{j=1}^{p}\beta_j x_{ij}}{\sigma},$

should be N(0, 1). We hope to see normality when we look at the diagnostic plots. Another way to think about this is:

$\dfrac{y_i - \beta_0 -\sum_{j=1}^{p}\beta_j x_{ij}}{\sigma} \approx \dfrac{y_i - \widehat{y}_i}{\sqrt{Var(y_i - \widehat{y}_i)}}.$

Now, with time series, things are very similar:

$\dfrac{x_t-\tilde{x_t}}{\sqrt{P^{t-1}_t}},$

where

$\tilde{x_t} = E(x_t|x_{t-1}, x_{t-2}, \dots) \text{ and } P^{t-1}_t = E\left[(x_t-\tilde{x_t})^2\right].$

This is where the standardized residuals come from. This is also essentially how a time series is fit using R. We want to minimize the sums of these squared values:

$\sum_{t=1}^{n}\left(\dfrac{x_t - \tilde{x_t}}{\sqrt{P^{t-1}_t}}\right)^2$

(In reality, it is slightly more complicated. The log-likelihood function is minimized, and this is one term of that function).

Consider the AR(2) model $x_t=\delta+\phi_1x_{t-1}+\phi_2x_{t-2}+w_t$. In this model, $x_{t}$ is a linear function of the values of x at the previous two times. Suppose that we have observed n data values and wish to use the observed data and estimated AR(2) model to forecast the value of $x_{n+1}$ and $x_{n+2}$, the values of the series at the next two times past the end of the series. The equations for these two values are

$x_{n+1}=\delta+\phi_1x_n+\phi_2x_{n-1}+w_{n+1}$

$x_{n+2}=\delta+\phi_1x_{n+1}+\phi_2x_n+w_{n+2}$

To use the first of these equations, we simply use the observed values of $x_{n }$and $x_{n-1 }$and replace $w_{n+1}$ by its expected value of 0 (the assumed mean for the errors).

The second equation for forecasting the value at time n + 2 presents a problem. It requires the unobserved value of $x_{n+1}$ (one time past the end of the series). The solution is to use the forecasted value of $x_{n+1}$ (the result of the first equation).

In general, the forecasting procedure, assuming a sample size of n, is as follows:

• For any $w_{j}$ with 1 ≤ j ≤ n, use the sample residual for time point j
• For any $w_{j}$ with j > n, use 0 as the value of $w_{j}$

• For any $x_{j}$ with 1 ≤ j ≤ n, use the observed value of $x_{j}$
• For any $x_{j }$with j > n, use the forecasted value of $x_{j}$
  
  

Our authors use the notation $x^n_{n+m}$ to represent a forecast m times past the end of the observed series. The "superscript" is to be read as "given data up to time n". Other authors use the notation $x_{n} \left(m \right)$ to denote a forecast m times past time n.

To understand the formula for the standard error of the forecast error, we first need to define the concept of psi-weights.

Any ARIMA model can be converted to an infinite order MA model:

\begin{align} x_t - \mu & =  w_t + \Psi_1w_{t-1} + \Psi_2w_{t-2} + \dots + \Psi_kw_{t-k} + \dots \\ & =  \sum_{j=0}^{\infty}\Psi_jw_{t-j} \text{ where} \Psi_0 \equiv 1 \end{align}

An important constraint so that the model doesn't "explode" as we go back in time is

$\sum_{j=0}^{\infty}|\Psi_j| < \infty .$

The process of finding the "psi-weight" representation can involve a few algebraic tricks. Fortunately, R has a routine, ARMAtoMA, that will do it for us. To illustrate how psi-weights may be determined algebraically, we'll consider a simple example.

Suppose that an AR(1) model is $x _ { t } = 40 + 0.6 x _ { t - 1 } + w _ { t }$

For an AR(1) model, the mean $\mu=\delta/(1-\phi_1)$ so in this case, $\mu = 40/(1 - .6) = 100$. We'll define $z _ { t } = x _ { t } - 100$ and rewrite the model as $z _ { t } = 0.6 z _ { t - 1 } + w _ { t }$. (You can do the algebra to check that things match between the two expressions of the model.)

To find the psi-weight expression, we'll continually substitute for the z on the right side in order to make the expression become one that only involves w values.

$z _ { t } = 0.6 z _ { t - 1 } + w _ { t } , \text { so } z _ { t - 1 } = 0.6 z _ { t - 2 } + w _ { t - 1 }$

Substitute the right side of the second expression for $z_{t - 1}$ in the first expression.

This gives

$z _ { t } = 0.6 \left( 0.6 z _ { t - 2 } + w _ { t - 1 } \right) + w _ { t } = 0.36 z _ { t - 2 } + 0.6 w _ { t - 1 } + w _ { t }$.

Next, note that $z _ { t - 2 } = 0.6 z _ { t - 3 } + w _ { t - 2 }$.

Substituting this into the equation gives

$z _ { t } = 0.216 z _ { t - 3 } + 0.36 w _ { t - 2 } + 0.6 w _ { t - 1 } + w _ { t }$.

If you keep going, you'll soon see that the pattern leads to

$z_t = x_t -100 = \sum_{j=0}^{\infty}(0.6)^jw_{t-j}$

Thus the psi-weights for this model are given by $\psi_j=(0.6)^j$ for $j=0,1,\ldots,\infty$.

In R, the command ARMAtoMA(ar = .6, ma=0, 12) gives the first 12 psi-weights. This will give the psi-weights $\psi_1$ to $\psi_{12}$ in scientific notation.

For the AR(1) with AR coefficient = 0.6, they are:

Remember that $\psi_0 \equiv 1$. R doesn't give this value. Its listing starts with $\psi_1$, which equals 0.6 in this case.

MA Models: The psi-weights are easy for an MA model because the model is already written in terms of errors. The psi-weights = 0 for lags past the order of the MA model and equal the coefficient values for lags of the errors that are in the model. Remember that we always have $\psi_0=1$.

Without proof, we'll state a result:

The variance of the difference between the forecasted value at time n + m and the (unobserved) value at time n + m is

Variance of $(x^{n}_{n+m} - x_{n+m}) = \sigma^2_w \sum_{j=0}^{m-1}\Psi^2_j.$

Thus the estimated standard deviation of the forecast error at time n + m is

Standard error of $(x^n_{n+m}-x_{n+m}) = \sqrt{\widehat{\sigma}_w^2\sum_{j=0}^{m-1}\Psi^2_j}.$

When forecasting m = 1 time past the end of the series, the standard error of the forecast error is

Standard error of $(x^n_{n+1}-x_{n+1}) = \sqrt{\widehat{\sigma}_w^2(1)}$

When forecasting the value m = 2 times past the end of the series, the standard error of the forecast error is

Standard error of $(x^n_{n+2}-x_{n+2}) = \sqrt{\widehat{\sigma}_w^2(1+\Psi^2_1)}.$

Notice that the variance will not be too big when m = 1. But, as you predict out farther in the future, the variance will increase. When m is very large, we will get the total variance. In other words, if you are trying to predict very far out, we will get the variance of the entire time series, as if you haven't even looked at what was going on previously.

95% Prediction Interval for $x_{n+m}$

With the assumption of normally distributed errors, a 95% prediction interval for $x_{n+m}$, the future value of the series at time n + m, is

$x^n_{n+m} \pm 1.96 \sqrt{\widehat{\sigma}_w^2\sum_{j=0}^{m-1}\Psi^2_j}.$

Suppose that an AR(1) model is estimated to be $x _ { t } = 40 + 0.6 x _ { t - 1 } + w _ { t }$. This is the same model used earlier in this handout, so the psi-weights we got there apply.

Suppose that we have n = 100 observations, $\widehat{\sigma}_w^2 = 4$ and $x_{100} = 80$. We wish to forecast the values at times 101 and 102 and create prediction intervals for both forecasts.

First, we forecast time 101.

$\begin{array}{lll}x_{101} & = & 40 + 0.6x_{100} + w_{101} \\ x^{100}_{101} & = & 40 +0.6 (80) + 0  =  88  \end{array}$

The standard error of the forecast error at time 101 is

$\sqrt{\widehat{\sigma}_w^2 \sum_{j=0}^{1-1}\psi^2_j} = \sqrt{4(1)} = 2.$

The 95% prediction interval for the value at time 101 is 88 ± 2(1.96), which is 84.08 to 91.96. We are, therefore, 95% confident that the observation at time 101 will be between 84.08 and 91.96. If we repeated this exact process many times, then 95% of the computed prediction intervals would contain the true value of x at time 101.

The forecast for time 102 is

$x^{100}_{102} = 40 + 0.6(88) + 0 = 92.8$

The relevant standard error is

$\sqrt{\widehat{\sigma}_w^2 \sum_{j=0}^{2-1}\psi^2_j} = \sqrt{4(1+0.6^2)} = 2.332$

A 95% prediction interval for the value at time 102 is 92.8 ± (1.96)(2.332).

To forecast using an ARIMA model in R, we recommend our textbook author's script called sarima.for. (It is part of the astsa library recommended previously).

In the homework for Lesson 2, problem 5 asked you to suggest a model for a time series of stride lengths measured every 30 seconds for a runner on a treadmill. 

From R, the estimated coefficients for an AR(2) model and the estimated variance are as follows for a similar data set with n = 90 observations:

The command

sarima.for(stridelength, 6, 2, 0, 0) # 6 forecasts with an AR(2) model for stridelength 

will give forecasts and standard errors of prediction errors for the next six times past the end of the series. Here's the output (slightly edited to fit here):

The forecasts are given in the first batch of values under \$pred, and the standard errors of the forecast errors are given in the last line in the batch of results under \$se.

The procedure also gave this graph, which shows the series followed by the forecasts as a red line and the upper and lower prediction limits as blue dashed lines:

If we wanted to verify the standard error calculations for the six forecasts past the end of the series, we would need to know the psi-weights. To get them, we need to supply the estimated AR coefficients for the AR(2) model to the ARMAtoMA command.

The R command in this case is ARMAtoMA(ar = list(1.148, -0.3359), ma = 0, 5)

This will give the psi-weights $\psi_1$ to $\psi_5$ in scientific notation. The answer provided by R is:

[1] 1.148000e+00 9.820040e-01 7.417274e-01 5.216479e-01 3.497056e-01

(Remember that $\psi_0\equiv 1$ in all cases)

The output for estimating the AR(2) included this estimate of the error variance:

sigma^2 estimated as 11.47

As an example, the standard error of the forecast error for 3 times past the end of the series is

$\sqrt{\widehat{\sigma}^2_w\sum_{j=0}^{3-1}\psi^2_j} = \sqrt{11.47(1+1.148^2+0.982^2)} = 6.1357$

which, except for the round-off error, matches the value of 6.135493 given as the third standard error in the sarima.for output above.

Where will the Forecasts End Up?

For a stationary series and model, the forecasts of future values will eventually converge to the mean and then stay there. Note below what happened with the stride length forecasts when we asked for 30 forecasts past the end of the series. [Command was sarima.for (stridelength, 30, 2, 0, 0)]. The forecast got to 48.74753 and then stayed there.

The graph showing the series and the six prediction intervals is the following

The model includes a non-seasonal MA(1) term, a seasonal MA(1) term, no differencing, no AR terms, and the seasonal period is S = 12.

The non-seasonal MA(1) polynomial is $\theta(B) = 1 + \theta_1 B$. 

The seasonal MA(1) polynomial is $\Theta(B^{12}) = 1 + \Theta_1B^{12}$. 

The model is $(x_t - \mu) = \Theta(B^{12})\theta(B)w_t = (1+\Theta_1B^{12})(1+\theta_1 B)w_t$.  

When we multiply the two polynomials on the right side, we get

$(x_t - \mu) = (1 + \theta_1 B + \Theta_1 B^{12} + \theta_1 \Theta_1 B^{13})w_t$

$= w_t + \theta_1 w_{t-1} + \Theta_1 w_{t-12} + \theta_1 \Theta_1 w_{t-13}$.

Thus the model has MA terms at lags 1, 12, and 13. This leads many to think that the identifying ACF for the model will have non-zero autocorrelations only at lags 1, 12, and 13. There's a slight surprise here. There will also be a non-zero autocorrelation at lag 11. We supply a proof in the Appendix.

We simulated n = 1000 values from an ARIMA $( 0,0,1 ) \times ( 0,0,1 ) _ { 12 }$. The non-seasonal MA (1) coefficient was $\theta_1=0.7$. The seasonal MA(1) coefficient was $\Theta_1=0.6$. The sample ACF for the simulated series was as follows:

The model includes a non-seasonal AR(1) term, a seasonal AR(1) term, no differencing, no MA terms, and the seasonal period is S = 12.

The non-seasonal AR(1) polynomial is $\phi(B) = 1-\phi_1B$. 

The seasonal AR(1) polynomial is $\Phi(B^{12}) = 1-\Phi_1B^{12}$.

The model is $(1-\phi_1B)(1-\Phi_1B^{12})(x_t-\mu) = w_t$.

If we let $z_t = x_t -\mu$ (for simplicity), multiply the two AR components and push all but z_t to the right side we get $z_{t} = \phi_1 z_{t-1} + \Phi_1 z_{t-12} + \left( - \phi_1\Phi_1 \right) z_{t-13} + w_t$.

This is an AR model with predictors at lags 1, 12, and 13.

R can be used to determine and plot the PACF for this model, with $\phi_1$=.6 and $\Phi_1$=.5. That PACF (partial autocorrelation function) is:

It's not quite what you might expect for an AR model, but it almost is. There are distinct spikes at lags 1, 12, and 13, with a bit of action coming before lag 12. Then, it cuts off after lag 13.

R commands were

thepacf=ARMAacf (ar = c(.6,0,0,0,0,0,0,0,0,0,0,.5,-.30),lag.max=30,pacf=T)

1. Step 1: Do a time series plot of the data.

   Examine it for features such as trend and seasonality. You'll know that you've gathered seasonal data (months, quarters, etc.,) so look at the pattern across those time units (months, etc.) to see if there is indeed a seasonal pattern.

2. Step 2: Do any necessary differencing.
3. If there is seasonality and no trend, then take a difference of lag S. For instance, take a 12th difference for monthly data with seasonality. Seasonality will appear in the ACF by tapering slowly at multiples of S. (View the TS plot and ACF/PACF plots for an example of data that requires a seasonal difference. Note that the TS plot shows a clear seasonal pattern that repeats over 12 time points. Seasonal differences are supported in the ACF/PACF of the original data because the first seasonal lag in the ACF is close to 1 and decays slowly over multiples of S=12. Once seasonal differences are taken, the ACF/PACF plots of twelfth differences support a seasonal AR(1) pattern).

1. If there is a linear trend and no obvious seasonality, then take a first difference. If there is a curved trend, consider a transformation of the data before differencing.
2. If there is both trend and seasonality, apply a seasonal difference to the data and then re-evaluate the trend. If a trend remains, then take first differences. For instance, if the series is called x, the commands in R would be:
   
   

   diff12=diff(x, 12)
   plot(diff12)
   acf2(diff12)
   diff1and12 = diff(diff12, 1) 

3. If there is neither an obvious trend nor seasonality, don't take any differences.
4. Step 3: Examine the ACF and PACF of the differenced data (if differencing is necessary).

   We're using this information to determine possible models. This can be tricky going involving some (educated) guessing. Some basic guidance:

   • Non-seasonal terms: Examine the early lags (1, 2, 3, …) to judge non-seasonal terms. Spikes in the ACF (at low lags) with a tapering PACF indicate non-seasonal MA terms. Spikes in the PACF (at low lags) with a tapering ACF indicate possible non-seasonal AR terms.

   • Seasonal terms: Examine the patterns across lags that are multiples of S. For example, for monthly data, look at lags 12, 24, 36, and so on (probably won't need to look at much more than the first two or three seasonal multiples). Judge the ACF and PACF at the seasonal lags in the same way you do for the earlier lags.

5. Step 4: Estimate the model(s) that might be reasonable on the basis of step 3.

   Don't forget to include any differencing that you did before looking at the ACF and PACF. In the software, specify the original series as the data and then indicate the desired differencing when specifying parameters in the arima command that you're using.

6. Step 5: Examine the residuals (with ACF, Box-Pierce, and any other means) to see if the model seems good.

   Compare AIC or BIC values to choose among several models.

   If things don't look good here, it's back to Step 3 (or maybe even Step 2).
   
   

The data series are a monthly series of a measure of the flow of the Colorado River, at a particular site, for n = 600 consecutive months.

Step 1

A time series plot is

With so many data points, it's difficult to judge whether there is seasonality. If it was your job to work on data like this, you probably would know that river flow is seasonal – perhaps likely to be higher in the late spring and early summer, due to snow runoff.

Without this knowledge, we might determine means by month of the year. Below is a plot of means for the 12 months of the year. It's clear that there are monthly differences (seasonality).

Looking back at the time series plot, it's hard to judge whether there's any long-run trend. If there is, it's slight.

Steps 2 and 3

We might try the idea that there is seasonality but no trend. To do this, we can create a variable that gives the 12^th differences (seasonal differences), calculated as $x _ { t } - x _ { t - 12}$. Then, we look at the ACF and the PACF for the 12^th difference series (not the original data). Here they are:

• Non-seasonal behavior: The PACF shows a clear spike at lag 1 and not much else until about lag 11. This is accompanied by a tapering pattern in the early lags of the ACF. A non-seasonal AR(1) may be a useful part of the model.
• Seasonal behavior: We look at what's going on around lags 12, 24, and so on. In the ACF, there's a cluster of (negative) spikes around lag 12 and then not much else. The PACF tapers in multiples of S; that is, the PACF has significant lags at 12, 24, 36, and so on. This is similar to what we saw for a seasonal MA(1) component in Example 1 of this lesson.

Remembering that we're looking at 12^th differences, the model we might try for the original series is ARIMA $( 1,0,0 ) \times ( 0,1,1 ) _ { 12 }$.

Step 4

R results for the ARIMA $( 1,0,0 ) \times ( 0,1,1 ) _ { 12 }$:

Step 5 (diagnostics)

The normality and Box-Pierce test results are shown in Lesson 4.2.  Besides normality, things look good.  The Box-Pierce statistics are all non-significant, and the estimated ARIMA coefficients are statistically significant.

The ACF of the residuals looks good too:

What doesn't look perfect is a plot of residuals versus fits. There's non-constant variance.

We've got three choices for what to do about the non-constant variance: (1) ignore it, (2) go back to step 1 and try a variance stabilizing transformation like log or square root, or (3) use an ARCH model that includes a component for changing variances. We'll get to ARCH models later in the course.

Lesson 4.2 for this week will give R guidance and an additional example or two.

Only those interested in theory need to look at the following.

In Example 4-1, we promised a proof that $\rho_{11}$ ≠ 0 for ARIMA $( 0,0,1 ) \times ( 0,0,1 ) _ { 12 }$.

A correlation is defined as Covariance/ product of standard deviations.

The covariance between $x_t$ and $x _ { t - 11 } = E \left( x _ { t } - \mu \right) \left( x _ { t - 11 } - \mu \right)$.

For the model in Example 1,

$x _ { t } - \mu = w _ { t } + \theta _ { 1 } w _ { t - 1 } + \Theta _ { 1 } w _ { t - 12 } + \theta _ { 1 } \Theta _ { 1 } w _ { t - 13 }$

$x _ { t - 11 } - \mu = w _ { t - 11 } + \theta _ { 1 } w _ { t - 12 } + \Theta _ { 1 } w _ { t - 23 } + \theta _ { 1 } \Theta _ { 1 } w _ { t - 24 }$

The covariance between$x_t$ and $x_{t-11}$

(2) $\mathrm { E } \left( w _ { t } + \theta _ { 1 } w _ { t - 1 } + \Theta _ { 1 } w _ { t - 12 } + \theta _ { 1 } \Theta _ { 1 } w _ { t - 13 } \right) \left( w _ { t - 11 } + \theta _ { 1 } w _ { t - 12 } + \Theta _ { 1 } w _ { t - 23 } + \theta _ { 1 } \Theta _ { 1 } w _ { t - 24 } \right)$

The w's are independent errors. The expected value of any product involving w's with different subscripts will be 0. A covariance between w's with the same subscripts will be the variance of w.

If you inspect all possible products in expression 2, there will be one product with matching subscripts. They have lag t – 12. Thus this expected value (covariance) will be different from 0.

This shows that the lag 11 autocorrelation will be different from 0. If you look at the more general problem, you can find that only lags 1, 11, 12, and 13 have non-zero autocorrelations for the ARIMA$( 0,0,1 ) \times ( 0,0,1 ) _ { 12 }$.

A seasonal ARIMA model incorporates both non-seasonal and seasonal factors in a multiplicative fashion.