Breadcrumb

8.1 Linear Regression Models with Autoregressive Errors

When we do regressions using time series variables, it is common for the errors (residuals) to have a time series structure. This violates the usual assumption of independent errors made in ordinary least squares regression. The consequence is that the estimates of coefficients and their standard errors will be wrong if the time series structure of the errors is ignored.

It is possible, though, to adjust estimated regression coefficients and standard errors when the errors have an AR structure. More generally, we will be able to make adjustments when the errors have a general ARIMA structure.

The Regression Model with AR Errors

Suppose that yt and xt are time series variables. A simple linear regression model with autoregressive errors can be written as

If we let \(\Phi(B)=1-\phi_{1}B- \phi_{2}B^2 - \cdots\), then we can write the AR model for the errors as

\(\Phi(B)\epsilon_{t}=w_{t}.\)

If we assume that an inverse operator, \(\Phi^{-1}(B)\), exists, then \(\epsilon_{t}=\Phi^{-1}(B)w_{t}\).

So, the model can be written

\(y_t =\beta_{0} +\beta_{1}x_{t} + \Phi^{-1}(B)w_{t},\)

where \(w_{t}\) is the usual white noise series.

Note!
This generalizes to the multiple linear regression structure as well. We can have more than one x-variable (time series) on the right side of the equation. Each x-variable is adjusted in the manner described below.

Examining Whether This Model May be Necessary

Start by doing an ordinary regression. Store the residuals.

Analyze the time series structure of the residuals to determine if they have an AR structure.

If the residuals from the ordinary regression appear to have an AR structure, estimate this model and diagnose whether the model is appropriate.

Note!
In the last of the equations just given \(\beta_{0}\) is an unknown constant that does not move in time. Thus it is not affected by any backshift operation. That’s why the B operations were not applied in that equation.

Using the items just defined, we can write the model as

\((2) \;\;\; y^*_t =\beta^*_0 +\beta_1x^*_t + w_{t}\)

Remember that \(w_t\) is a white noise series, so this is just the usual simple regression model.

We use results from model (2) to iteratively adjust our estimates of coefficients in model (1).

The estimated slope \(\widehat{\beta}_1\) from model (2) will be the adjusted estimate of the slope in model (1) (and its standard error from this model will be correct as well).

To adjust our estimate of the intercept \(\beta_0\) in model (1), we use the relationship \(\beta^*_0=(1 - \phi_{1} - \cdots - \phi_{p})\beta_{0}\) to determine that

This procedure is attributed to Cochrane and Orcutt (1949) and is repeated until the estimates converge, that is we observe a very small difference in our estimates between iterations. When the errors exhibit an AR(1) pattern, the cochrane.orcutt function found within the orcutt package in R iterates this procedure.

Additional Comment

For a higher order AR, the adjustment variables are calculated in the same manner with more lags. For instance, suppose the residuals were found to have an AR(2) with estimated coefficients 0.9 and -0.2. Then the y- and x- variables for the adjustment regression would be

There are n = 76 annual observations in sequence. Suppose x is a leading economic indicator (predictor) for a country and y = a measure of the state of the economy. The following plot shows the relationship between x and y for 76 years.

Suppose that we want to estimate the linear regression relationship between y and x at concurrent times.

Step 1

Estimate the usual regression model. Results from R are:

Coefficients:

Estimate

Std. Error

t value

Pr(>|t|)

(Intercept)

3.600e+04

7.152e+01

503.373

< 2e-16 ***

x

1.610e+00

2.023e-01

7.956

1.57e-11 ***

(It is not shown, but worthy to note that trend=time(x) is not significant in this regression and x is not changing over time.)

Step 2

Examine the AR structure of the residuals. Following are the ACF and PACF of the residuals. It looks like the errors from Step 1 have an AR(1) structure.

Step 3

Estimate the AR coefficients (and make sure that the AR model actually fits the residuals). For this example, the R estimate of the AR(1) coefficient is:

Coefficients:

ar1

0.6488

s.e.

0.0875

Model diagnostics (not shown here) were okay.

