Prais–Winsten Algorithm for Regression with Second or Higher Order Autoregressive Errors

Abstract

There is no available Prais–Winsten algorithm for regression with AR(2) or higher order errors, and the one with AR(1) errors is not fully justified or is implemented incorrectly (thus being inefficient). This paper addresses both issues, providing an accurate, computationally fast, and inexpensive generic zig-zag algorithm.

1. Introduction

Both simulation and theoretical evidence show that the estimation of a regression with an autoregressive of order one (AR(1)) errors via exact (conditional) maximum likelihood (ML) is inferior to an exact nonlinear least squares (NLS) estimation (the correctly implemented Prais and Winsten (1954) (PW) method), at least for trending data. Park and Mitchell (1980) show via simulation that, for trending regressors, the PW method is more efficient than the exact/conditional ML (the (Beach and MacKinnon 1978a) (BM) method). This is theoretically confirmed by (Magge 1985; Magee 1989) who approximates biases of various two-step and iterative estimators. Correctly implemented, the PW method delivers exact NLS (equivalent to unconditional ML under normality) estimators. No normality requirement is needed for the PW method. Gurland (1954) formalises and resolves criticism against the estimation method of Cochrane and Orcutt (1949) (CO, hereafter) for losing the first observation,1 which is employed by PW or re-proposed by Kadiyala (1968). Koopmans (1942) was the first to examine the stationary AR(1) correctly. The PW algorithm is a fast zig-zag (computationally inexpensive) algorithm, which retains the first observation and requires no numerical optimiser. However, in order to obtain an (econometrically efficient) exact NLS estimator, the correct (exact) closed form formula must be employed in the iterations for the update of the AR parameter. This is not the case for all available implementations of the PW method (see Magee (1989)). For example, (Judge et al. 1985, chp. 8) propose using the autocorrelation coefficient formula for the update, but this is not correct. Exact/conditional ML fast zig-zag algorithms are provided by (Beach and MacKinnon 1978a, 1978b) for AR(1)2 and AR(2)3 errors, respectively. None is available for higher than two AR orders. The PW algorithm for a regression with AR(1) errors is relatively well known, although usually incorrectly implemented. This is because an incorrect closed form estimator is used for the update. Park and Mitchell (1980) implement the PW algorithm correctly, although they do not give all the details. Furthermore, there is no zig-zag PW algorithm for a regression with AR(2) or higher order errors, and this literature gap is filled in this paper. Reliable standard errors are also calculated that are not affected by the presence of the lagged dependent variable as a regressor. In fact, our proposed method corrects the inefficient CO method and poor estimation via bad implementation of the Gauss–Newton algorithm.

The paper is organised as follows: Section 2 fully discusses the closed form exact NLS estimation of pure AR(1), AR(2), and AR(p) coefficients, while Section 3 examines the iterative exact joint estimation of regression and autoregressive coefficients. Finally, Section 4 concludes.

2. Closed Form Exact NLS: AR(1,2,…,p)

The main ingredient for the correct PW algorithm is the correct4 closed form update(s) for the autoregressive parameter(s), along with exact generalised least squares (GLS) iterations. Let (observed) $y_t$ be generated via the stationary AR(1)

$$y_t=\theta y_{t-1}+v_t,t=2,\cdots ,n,$$

(1)

with

$$y_1=\frac{v_1}{\sqrt{1-{\theta }^2}},$$

(2)

$v_t\sim iid(0,{\sigma }_v^2)$, $t=1,\cdots ,n$. This assumption can be relaxed, and $v_t$ can be a martingale difference sequence with finite fourth moment (see (Stock 1991), for example). Assuming normal $v_t$, the exact/conditional ML estimator of (1), conditional on (2), in closed form, is the solution to their cubic equation that BM provides. PW/exact NLS (eventually in closed form) requires the exact sum of squares to be minimised, that is,

$$\underset{\theta }{min}\{y_1^2(1-{\theta }^2)+\sum _2^n{(y_t-\theta y_{t-1})}^2\},$$

(3)

