Prais-Winsten estimator with AR(1) errors

When to use the estimator

If the errors of a linear regression follow an AR(1) process, the coefficients of ordinary least squares are still unbiased, but they are no longer efficient and their standard errors are wrong. The Prais-Winsten estimator (Prais and Winsten, 1954) corrects for this by transforming the data so that the transformed errors are serially uncorrelated, and by estimating the coefficients and the autocorrelation of the errors in turn until the latter converges.

The estimator assumes that the regressors are strictly exogenous. A model that contains a lagged dependent variable does not satisfy this, and its estimates would be inconsistent. The transformation also assumes an AR(1) process, not a higher order one.

The Cochrane-Orcutt estimator solves the same problem by dropping the first observation of the sample. Prais-Winsten keeps it and rescales it instead, which matters in short samples.

library(prais)
#> Loading required package: sandwich
#> Loading required package: pcse
#> 
#> Attaching package: 'pcse'
#> The following object is masked from 'package:sandwich':
#> 
#>     vcovPC

A simple time series model

The package ships the barium data, monthly observations from the U.S. chemical industry. The column t is the time index.

data("barium")

str(barium[, c("lchnimp", "lchempi", "lgas", "lrtwex", "t")])
#> 'data.frame':    131 obs. of  5 variables:
#>  $ lchnimp: num  5.4 4.55 5.39 5.76 4.74 ...
#>  $ lchempi: num  4.61 4.61 4.62 4.63 4.65 ...
#>  $ lgas   : num  22.8 22.9 22.9 22.9 22.9 ...
#>  $ lrtwex : num  4.46 4.45 4.45 4.47 4.46 ...
#>  $ t      : num  1 2 3 4 5 6 7 8 9 10 ...

The model is specified as it would be for lm. Argument index names the time variable, by which the observations are ordered before the transformation is applied. Without it the estimator would have to assume that the rows already are in the right order.

pw <- prais_winsten(lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6,
                    data = barium, index = "t")
#> Iteration 0: rho = 0
#> Iteration 1: rho = 0.2708
#> Iteration 2: rho = 0.291
#> Iteration 3: rho = 0.293
#> Iteration 4: rho = 0.2932
#> Iteration 5: rho = 0.2932
#> Iteration 6: rho = 0.2932
#> Iteration 7: rho = 0.2932

The iterations are reported as they happen, which shows how the estimate of the AR(1) coefficient rho settles down. They are ordinary messages and can be silenced with suppressMessages.

Reading the output

summary(pw)
#> 
#> Call:
#> prais_winsten(formula = lchnimp ~ lchempi + lgas + lrtwex + befile6 + 
#>     affile6 + afdec6, data = barium, index = "t")
#> 
#> Residuals of the transformed model:
#>      Min       1Q   Median       3Q      Max 
#> -2.01146 -0.39152  0.06758  0.35063  1.35021 
#> 
#> AR(1) coefficient rho after 7 iterations: 0.2932
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -37.07582   22.77843  -1.628   0.1061    
#> lchempi       2.94096    0.63284   4.647 8.46e-06 ***
#> lgas          1.04630    0.97734   1.071   0.2864    
#> lrtwex        1.13277    0.50666   2.236   0.0272 *  
#> befile6      -0.01648    0.31938  -0.052   0.9589    
#> affile6      -0.03316    0.32181  -0.103   0.9181    
#> afdec6       -0.57681    0.34199  -1.687   0.0942 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.5733 on 124 degrees of freedom
#> Multiple R-squared:  0.2021, Adjusted R-squared:  0.1635 
#> F-statistic: 5.235 on 6 and 124 DF,  p-value: 7.764e-05
#> 
#> Durbin-Watson statistic (original): 1.458 
#> Durbin-Watson statistic (transformed): 2.087

Three parts of the output are specific to this estimator.

  • The AR(1) coefficient and the number of iterations it took to converge. A value clearly different from zero is what justifies using the estimator in the first place.
  • The residuals of the transformed model, rather than those of the original one. The transformed model is what is actually fitted by least squares, so its residuals are the ones whose properties matter.
  • Two Durbin-Watson statistics. The first belongs to the original model and the second to the transformed one. A value near 2 indicates no first order serial correlation, so the pair shows whether the transformation did its job.

It is worth comparing the result with ordinary least squares.

ols <- lm(lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6,
          data = barium)

round(cbind("OLS" = coef(ols),
            "PW" = coef(pw),
            "SE (OLS)" = sqrt(diag(vcov(ols))),
            "SE (PW)" = sqrt(diag(vcov(pw)))), 4)