Step 4

Calculate variables to use in the adjustment regression:

\(x^*_t = x_t - 0.6488x_{t-1}\)

\(y^*_t = y_t - 0.6488y_{t-1}\)

Step 5

Use ordinary regression to estimate the model \(y^*_t = \beta^*_0 + \beta_1 x^*_t +w_t\). For this example, the results are

Coefficients:

Estimate

Std. Error

t value

Pr(>|t|)

(Intercept)

1.264e+04

1.476e+01

856.50

< 2e-16 ***

xnew

1.635e+00

1.174e-01

13.93

< 2e-16 ***

The slope estimate (1.64) and its standard error (0.1174) are the adjusted estimates for the original model.

The adjusted estimate of the intercept of the original model is 12640/(1-0.6488) = 35990. The estimated standard error of the intercept is 14.76/(1-0.6488) = 42.03.

This procedure is iterated until the estimates converge. In R, the Cochrane.orcutt function iterates the steps:

Estimating the Coefficients of the Adjusted Regression Model with Maximum Likelihood

The method used here depends upon what program you're using. In R (with gls and arima) and in SAS (with PROC AUTOREG) it's possible to specify a regression model with errors that have an ARIMA structure. With a package that includes regression and basic time series procedures, it's relatively easy to use an iterative procedure to determine adjusted regression coefficient estimates and their standard errors. Remember, the purpose is to adjust "ordinary" regression estimates for the fact that the residuals have an ARIMA structure.

Carrying out the Procedure

Note!We are modeling a potential trend over time with \(\beta_0 + \beta_1t\) and may need to de-trend x if x is strongly correlated with time, \(t\). Please see Section 2.2 of the text for a discussion on de-trending.

Step 2

Examine the ARIMA structure (if any) of the sample residuals from the model in step 1.

Step 3

If the residuals do have an ARIMA structure, use maximum likelihood to simultaneously estimate the regression model using ARIMA estimation for the residuals.

Step 4

Examine the ARIMA structure (if any) of the sample residuals from the model in step 3. If white noise is present, then the model is complete. If not, continue to adjust the ARIMA model for the errors until the residuals are white noise.

Trend is not significant, however we notice that x is increasing over time as well and is strongly correlated with time, r = .998. To reduce the multicollinearity between x and t, we can detrend x and examine the regression for y on the de-trended x, dtx, and time, t. The regression results are now:

Following are the ACF and PACF of the residuals. Because both the ACF and PACF spike and then cut off, we should compare AR(1), MA(1), and ARIMA(1,0,1). We will continue with the MA(1) model in the notes.

Step 3: Estimate the adjusted model with a MA(1) structure for the residuals (and make sure that the MA model actually fits the residuals)

Note that in this example it might work better to use an ARIMA model as we have a univariate time series, but we’ll use the approach of these notes for illustrative purposes. We analyze the glacial varve data described in Example 2.5, page 62 of the text. The response is a measure of the thickness of deposits of sand and silt (varve) left by spring melting of glaciers about 11,800 years ago. The data are annual estimates of varve thickness at a location in Massachusetts for 455 years beginning 11,834 years ago. There are nonconstant variance and outlier problems in the data, so we might try a log transformation with hopes of stabilizing the variance and diminishing the effects of outliers. Here’s a time series plot of the log10 series.

The autocorrelation and partial autocorrelation functions of the residuals from this model follow.

Similar to example 1, we might interpret the patterns either as an ARIMA(1,0,1), an AR(1), or a MA(1). We’ll pick the AR(1) – in large part to show an alternative to the MA(1) in Example 2.

The ARIMA results for a AR(1):

Coefficients:

ar1

intercept

trend

trend2

0.2810

1.2200

9e-04

8e-06

s.e.

0.0466

0.0202

1e-04

4e-05

sigma^2 estimated as 0.04176: log likelihood = 76.86, aic = -143.72

Check diagnostics:

The autocorrelation and partial autocorrelation functions of the residuals from this estimated model include no significant values. The fit seems is satisfactory, so we’ll use these results as our final model. The estimated model is