see (Kmenta 1986, p. 319 (full proof in the present paper)). Here, we make use of ${\sigma }_y/{\sigma }_v=\sqrt{1-{\theta }^2}$. Phillips (1977,1978) and Sawa (1978) (among others) employ an inexact version of (3), ignoring $y_1^2(1-{\theta }^2)$. This results in the so-called LS (also incorrectly called the ML) estimator (${\widehat{\theta }}_{LS}$), which is more biased and less efficient than the exact NLS estimator of (3) (also called the (original) PW estimator by Park and Mitchell (1980); PW2 in (Magee 1989, p. 662)). Minimisation is via

$$-2\theta y_1^2-2\sum _2^n(y_t-\theta y_{t-1})y_{t-1}=0,$$

(4)

which results in the closed form of the (genuine) PW estimator5

$${\widehat{\theta }}_{PW}=\frac{{\sum }_2^n{y}_t{y}_{t-1}}{{\sum }_3^n{y}_{t-1}^2},$$

(5)

which is better than the OLS estimator ${\widehat{\theta }}_{LS}={\sum }_2^n{y}_t{y}_{t-1}/{\sum }_2^n{y}_{t-1}^2$.6 An even worse estimator is the autocorrelation coefficient ${\sum }_2^n{y}_t{y}_{t-1}/{\sum }_1^n{y}_t^2$.7

For the stationary AR(2) model

$$y_t={\theta }_1{y}_{t-1}+{\theta }_2{y}_{t-2}+v_t,t=3,\cdots ,n,$$

(6)

$$y_1=v_1/(\sqrt{1-{\theta }_2^2}\sqrt{1-{\rho }_1^2}),y_2={\rho }_1{y}_1+v_2/\sqrt{1-{\theta }_2^2},{\rho }_1={\theta }_1/(1-{\theta }_2),$$

(7)

also ${\rho }_2={\theta }_1{\rho }_1+{\theta }_2$,8 no closed-form exact NLS/PW estimators of the autoregressive parameters are available in the literature. To derive them, the exact sum of squares to be minimised is

$$\underset{{\theta }_1,{\theta }_2}{min}\{(1-{\theta }_2^2)(1-{\rho }_1^2)y_1^2+(1-{\theta }_2^2){(y_2-{\rho }_1{y}_1)}^2+\sum _3^n{(y_t-{\theta }_1{y}_{t-1}-{\theta }_2{y}_{t-2})}^2\}.$$

(8)

Use has been made of the fact that in this case, ${\sigma }_y/{\sigma }_v=\sqrt{1-{\theta }_2^2}\sqrt{1-{\rho }_1^2}$. The required canonical equations (with the help of MATHEMATICA™ and imposing ${\rho }_1={\theta }_1/(1-{\theta }_2)$) are

$$\begin{array}{ccc}\hfill (1+{\theta }_2)y_1{y}_2+\sum _3^n{y}_t{y}_{t-1}-{\theta }_1\sum _3^n{y}_{t-1}^2-{\theta }_2\sum _3^n{y}_{t-1}{y}_{t-2}& =\hfill & \hfill 0\\ \hfill {\theta }_1{y}_1{y}_2+{\theta }_2(y_1^2+y_2^2)+\sum _3^n{y}_t{y}_{t-2}-{\theta }_1\sum _3^n{y}_{t-1}{y}_{t-2}-{\theta }_2\sum _3^n{y}_{t-2}^2& =\hfill & \hfill 0\end{array}$$

(9)

that become (after manipulations)

$$\begin{array}{ccc}\hfill (y_1{y}_2+\sum _3^n{y}_t{y}_{t-1})-{\theta }_1\sum _3^n{y}_{t-1}^2-{\theta }_2(-y_1{y}_2+\sum _3^n{y}_{t-1}{y}_{t-2})& =\hfill & \hfill 0\\ \hfill \sum _3^n{y}_t{y}_{t-2}-{\theta }_1(-y_1{y}_2+\sum _3^n{y}_{t-1}{y}_{t-2})-{\theta }_2(-y_1^2-y_2^2+\sum _3^n{y}_{t-2}^2)& =\hfill & \hfill 0.\end{array}$$

(10)

Re-arranging (10), the system of equations can be solved for the brand new (efficient/genuine) PW estimators of the AR(2) parameters. That is,