#>                  OLS       PW SE (OLS) SE (PW)
#> (Intercept) -17.8020 -37.0758  21.0455 22.7784
#> lchempi       3.1172   2.9410   0.4792  0.6328
#> lgas          0.1963   1.0463   0.9066  0.9773
#> lrtwex        0.9830   1.1328   0.4002  0.5067
#> befile6       0.0596  -0.0165   0.2610  0.3194
#> affile6      -0.0324  -0.0332   0.2643  0.3218
#> afdec6       -0.5652  -0.5768   0.2858  0.3420

Some coefficients move a good deal, lgas most of all. Every standard error is larger than its counterpart from least squares, here by between eight and thirty per cent, because least squares treats positively correlated errors as though they carried more independent information than they do. Which of the two effects dominates depends on the data; the point is that the corrected standard errors are the ones to report.

Confidence intervals and predictions

confint uses the t distribution with the residual degrees of freedom of the model, so the intervals agree with the p-values of the summary.

confint(pw)
#>                   2.5 %    97.5 %
#> (Intercept) -82.1607236 8.0090737
#> lchempi       1.6883996 4.1935273
#> lgas         -0.8881332 2.9807305
#> lrtwex        0.1299595 2.1355894
#> befile6      -0.6486200 0.6156637
#> affile6      -0.6701089 0.6037933
#> afdec6       -1.2536974 0.1000749

predict returns the conditional mean of the model. With interval = "confidence" the interval of that mean is added. The AR(1) structure of the errors is not used for the prediction, which is why an interval for an individual future observation is not available.

head(predict(pw, interval = "confidence"))
#>           fit      lwr      upr
#> [1,] 5.362201 5.017814 5.706587
#> [2,] 5.495592 5.186166 5.805019
#> [3,] 5.453788 5.151540 5.756037
#> [4,] 5.612285 5.290419 5.934151
#> [5,] 5.638605 5.338288 5.938921
#> [6,] 5.673724 5.334911 6.012538

Passing newdata predicts for observations that were not part of the estimation. The variables of the formula do not have to be transformed beforehand.

Robust standard errors

If the errors are heteroskedastic on top of being serially correlated, the covariance matrix of the transformed model can be made robust with vcovHC.

sqrt(diag(vcovHC(pw, type = "HC1")))
#> (Intercept)     lchempi        lgas      lrtwex     befile6     affile6 
#>  20.8974254   0.5995493   0.9251510   0.4951266   0.3277787   0.2772969 
#>      afdec6 
#>   0.4225520

A panel model

For panel data, index takes two names, the ID variable first and the time variable second. The transformation is then applied within each panel, so that the last observation of one panel is never differenced against the first observation of the next.

A small balanced panel with a known AR(1) coefficient makes it easy to see what the estimator recovers.

set.seed(123)

n_id <- 8
n_time <- 25
rho <- .6

panel <- do.call(rbind, lapply(seq_len(n_id), function(i) {
  u <- rnorm(n_time, 0, 2)
  for (j in 2:n_time) {
    u[j] <- u[j] + rho * u[j - 1]
  }
  x <- rnorm(n_time, 10, 3)
  data.frame(id = i, year = seq_len(n_time), x = x, y = 5 + 2 * x + u)
}))

head(panel)
#>   id year         x        y
#> 1  1    1  4.939920 13.75889
#> 2  1    2 12.513361 28.89380
#> 3  1    3 10.460119 28.35790
#> 4  1    4  6.585589 19.77479
#> 5  1    5 13.761445 33.74363
#> 6  1    6 11.279393 31.72136
pw_panel <- prais_winsten(y ~ x, data = panel, index = c("id", "year"))
#> Iteration 0: rho = 0
#> Iteration 1: rho = 0.5253
#> Iteration 2: rho = 0.5347
#> Iteration 3: rho = 0.5348
#> Iteration 4: rho = 0.5348

summary(pw_panel)
#> 
#> Call:
#> prais_winsten(formula = y ~ x, data = panel, index = c("id", 
#>     "year"))
#> 
#> Residuals of the transformed model:
#>     Min      1Q  Median      3Q     Max 
#> -5.0152 -1.2660 -0.1195  1.2123  6.2139 
#> 
#> AR(1) coefficient rho after 4 iterations: 0.5348
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  5.20326    0.51232   10.16   <2e-16 ***
#> x            1.99710    0.04285   46.61   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.967 on 198 degrees of freedom
#> Multiple R-squared:  0.9207, Adjusted R-squared:  0.9203 
#> F-statistic:  2298 on 1 and 198 DF,  p-value: < 2.2e-16
#> 
#> Durbin-Watson statistic (original): 0.9212 
#> Durbin-Watson statistic (transformed): 1.877