$$\left[\begin{array}{c}{\widehat{\theta }}_{1PW}\\ {\widehat{\theta }}_{2PW}\end{array}\right]={\left[\begin{array}{cc}\sum _3^n{y}_{t-1}^2& \sum _4^n{y}_{t-1}{y}_{t-2}\\ \sum _4^n{y}_{t-1}{y}_{t-2}& \sum _5^n{y}_{t-2}^2\end{array}\right]}^{-1}\left[\begin{array}{c}\sum _2^n{y}_t{y}_{t-1}\\ \sum _3^n{y}_t{y}_{t-2}\end{array}\right].$$

(11)

In view of (11), we can guess the brand new closed-form PW autoregressive estimators vector of any AR(p)

$$y_t={\theta }_1{y}_{t-1}+{\theta }_2{y}_{t-2}+\cdots +{\theta }_p{y}_{t-p}+v_t,t=p+1,\cdots ,n,$$

along with the required equations for observations 1,…,p (not stated as they are difficult). That is,

$$\left[\begin{array}{c}{\widehat{\theta }}_{1PW}\\ {\widehat{\theta }}_{2PW}\\ ⋮\\ {\widehat{\theta }}_{p-1PW}\\ {\widehat{\theta }}_{pPW}\end{array}\right]=$$

$${\left[\begin{array}{ccccc}\sum _3^n{y}_{t-1}^2& \sum _4^n{y}_{t-1}{y}_{t-2}& \dots & \sum _{1+p}^n{y}_{t-1}{y}_{t-p+1}& \sum _{2+p}^n{y}_{t-1}{y}_{t-p}\\ \sum _4^n{y}_{t-1}{y}_{t-2}& \sum _5^n{y}_{t-2}^2& \dots & \sum _{2+p}^n{y}_{t-2}{y}_{t-p+1}& \sum _{3+p}^n{y}_{t-2}{y}_{t-p}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \sum _{1+p}^n{y}_{t-1}{y}_{t-p+1}& \sum _{2+p}^n{y}_{t-2}{y}_{t-p+1}& \dots & \sum _{2p-1}^n{y}_{t-p+1}^2& \sum _{2p}^n{y}_{t-p+1}{y}_{t-p}\\ \sum _{2+p}^n{y}_{t-1}{y}_{t-p}& \sum _{3+p}^n{y}_{t-2}{y}_{t-p}& \dots & \sum _{2p}^n{y}_{t-p+1}{y}_{t-p}& \sum _{2p+1}^n{y}_{t-p}^2\end{array}\right]}^{-1}$$

$$\times \left[\begin{array}{c}\sum _2^n{y}_t{y}_{t-1}\\ \sum _3^n{y}_t{y}_{t-2}\\ ⋮\\ \sum _p^n{y}_t{y}_{t-p+1}\\ \sum _{p+1}^n{y}_t{y}_{t-p}\end{array}\right].$$

This expression is easily programmable in a matrix programming language.

3. Iteration

The closed-form estimators of the previous section are to be used in iterations for updating. Let observed $Y_t$ be generated via

$$Y_t=x_t^{\prime }\beta +y_t,$$

(12)

with (unobserved) $y_t$ following (1) and (2); $x_t$ is a $k\times 1$ vector of regressors (the t-th row of regressor matrix X), and $\beta$ is the corresponding coefficient vector. The PW iterative algorithm has the following steps: (I) Apply OLS to (12) and derive residual ${\widehat{y}}_t$ and ${\widehat{\theta }}_{PW}$ from (5), replacing $y_t$ with ${\widehat{y}}_t$. (II) Re-estimate $\beta$ from

$$Y_1\sqrt{1-{\widehat{\theta }}_{PW}^2}=x_1^{\prime }\beta \sqrt{1-{\widehat{\theta }}_{PW}^2}+v_1$$

(13)

$$Y_t-{\widehat{\theta }}_{PW}^2{Y}_{t-1}={(x_t-{\widehat{\theta }}_{PW}^2{x}_{t-1})}^{\prime }\beta +v_t,t=2,\cdots ,n.$$

(14)

For a PW algorithm that delivers exact NLS estimators, the first observation must be employed (see (13)), and the $\theta$ update must be via (5). According to (Magee 1989, p. 663), (5) guarantees convergence. (III) Use the new estimate of $\beta$ and obtain new residual ${\widehat{y}}_t$ from (12) and new estimate of $\theta$, and proceed to step (II) if convergence criterion is not met. In addition, estimates of $\theta$ must be forced to be in $(-1,1)$. $P_1\left(\theta \right)$ is the exact GLS matrix of first order with the main diagonal $\sqrt{1-{\theta }^2},1,\cdots ,1$, first sub-diagonal $-\theta ,\cdots ,-\theta$, and zero elsewhere. Let $\tilde{v}$ be the converged innovation residual vector from (13) and (14), corresponding to converged ${\tilde{\theta }}_{PW}$; we define the innovation variance estimate ${\tilde{\sigma }}_v^2={\tilde{v}}^{\prime }\tilde{v}/(n-k)$, and $co{v}_{\tilde{\beta }}={\tilde{\sigma }}_v^2{(X^{\prime }{P}_1{({\tilde{\theta }}_{PW})}^{\prime }{P}_1({\tilde{\theta }}_{PW})X)}^{-1}$ for the converged estimate of $\beta$, $\tilde{\beta }$. For ${\tilde{\theta }}_{PW}$, we rely on the fact that exact NLS estimation is identical to unconditional ML under normality and calculate the quasi-ML covariance $co{v}_{{\tilde{\theta }}_{PW}}$ as the inverse of minus the Hessian of the concentrated unconditional loglikelihood $L\left(\theta \right)=-(n/2)ln{(Y^{*}-X^{*}{\widehat{\beta }}_{GLS,\theta })}^{\prime }(Y^{*}-X^{*}{\widehat{\beta }}_{GLS,\theta })$ evaluated at ${\tilde{\theta }}_{PW}$, with ${\widehat{\beta }}_{GLS,\theta }={(X^{*}{}^{\prime }{X}^{*})}^{-1}(X^{*}{}^{\prime }{Y}^{*})$, $Y^{*}=P_1\left(\theta \right)Y$ and $X^{*}=P_1\left(\theta \right)X$, and the MLE innovation variance estimate ${\widehat{\sigma }}_{vMLE,\theta }^2=\widehat{v}{\left(\theta \right)}^{\prime }\widehat{v}\left(\theta \right)/n$, with “residual” $\widehat{v}\left(\theta \right)$ for a given $\theta$, $\widehat{v}\left(\theta \right)=Y-X{\widehat{\beta }}_{GLS,\theta }$. Alternatively, we may use the asymptotic covariance for ${\tilde{\theta }}_{PW}$, requiring no normality (see below). A third option could be to calculate the sandwich covariance. Note that all these covariances are not affected by the presence of the lagged dependent variable as a regressor.

Similarly, when $y_t$ in (12) follows the AR(2) model, the genuine PW iterative algorithm has the following steps: (I) Apply OLS to (12) and derive residual ${\widehat{y}}_t$ and ${\widehat{\theta }}_{1PW}$ and ${\widehat{\theta }}_{2PW}$ from (11), replacing $y_t$ with ${\widehat{y}}_t$. (II) Re-estimate $\beta$ from

$$Y_1\sqrt{1-{\widehat{\rho }}_1^2}\sqrt{1-{\widehat{\theta }}_{2PW}^2}=x_1^{\prime }\beta \sqrt{1-{\widehat{\rho }}_1^2}\sqrt{1-{\widehat{\theta }}_{2PW}^2}+v_1$$

(15)

$$(Y_2-{\widehat{\rho }}_1{Y}_1)\sqrt{1-{\widehat{\theta }}_{2PW}^2}={(x_2-{\widehat{\rho }}_1{x}_1)}^{\prime }\beta \sqrt{1-{\widehat{\theta }}_{2PW}^2}+v_2$$

(16)

$$\begin{array}{ccc}\hfill Y_t-{\widehat{\theta }}_{1PW}{Y}_{t-1}-{\widehat{\theta }}_{2PW}{Y}_{t-2}& =\hfill & \\ \hfill {(x_t-{\widehat{\theta }}_{1PW}{x}_{t-1}-{\widehat{\theta }}_{2PW}{x}_{t-2})}^{\prime }\beta +v_t,t& =\hfill & \hfill 3,\cdots ,n.\end{array}$$