The coefficient on x is recovered almost exactly, and the Durbin-Watson statistic of the transformed model is much closer to 2 than that of the original one. The estimate of rho lands below the 0.6 the errors were generated with. That is the usual downward bias of an autocorrelation estimated from the residuals of a short series, which is of the order of one divided by the number of periods per panel, and it shrinks as the panels get longer. The correction does not depend on hitting rho exactly.

One coefficient per panel

By default a single rho is estimated from all panels together. With panelwise = TRUE every panel gets its own, which is appropriate if the panels differ in how persistent their errors are.

pw_pw <- prais_winsten(y ~ x, data = panel, index = c("id", "year"),
                       panelwise = TRUE)

round(pw_pw$rho[nrow(pw_pw$rho), ], 3)
#>     1     2     3     4     5     6     7     8 
#> 0.601 0.415 0.785 0.041 0.618 0.610 0.652 0.191

Panel specific estimates come from short series and are therefore noisy. Argument rhoweight averages them into a single coefficient, weighted by the length of each panel, which trades some flexibility for a more stable estimate.

pw_t1 <- prais_winsten(y ~ x, data = panel, index = c("id", "year"),
                       panelwise = TRUE, rhoweight = "T1")

pw_t1$rho[nrow(pw_t1$rho), ]
#> [1] 0.4988328

Panel-corrected standard errors

Errors are often correlated across panels within the same period. vcovPC implements the panel-corrected standard errors of Beck and Katz (1995) for that case.

sqrt(diag(vcovPC(pw_panel)))
#> (Intercept)           x 
#>  0.48829896  0.04036159

pairwise = FALSE, the default, uses only the periods that all panels have in common. pairwise = TRUE uses every period that a given pair of panels shares, which matters for unbalanced panels.

Gaps in the time variable

Two observations that lie k periods apart are correlated by rho^k under an AR(1) process. The transformation accounts for this, so a panel with missing periods is not treated as if its observations were consecutive.

gapped <- panel[panel$year != 10, ]

pw_gap <- prais_winsten(y ~ x, data = gapped, index = c("id", "year"))
#> Iteration 0: rho = 0
#> Iteration 1: rho = 0.5182
#> Iteration 2: rho = 0.5312
#> Iteration 3: rho = 0.5313
#> Iteration 4: rho = 0.5313
#> Iteration 5: rho = 0.5313

pw_gap$rho[nrow(pw_gap$rho), ]
#> [1] 0.5312979

The periods are counted in steps of the greatest common divisor of the differences of the time variable, so quarterly, monthly or biennial data need no further arguments, and Date variables are supported.

Models estimated with plm

A pooled panel model that was fitted with plm::plm can be handed to prais_winsten in place of the formula, the data and the index, which are then taken from the model.

data("Grunfeld", package = "plm")

pooled <- plm::plm(inv ~ value + capital, data = Grunfeld,
                   index = c("firm", "year"), model = "pooling")

coef(suppressMessages(prais_winsten(pooled)))
#>  (Intercept)        value      capital 
#> -29.80659971   0.09120477   0.29383533

Only models estimated with model = "pooling" are supported, because the Prais-Winsten estimator obtains all of its estimates by ordinary least squares. The within and the random effects estimator apply a transformation of their own, which does not leave the AR(1) structure of the errors intact. Fixed effects can be obtained by putting the dummy variables into the formula, as in y ~ x + factor(id).

Tidy output

Methods for the generics of broom are available, which packages for publishable regression tables use to read a model.

broom::tidy(pw_panel, conf.int = TRUE)
#> # A tibble: 2 × 7
#>   term        estimate std.error statistic   p.value conf.low conf.high
#>   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
#> 1 (Intercept)     5.20    0.512       10.2 8.95e- 20     4.19      6.21
#> 2 x               2.00    0.0428      46.6 1.08e-108     1.91      2.08

References

Beck, N. L. and Katz, J. N. (1995): What to do (and not to do) with time-series cross-section data. American Political Science Review 89, 634-647.

Prais, S. J. and Winsten, C. B. (1954): Trend Estimators and Serial Correlation. Cowles Commission Discussion Paper, 383 (Chicago).

Wooldridge, J. M. (2016): Introductory Econometrics. A Modern Approach. 6th ed. Mason, OH: South-Western Cengage Learning.