(17)

Note that ${\widehat{\rho }}_1={\widehat{\theta }}_{1PW}/(1-{\widehat{\theta }}_{2PW})$. (In addition, estimates of ${\theta }_1$ and ${\theta }_2$ must be forced to be inside the stationarity triangle.) (III) Use the new estimator of $\beta$ and obtain new residual ${\widehat{y}}_t$ from (12) and new estimators of ${\theta }_1$ and ${\theta }_2$ from (11), and proceed to step (II) if convergence criterion is not met. $P_2({\theta }_1,{\theta }_2)$ is the exact GLS matrix of second order with the main diagonal $\sqrt{1-{\theta }_2^2}\sqrt{1-{\rho }_1^2},\sqrt{1-{\theta }_2^2},1,\cdots ,1$, first sub-diagonal $-{\rho }_1\sqrt{1-{\theta }_2^2},-{\theta }_1,\cdots ,-{\theta }_1$, second subdiagonal $-{\theta }_2,\cdots ,-{\theta }_2$, and zero elsewhere. Let $\tilde{v}$ be the converged innovation residual from (15)–(17), and define ${\tilde{\sigma }}_v^2={\tilde{v}}^{\prime }\tilde{v}/(n-k)$. For the converged estimate of $\beta$, $\tilde{\beta }$, we define $co{v}_{\tilde{\beta }}={\tilde{\sigma }}_v^2{(X^{\prime }{P}_2{({\tilde{\theta }}_{1PW},{\tilde{\theta }}_{2PW})}^{\prime }{P}_2({\tilde{\theta }}_{1PW},{\tilde{\theta }}_{2PW})X)}^{-1}$, and ${\tilde{\theta }}_{1PW}$ and ${\tilde{\theta }}_{2PW}$ are the converged estimates of ${\theta }_1$ and ${\theta }_2$, respectively. For $\theta =\{{\theta }_1,{\theta }_2\}$, under normality, we calculate the quasi-ML covariance $co{v}_{{\tilde{\theta }}_{PW}}$ from the inverse of minus the Hessian of the concentrated unconditional loglikelihood $L({\theta }_1,{\theta }_2)=-(n/2)ln{(Y^{*}-X^{*}{\widehat{\beta }}_{GLS,\theta })}^{\prime }(Y^{*}-X^{*}{\widehat{\beta }}_{GLS,\theta })$ evaluated at ${\tilde{\theta }}_{1PW}$ and ${\tilde{\theta }}_{2PW}$, with ${\widehat{\beta }}_{GLS,\theta }={(X^{*}{}^{\prime }{X}^{*})}^{-1}(X^{*}{}^{\prime }{Y}^{*})$, $Y^{*}=P_2({\theta }_1,{\theta }_2)Y$ and $X^{*}=P_2({\theta }_1,{\theta }_2)X$, and the MLE innovation variance estimate ${\widehat{\sigma }}_{vMLE,\theta }^2=\widehat{v}{\left(\theta \right)}^{\prime }\widehat{v}\left(\theta \right)/n$, with “residual” $\widehat{v}\left(\theta \right)$ for a given vector $\theta$, $\widehat{v}\left(\theta \right)=Y-X{\widehat{\beta }}_{GLS,\theta }$. Alternatively, we may use the asymptotic covariance for ${\tilde{\theta }}_{PW}=\{{\tilde{\theta }}_{1PW},{\tilde{\theta }}_{2PW}\}$, requiring no normality (see below). A third option could be to calculate the sandwich covariance. Note that all these covariances are not affected by the presence of the lagged dependent variable as a regressor.

Finally, for the model $Y=X\beta +y$ when y follows AR(p) with $\theta ={\{{\theta }_1,\cdots ,{\theta }_p\}}^{\prime }$, OLS provides $\widehat{\beta }$ and $\widehat{y}$, and the first step ${\widehat{\theta }}_{PW}={\{{\widehat{\theta }}_{1PW},\cdots ,{\widehat{\theta }}_{pPW}\}}^{\prime }$, using the generic formula above. Then the OLS (in fact GLS) regression of $P({\widehat{\theta }}_{PW})Y$ on $P({\widehat{\theta }}_{PW})X$ results in a new $\widehat{\beta }$, new $\widehat{y}=Y-X\widehat{\beta }$, and new ${\widehat{\theta }}_{PW}$, repeating until convergence and restricting the elements of $\theta$ accordingly. The exact GLS matrix P has the correct Cholesky decomposition9 of $V_p^{-1}$ in positions $1:p$ and $1:p$, $P[j,j]=1$ for $j=p+1,\cdots ,n$, $P[j+i,j]=-{\theta }_i$ for $j=p+1,\cdots ,n$ and $i=1,\cdots ,p$ with the restriction $j+i\le n$, and zeros elsewhere. The fast calculation of $P\left(\theta \right)Y$ and $P\left(\theta \right)X$ may rely on the convolution procedure, except for the first p rows where it has to be implemented manually. The typical element of $V_p^{-1}$, $v_{ij}$, is given in (Hamilton 1994, p. 125) or (Galbraith and Galbraith 1974, p. 70). For the converged $\beta$, $\tilde{\beta }$, $co{v}_{\tilde{\beta }}={\tilde{\sigma }}_v^2{(X^{\prime }P{({\tilde{\theta }}_{PW})}^{\prime }P({\tilde{\theta }}_{PW})X)}^{-1}$ is used, where $\tilde{v}$ is the converged innovation residual and ${\tilde{\sigma }}_v^2={\tilde{v}}^{\prime }\tilde{v}/(n-k)$. Assuming normality, we can calculate the quasi-ML covariance $co{v}_{{\tilde{\theta }}_{PW}}$ from the inverse of minus the Hessian of the concentrated unconditional loglikelihood $L\left(\theta \right)=-(n/2)ln{(Y^{*}-X^{*}{\widehat{\beta }}_{GLS,\theta })}^{\prime }(Y^{*}-X^{*}{\widehat{\beta }}_{GLS,\theta })$ evaluated at converged ${\tilde{\theta }}_{PW}$ with ${\widehat{\beta }}_{GLS,\theta }={(X^{*}{}^{\prime }{X}^{*})}^{-1}(X^{*}{}^{\prime }{Y}^{*})$, $Y^{*}=P\left(\theta \right)Y$ and $X^{*}=P\left(\theta \right)X$, and the MLE innovation variance estimate ${\widehat{\sigma }}_{vMLE,\theta }^2=\widehat{v}{\left(\theta \right)}^{\prime }\widehat{v}\left(\theta \right)/n$, with “residual” $\widehat{v}\left(\theta \right)$ for a given vector $\theta$, $\widehat{v}\left(\theta \right)=Y-X{\widehat{\beta }}_{GLS,\theta }$. An alternative covariance for ${\tilde{\theta }}_{PW}$ is the asymptotic covariance ${\tilde{V}}_p^{-1}/n$ (${\tilde{V}}_p^{-1}$ is $V_p^{-1}$ evaluated at ${\tilde{\theta }}_{PW}={\{{\tilde{\theta }}_{1PW},\cdots ,{\tilde{\theta }}_{pPW}\}}^{\prime }$). A third covariance option is the sandwich covariance. Again, all these covariances are not affected by the presence of the lagged dependent variable as a regressor.

GAUSS™/MATLAB™/GNU Octave programmes for the method are available upon request.

4. Conclusions

We have proposed the correct, fast, and generic PW algorithm for exact NLS estimation of a regression with AR errors of any order (requiring no normality assumption). However, if the innovations are normal, then our proposed estimators are identical to the unconditional ML estimators. No numerical optimiser is required either, although we have checked the results of our algorithm with the results of numerical optimisation. Quasi-normality may be assumed in constructing a covariance (quasi-ML or sandwich) for the AR parameters. Alternatively, the asymptotic covariance can also be calculated (with no normality requirement). The covariance of the regressor parameters is obtained from the GLS formulae. An information criterion may be employed to select the unknown AR order in practise, approximating any underlying ARMA process. For a recent paper on AR lag order selection, see Butler and Paolella (2017). The estimation method is suitable and efficient for trending data. It is also applicable to dynamic regression with a lagged dependent variable included as a regressor. It provides reliable standard errors in the presence of the lagged dependent variable as a regressor. In addition, it can be used to derive efficient Wald-type unit root tests.