Violations of Classical Assumptions: Heteroscedasticity, Weighted Least Squares, GLS & Autocorrelation

Consider the generalized linear regression model in matrix notation for $n$ observations and $p$ parameters (consisting of an intercept and $k = p - 1$ explanatory variables):

where $\mathbf{Y} \in \mathbb{R}^{n \times 1}$ is the vector of observed response variables, $\mathbf{X} \in \mathbb{R}^{n \times p}$ is a full-column-rank design matrix ($\operatorname{rank}(\mathbf{X}) = p < n$), $\boldsymbol{\beta} \in \mathbb{R}^{p \times 1}$ is the vector of unknown population regression parameters, and $\boldsymbol{\varepsilon} \in \mathbb{R}^{n \times 1}$ is the unobserved random disturbance vector. We maintain the strict exogeneity assumption $\mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] = \mathbf{0}$, which guarantees that the linear predictor $\mathbf{X}\boldsymbol{\beta}$ correctly specifies the conditional expectation function $\mathbb{E}[\mathbf{Y} \mid \mathbf{X}] = \mathbf{X}\boldsymbol{\beta}$. However, we relax the spherical error assumption and specify a general positive definite conditional covariance structure:

Here, $\sigma^2 > 0$ is a positive scaling constant (often normalized so that $\operatorname{tr}(\boldsymbol{\Omega}) = n$), and $\boldsymbol{\Omega}$ is a symmetric, positive definite $n \times n$ matrix ($\boldsymbol{\Omega}^\top = \boldsymbol{\Omega}$ and $\mathbf{a}^\top \boldsymbol{\Omega} \mathbf{a} > 0$ for all non-zero $\mathbf{a} \in \mathbb{R}^n$). This general parameterization encompasses three primary empirical regimes:

1. Classical Spherical Disturbances: $\boldsymbol{\Omega} = \mathbf{I}_n$, so $\boldsymbol{\Sigma} = \sigma^2 \mathbf{I}_n$. All disturbances have constant variance $\sigma^2$ and zero pairwise covariance.

2. Pure Heteroscedasticity: $\boldsymbol{\Omega} = \operatorname{diag}(\omega_1, \omega_2, \dots, \omega_n)$ with $\omega_i > 0$ not all equal. The disturbances are mutually uncorrelated ($\operatorname{Cov}(\varepsilon_i, \varepsilon_j \mid \mathbf{X}) = 0$ for $i \neq j$), but the conditional variance varies across observations as $\operatorname{Var}(\varepsilon_i \mid \mathbf{X}) = \sigma_i^2 = \sigma^2 \omega_i$.

3. Serial Correlation / Autocorrelation and Spatial Dependence: $\boldsymbol{\Omega}$ contains non-zero off-diagonal entries $\omega_{ij} \neq 0$ for $i \neq j$, reflecting temporal correlation across time-series observations or spatial correlation across geographic units.

Unbiasedness and Consistency of OLS Under Nonspherical Errors

Suppose an investigator applies Ordinary Least Squares to estimate $\boldsymbol{\beta}$ without accounting for the non-identity structure of $\boldsymbol{\Omega}$:

Taking the conditional expectation given $\mathbf{X}$, we observe that unbiasedness depends exclusively on the first-moment condition $\mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] = \mathbf{0}$ and requires no assumptions regarding the second-moment matrix $\boldsymbol{\Omega}$:

Thus, OLS remains unconditionally and conditionally unbiased under arbitrary heteroscedasticity and autocorrelation, provided that the regressors are exogenous. Furthermore, if the sample size grows ($n \to \infty$) such that the empirical Gram matrix converges in probability to a finite positive definite matrix, $\operatorname{plim}_{n \to \infty} \left(\frac{1}{n}\mathbf{X}^\top \mathbf{X}\right) = \mathbf{Q}_{XX}$, and the scaled error covariance satisfies $\lim_{n \to \infty} \frac{1}{n^2} \mathbf{X}^\top (\sigma^2 \boldsymbol{\Omega}) \mathbf{X} = \mathbf{0}$, then $\operatorname{Var}(\hat{\boldsymbol{\beta}}_{\text{OLS}} \mid \mathbf{X}) \to \mathbf{0}$ in matrix norm. By Chebyshev's inequality, $\hat{\boldsymbol{\beta}}_{\text{OLS}} \xrightarrow{P} \boldsymbol{\beta}$, confirming that OLS remains weakly consistent in large samples under mild regularity conditions on the growth of $\boldsymbol{\Omega}$.

True Sampling Variance of OLS Versus the Naive Covariance Matrix

Despite preserving unbiasedness, the sampling precision of $\hat{\boldsymbol{\beta}}_{\text{OLS}}$ is profoundly altered by $\boldsymbol{\Omega}$. We derive the exact conditional variance-covariance matrix of the OLS estimator by evaluating the outer product of the estimation error $\hat{\boldsymbol{\beta}}_{\text{OLS}} - \boldsymbol{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \boldsymbol{\varepsilon}$:

Substituting $\mathbb{E}[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\top \mid \mathbf{X}] = \sigma^2 \boldsymbol{\Omega} = \boldsymbol{\Sigma}$, we obtain the true sandwich variance-covariance matrix of OLS:

Notice that the outer factors $(\mathbf{X}^\top \mathbf{X})^{-1}$ form the "bread" of the sandwich, while the middle matrix $\mathbf{X}^\top \boldsymbol{\Sigma} \mathbf{X}$ forms the "meat." Under the classical spherical assumption ($\boldsymbol{\Omega} = \mathbf{I}_n$), the middle factor collapses to $\mathbf{X}^\top \mathbf{I}_n \mathbf{X} = \mathbf{X}^\top \mathbf{X}$, cancelling one inverse factor and reducing to the familiar formula $\sigma^2 (\mathbf{X}^\top \mathbf{X})^{-1}$. When $\boldsymbol{\Omega} \neq \mathbf{I}_n$, no such cancellation occurs, and the true variance $(\mathbf{X}^\top \mathbf{X})^{-1}(\mathbf{X}^\top \boldsymbol{\Sigma} \mathbf{X})(\mathbf{X}^\top \mathbf{X})^{-1}$ differs structurally from $\sigma^2(\mathbf{X}^\top \mathbf{X})^{-1}$.

Mathematical Proof of the Bias of the Classical Residual Variance Estimator $s^2$

Even more subtle is the behavior of the classical residual variance estimator $s^2 = \frac{\mathbf{e}^\top \mathbf{e}}{n - p}$, where $\mathbf{e} = \mathbf{Y} - \mathbf{X}\hat{\boldsymbol{\beta}}_{\text{OLS}} = (\mathbf{I}_n - \mathbf{H})\boldsymbol{\varepsilon}$ is the OLS residual vector and $\mathbf{H} = \mathbf{X}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top$ is the orthogonal projection (hat) matrix. Taking the expectation of the Residual Sum of Squares ($\text{SSE} = \mathbf{e}^\top \mathbf{e}$) using the cyclic permutation invariance of the matrix trace operator $\operatorname{tr}(\mathbf{A}\mathbf{B}) = \operatorname{tr}(\mathbf{B}\mathbf{A})$, we obtain:

Expanding the trace yields:

Consequently, the naive OLS covariance matrix estimator $\widehat{\operatorname{Var}}_{\text{naive}}(\hat{\boldsymbol{\beta}}_{\text{OLS}}) = s^2 (\mathbf{X}^\top \mathbf{X})^{-1}$ has expected value:

To appreciate the practical direction of this bias, consider simple linear regression $Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$ under pure heteroscedasticity ($\operatorname{Var}(\varepsilon_i) = \sigma_i^2$, $\operatorname{Cov}(\varepsilon_i, \varepsilon_j) = 0$). The true sampling variance of the OLS slope estimator $\hat{\beta}_{1,\text{OLS}} = \frac{\sum_{i=1}^n (x_i - \bar{x})Y_i}{S_{xx}}$ (where $S_{xx} = \sum_{i=1}^n (x_i - \bar{x})^2$) is:

Meanwhile, in large samples, the naive estimator $s^2 / S_{xx}$ converges approximately to $\frac{\bar{\sigma}^2}{S_{xx}}$, where $\bar{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n \sigma_i^2$ is the average error variance. Subtracting the approximate naive variance from the true variance reveals a covariance relationship between squared predictor deviations and error variances:

Whenever observations with extreme predictor values (large $(x_i - \bar{x})^2$, i.e., high leverage $h_{ii}$) also exhibit high error variance $\sigma_i^2$—the standard empirical case when variance increases with $x_i$—the numerator is strictly positive. In this prevalent scenario, the naive OLS standard error systematically underestimates the true standard error of $\hat{\beta}_{1,\text{OLS}}$, producing artificially inflated $t$-statistics and confidence intervals whose actual coverage probability falls substantially below $1 - \alpha$.

---

When the positive definite error covariance structure $\boldsymbol{\Omega}$ is known, Alexander Aitken (1935) demonstrated that optimal estimation can be fully restored by transforming the regression system into an equivalent homoscedastic, uncorrelated representation—a procedure known in signal processing and econometrics as pre-whitening.

Matrix Spectral and Cholesky Factorization of $\boldsymbol{\Omega}$

Because $\boldsymbol{\Omega} \in \mathbb{R}^{n \times n}$ is symmetric and strictly positive definite, all of its eigenvalues $\lambda_1, \dots, \lambda_n$ are real and strictly positive, and its inverse $\boldsymbol{\Omega}^{-1}$ is likewise symmetric and positive definite. By the spectral decomposition theorem (or via Cholesky factorization), there exists a non-singular $n \times n$ matrix $\mathbf{P}$ (often called the whitening matrix or symmetric square root $\boldsymbol{\Omega}^{-1/2}$) such that:

Pre-multiplying the structural regression equation $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ on the left by $\mathbf{P}$ yields the transformed linear model:

where $\mathbf{Y}^* = \mathbf{P}\mathbf{Y}$, $\mathbf{X}^* = \mathbf{P}\mathbf{X}$, and $\boldsymbol{\varepsilon}^* = \mathbf{P}\boldsymbol{\varepsilon}$. Let us verify that the transformed model satisfies all classical Gauss–Markov assumptions:

1. Full Column Rank: Since $\mathbf{P}$ is non-singular ($\operatorname{rank}(\mathbf{P}) = n$) and $\operatorname{rank}(\mathbf{X}) = p$, Sylvester's rank inequality guarantees $\operatorname{rank}(\mathbf{X}^*) = \operatorname{rank}(\mathbf{P}\mathbf{X}) = p$.

2. Zero Conditional Mean: $\mathbb{E}[\boldsymbol{\varepsilon}^* \mid \mathbf{X}] = \mathbf{P}\mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] = \mathbf{P}\mathbf{0} = \mathbf{0}$.

3. Spherical Transformed Error Covariance: Evaluating the conditional covariance matrix of $\boldsymbol{\varepsilon}^*$ gives:

Derivation of the Generalized Least Squares (GLS) Estimator

Because the transformed model $\mathbf{Y}^* = \mathbf{X}^*\boldsymbol{\beta} + \boldsymbol{\varepsilon}^*$ satisfies the classical spherical assumptions, we obtain the optimal estimator of $\boldsymbol{\beta}$ by applying Ordinary Least Squares to the transformed quantities $(\mathbf{Y}^*, \mathbf{X}^*)$. Minimizing the transformed sum of squared residuals corresponds to minimizing the generalized quadratic form in the Mahalanobis metric induced by $\boldsymbol{\Omega}^{-1}$:

Differentiating $S_{\text{GLS}}(\boldsymbol{\beta})$ with respect to $\boldsymbol{\beta}$ and setting the gradient vector to zero yields the Generalized Normal Equations:

Since $\mathbf{X}$ has full column rank $p$ and $\boldsymbol{\Omega}^{-1}$ is strictly positive definite, for any non-zero vector $\mathbf{c} \in \mathbb{R}^p$ the vector $\mathbf{z} = \mathbf{X}\mathbf{c}$ is non-zero, so $\mathbf{c}^\top (\mathbf{X}^\top \boldsymbol{\Omega}^{-1}\mathbf{X})\mathbf{c} = \mathbf{z}^\top \boldsymbol{\Omega}^{-1}\mathbf{z} > 0$. Thus $\mathbf{X}^\top \boldsymbol{\Omega}^{-1}\mathbf{X}$ is positive definite and invertible. Solving for $\hat{\boldsymbol{\beta}}_{\text{GLS}}$ gives Aitken's Generalized Least Squares (GLS) Estimator:

The variance-covariance matrix of $\hat{\boldsymbol{\beta}}_{\text{GLS}}$ follows directly from substituting $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$:

Moreover, the unbiased estimator of the scalar scale parameter $\sigma^2$ is the normalized generalized residual sum of squares:

Theorem and Proof: The Aitken–Gauss–Markov Theorem

[!IMPORTANT] Theorem (Aitken, 1935): In the generalized linear model $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ with $\mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] = \mathbf{0}$ and $\operatorname{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \sigma^2 \boldsymbol{\Omega}$ where $\boldsymbol{\Omega}$ is a known positive definite matrix and $\operatorname{rank}(\mathbf{X}) = p$, the GLS estimator $\hat{\boldsymbol{\beta}}_{\text{GLS}} = (\mathbf{X}^\top \boldsymbol{\Omega}^{-1}\mathbf{X})^{-1}\mathbf{X}^\top \boldsymbol{\Omega}^{-1}\mathbf{Y}$ is the unique Best Linear Unbiased Estimator (BLUE) of $\boldsymbol{\beta}$. That is, for any other linear unbiased estimator $\tilde{\boldsymbol{\beta}} = \mathbf{C}\mathbf{Y}$, the covariance matrix difference $\operatorname{Var}(\tilde{\boldsymbol{\beta}} \mid \mathbf{X}) - \operatorname{Var}(\hat{\boldsymbol{\beta}}_{\text{GLS}} \mid \mathbf{X})$ is positive semi-definite.

Proof: Let $\tilde{\boldsymbol{\beta}} = \mathbf{C}\mathbf{Y}$ be an arbitrary linear estimator of $\boldsymbol{\beta}$, where $\mathbf{C} \in \mathbb{R}^{p \times n}$ is a non-stochastic matrix (or measurable with respect to $\mathbf{X}$). Define the discrepancy matrix $\mathbf{D} = \mathbf{C} - (\mathbf{X}^\top \boldsymbol{\Omega}^{-1}\mathbf{X})^{-1}\mathbf{X}^\top \boldsymbol{\Omega}^{-1}$, so that $\mathbf{C} = (\mathbf{X}^\top \boldsymbol{\Omega}^{-1}\mathbf{X})^{-1}\mathbf{X}^\top \boldsymbol{\Omega}^{-1} + \mathbf{D}$. Taking the conditional expectation of $\tilde{\boldsymbol{\beta}}$ gives:

For $\tilde{\boldsymbol{\beta}}$ to be unbiased for all parameter vectors $\boldsymbol{\beta} \in \mathbb{R}^p$, it is necessary and sufficient that $\mathbf{D}\mathbf{X} = \mathbf{0}_{p \times p}$. Now compute the variance-covariance matrix of $\tilde{\boldsymbol{\beta}}$:

Expanding the quadratic product yields four terms:

Observe that the two cross-product terms simplify via $\boldsymbol{\Omega}^{-1}\boldsymbol{\Omega} = \mathbf{I}_n$ and vanish identically by the unbiasedness orthogonality condition $\mathbf{D}\mathbf{X} = \mathbf{0}$ (and $\mathbf{X}^\top \mathbf{D}^\top = (\mathbf{D}\mathbf{X})^\top = \mathbf{0}$):

Because $\boldsymbol{\Omega}$ is positive definite, for any vector $\mathbf{a} \in \mathbb{R}^p$, the quadratic form $\mathbf{a}^\top (\sigma^2 \mathbf{D}\boldsymbol{\Omega}\mathbf{D}^\top)\mathbf{a} = \sigma^2 (\mathbf{D}^\top \mathbf{a})^\top \boldsymbol{\Omega} (\mathbf{D}^\top \mathbf{a}) \ge 0$, proving that $\operatorname{Var}(\tilde{\boldsymbol{\beta}} \mid \mathbf{X}) - \operatorname{Var}(\hat{\boldsymbol{\beta}}_{\text{GLS}} \mid \mathbf{X}) = \sigma^2 \mathbf{D}\boldsymbol{\Omega}\mathbf{D}^\top$ is positive semi-definite. In particular, selecting $\tilde{\boldsymbol{\beta}} = \hat{\boldsymbol{\beta}}_{\text{OLS}}$ reveals that $\operatorname{Var}(\hat{\boldsymbol{\beta}}_{\text{OLS}} \mid \mathbf{X}) - \operatorname{Var}(\hat{\boldsymbol{\beta}}_{\text{GLS}} \mid \mathbf{X})$ is positive semi-definite (and strictly positive definite unless the column space of $\boldsymbol{\Omega}\mathbf{X}$ coincides with the column space of $\mathbf{X}$, known as Zyskind's condition or Kruskal's theorem). $\blacksquare$

---

An important specialization of Generalized Least Squares arises when the disturbances are mutually uncorrelated across observations ($\operatorname{Cov}(\varepsilon_i, \varepsilon_j \mid \mathbf{X}) = 0$ for $i \neq j$), but exhibit non-constant conditional variance $\operatorname{Var}(\varepsilon_i \mid \mathbf{X}) = \sigma_i^2 = \sigma^2 \omega_i$. In this pure heteroscedasticity regime, the covariance matrix $\boldsymbol{\Omega}$ and its inverse $\boldsymbol{\Omega}^{-1}$ are diagonal matrices:

The symmetric whitening matrix is simply the diagonal matrix of inverse standard deviations, $\mathbf{P} = \mathbf{W}^{1/2} = \operatorname{diag}(\sqrt{w_1}, \sqrt{w_2}, \dots, \sqrt{w_n})$. Applying $\mathbf{P}$ to the $i$-th row of the regression system $Y_i = \beta_0 + \beta_1 X_{i1} + \dots + \beta_{p-1} X_{i,p-1} + \varepsilon_i$ divides every term in the $i$-th observation (including the intercept column of ones!) by $\sqrt{\omega_i}$:

Each transformed error $\varepsilon_i^* = \varepsilon_i / \sqrt{\omega_i}$ now has constant variance $\operatorname{Var}(\varepsilon_i^* \mid \mathbf{X}) = \frac{1}{\omega_i}\operatorname{Var}(\varepsilon_i \mid \mathbf{X}) = \frac{\sigma^2 \omega_i}{\omega_i} = \sigma^2$. Notice that the transformed equation has no constant intercept column unless one of the original regressors was already proportional to $\sqrt{\omega_i}$; instead, $\beta_0$ is the slope coefficient attached to the reciprocal square-root weight variable $1/\sqrt{\omega_i}$.

The WLS Objective Function and Normal Equations

The GLS criterion under diagonal $\boldsymbol{\Omega}$ reduces to minimizing the weighted sum of squared residuals, where each squared residual is weighted by $w_i = 1/\omega_i \propto 1/\sigma_i^2$:

This yields the Weighted Least Squares (WLS) Estimator:

Intuitively, WLS assigns high analytical weight $w_i$ to observations drawn from low-variance, high-precision sub-populations (small $\sigma_i^2$), forcing the fitted hyperplane to pass closely through reliable data points, while down-weighting noisy, volatile observations (large $\sigma_i^2$) so that erratic outliers do not pull the regression surface off course.

In the special case of Simple Linear Regression $Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$ with weights $w_i > 0$, minimizing $\sum_{i=1}^n w_i (Y_i - \beta_0 - \beta_1 x_i)^2$ yields closed-form expressions that parallel classical OLS except that unweighted sample means are replaced by weighted sample means $\bar{x}_w = \frac{\sum_{i=1}^n w_i x_i}{\sum_{i=1}^n w_i}$ and $\bar{Y}_w = \frac{\sum_{i=1}^n w_i Y_i}{\sum_{i=1}^n w_i}$:

with exact sampling variance $\operatorname{Var}(\hat{\beta}_{1,\text{WLS}}) = \frac{\sigma^2}{\sum_{i=1}^n w_i (x_i - \bar{x}_w)^2}$. By the Cauchy–Schwarz inequality, for any positive sequence $\omega_i$:

which immediately implies $\operatorname{Var}(\hat{\beta}_{1,\text{WLS}}) \le \operatorname{Var}(\hat{\beta}_{1,\text{OLS}})$, with strict inequality whenever $\omega_i$ varies across observations with $x_i \neq \bar{x}$.

Common Empirical Sources of Known vs. Unknown Weights

1. Grouped / Cell-Mean Aggregated Data: Suppose micro-level observations $Y_{ij} = \mathbf{x}_i^\top \boldsymbol{\beta} + \varepsilon_{ij}$ ($j = 1, \dots, n_i$) within region or demographic cell $i$ satisfy homoscedastic individual variance $\operatorname{Var}(\varepsilon_{ij}) = \sigma^2$, but the investigator only observes cell averages $\bar{Y}_i = \frac{1}{n_i}\sum_{j=1}^{n_i} Y_{ij}$ based on group sizes $n_i$. Then $\bar{Y}_i = \mathbf{x}_i^\top \boldsymbol{\beta} + \bar{\varepsilon}_i$, where the averaged error has variance $\operatorname{Var}(\bar{\varepsilon}_i) = \frac{\sigma^2}{n_i}$. Here $\omega_i = 1/n_i$ is known exactly from the sample sizes, so WLS uses weights $w_i = n_i$, multiplying each cell-mean equation by $\sqrt{n_i}$.

2. Proportional Variance Models: In cross-sectional expenditure models, standard deviation often scales proportionally with a size variable $Z_i$ (e.g., firm assets or household income): $\operatorname{SD}(\varepsilon_i \mid Z_i) = \sigma Z_i$, so $\operatorname{Var}(\varepsilon_i \mid Z_i) = \sigma^2 Z_i^2$. Here $\omega_i = Z_i^2$ and WLS weights are $w_i = 1/Z_i^2$, dividing the regression equation through by $Z_i$.

Feasible Generalized Least Squares (FGLS / Two-Step Estimation)

In most econometric applications, the covariance matrix $\boldsymbol{\Omega}$ is not known a priori. If $\boldsymbol{\Omega}$ contained $n(n+1)/2$ unrestricted unknown parameters, total parameters would exceed the sample size $n$, rendering estimation impossible (an incidental parameters problem). Instead, Feasible Generalized Least Squares (FGLS) specifies $\boldsymbol{\Omega} = \boldsymbol{\Omega}(\boldsymbol{\gamma})$ as a smooth parametric function of a low-dimensional parameter vector $\boldsymbol{\gamma} \in \mathbb{R}^m$ ($m \ll n$). A standard specification for heteroscedasticity is the multiplicative exponential skedastic function (Harvey, 1976):

where $\mathbf{z}_i \in \mathbb{R}^m$ is a vector of exogenous variance-driving variables (often a subset of $\mathbf{x}_i$), and $\sigma^2 = e^{\alpha_0}$. The exponential form guarantees strictly positive predicted variances $\hat{\sigma}_i^2 > 0$ for any parameter vector $\hat{\boldsymbol{\gamma}}$. The two-step FGLS procedure proceeds as follows:

• Step 1 (Consistent Initial OLS and Residual Extraction): Estimate $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ by OLS to obtain consistent residuals $e_i = Y_i - \mathbf{x}_i^\top \hat{\boldsymbol{\beta}}_{\text{OLS}}$.
• Step 2 (Auxiliary Skedastic Regression): Take the natural logarithm of the squared OLS residuals $\ln(e_i^2) = \alpha_0 + \mathbf{z}_i^\top \boldsymbol{\gamma} + v_i$, where $v_i = \ln(e_i^2 / \sigma_i^2)$, and estimate $(\hat{\alpha}_0, \hat{\boldsymbol{\gamma}})$ via OLS. Compute fitted variance weights $\hat{\omega}_i = \exp(\mathbf{z}_i^\top \hat{\boldsymbol{\gamma}})$ and form $\hat{\boldsymbol{\Omega}} = \operatorname{diag}(\hat{\omega}_1, \dots, \hat{\omega}_n)$.
• Step 3 (Feasible WLS Update): Compute the FGLS estimator by replacing $\boldsymbol{\Omega}$ with its consistent plug-in estimate $\hat{\boldsymbol{\Omega}}$:

Under standard regularity conditions—specifically, if $\hat{\boldsymbol{\gamma}} \xrightarrow{P} \boldsymbol{\gamma}$ such that $\operatorname{plim}_{n \to \infty} \frac{1}{n}\mathbf{X}^\top (\hat{\boldsymbol{\Omega}}^{-1} - \boldsymbol{\Omega}^{-1})\mathbf{X} = \mathbf{0}$ and $\operatorname{plim}_{n \to \infty} \frac{1}{\sqrt{n}}\mathbf{X}^\top (\hat{\boldsymbol{\Omega}}^{-1} - \boldsymbol{\Omega}^{-1})\boldsymbol{\varepsilon} = \mathbf{0}$—the estimation error in $\hat{\boldsymbol{\Omega}}$ vanishes asymptotically at a faster rate than the $\sqrt{n}$-scaled regression score. Consequently, $\hat{\boldsymbol{\beta}}_{\text{FGLS}}$ is asymptotically equivalent to the infeasible exact GLS estimator:

[!WARNING] Finite-Sample Caveat of FGLS: Although $\hat{\boldsymbol{\beta}}_{\text{FGLS}}$ is asymptotically efficient when the skedastic function $\sigma_i^2 = \sigma^2 \omega(\mathbf{z}_i; \boldsymbol{\gamma})$ is correctly specified, $\hat{\boldsymbol{\beta}}_{\text{FGLS}}$ is a non-linear function of $\mathbf{Y}$ (because $\hat{\boldsymbol{\Omega}}$ depends on the data $\mathbf{Y}$ via the first-stage OLS residuals). Consequently, in finite samples, FGLS is generally neither exactly unbiased nor guaranteed to have smaller variance than OLS if $\hat{\boldsymbol{\gamma}}$ is noisily estimated.

---

Before abandoning classical OLS inference or specifying an FGLS variance model, an empirical researcher must test the null hypothesis of homoscedasticity $H_0: \sigma_1^2 = \sigma_2^2 = \dots = \sigma_n^2 = \sigma^2$ against structured or unstructured heteroscedastic alternatives. Because the true disturbances $\varepsilon_i$ are unobserved, all diagnostic tests use the OLS residuals $e_i = Y_i - \mathbf{x}_i^\top \hat{\boldsymbol{\beta}}_{\text{OLS}}$ as empirical proxies for $\varepsilon_i$. ### 1. The Goldfeld–Quandt Exact Finite-Sample $F$-Test (1965) When error variance is suspected to be a monotonic function of a single known continuous explanatory variable $X_j$ (e.g., $\sigma_i^2 = \sigma^2 X_{ij}^2$) and the disturbances are normally distributed ($\boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma})$), Stephen Goldfeld and Richard Quandt constructed an exact finite-sample $F$-test by partitioning the sample into low-variance and high-variance regimes: 1. **Order Observations**: Sort the $n$ observations in ascending order of the suspect variable $X_{ij}$. 2. **Omit Central Observations**: Discard $c$ central observations (typically $c \approx n/5$ to $n/4$) to sharpen the variance contrast between the lower and upper extremes, leaving two subsamples of size $n_1 = \frac{n - c}{2}$ (low $X_j$) and $n_2 = \frac{n - c}{2}$ (high $X_j$), requiring $n_1 > p$. 3. **Separate OLS Regressions**: Fit the linear regression model separately to the first $n_1$ observations and the last $n_2$ observations, obtaining residual sums of squares $\text{SSE}_1$ (with $n_1 - p$ degrees of freedom) and $\text{SSE}_2$ (with $n_2 - p$ degrees of freedom). 4. **Test Statistic**: Under $H_0: \sigma_i^2 = \sigma^2$ and normal independent errors, the two subsamples are completely independent, so $\frac{\text{SSE}_1}{\sigma^2} \sim \chi^2_{n_1 - p}$ and $\frac{\text{SSE}_2}{\sigma^2} \sim \chi^2_{n_2 - p}$ independently. The Goldfeld–Quandt statistic follows an exact Snedecor $F$-distribution under $H_0$: $$F_{\text{GQ}} = \frac{\text{SSE}_2 / (n_2 - p)}{\text{SSE}_1 / (n_1 - p)} = \frac{s_2^2}{s_1^2} \;\underset{H_0}{\sim}\; F_{n_2 - p, \, n_1 - p}$$ If error variance increases with $X_j$, $s_2^2$ will systematically exceed $s_1^2$, and we reject $H_0$ at significance level $\alpha$ if $F_{\text{GQ}} > F_{\alpha, \, n_2 - p, \, n_1 - p}$. ### 2. The Breusch–Pagan–Godfrey Lagrange Multiplier Test and Koenker–Bassett Studentization To test whether error variance depends on a vector of $m$ explanatory variables $\mathbf{z}_i = (z_{i1}, \dots, z_{im})^\top$ via any twice-differentiable skedastic function $\sigma_i^2 = h(\alpha_0 + \mathbf{z}_i^\top \boldsymbol{\gamma})$ where $h(\cdot) > 0$ and $h'(\cdot) \neq 0$, Trevor Breusch and Adrian Pagan (1979) and Leslie Godfrey (1978) derived the Lagrange Multiplier (Score) test of $H_0: \boldsymbol{\gamma} = \mathbf{0}$ (since $\boldsymbol{\gamma} = \mathbf{0} \implies \sigma_i^2 = h(\alpha_0) = \sigma^2$). - **Original Gaussian Breusch–Pagan Test**: Under normal errors $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$, the fourth central moment is $\mathbb{E}[\varepsilon_i^4] = 3\sigma^4$, so the variance of the squared error is $\operatorname{Var}(\varepsilon_i^2) = \mathbb{E}[\varepsilon_i^4] - (\mathbb{E}[\varepsilon_i^2])^2 = 3\sigma^4 - \sigma^4 = 2\sigma^4$, which implies $\operatorname{Var}(\varepsilon_i^2 / \sigma^2) = 2$. Defining the maximum likelihood variance estimate $\hat{\sigma}_{\text{MLE}}^2 = \frac{1}{n}\sum_{i=1}^n e_i^2$ and standardized squared residuals $g_i = e_i^2 / \hat{\sigma}_{\text{MLE}}^2$, the Gaussian score statistic equals one-half the Explained Sum of Squares ($\text{SSR}_{\text{aux}}$) from regressing $g_i$ on $(1, \mathbf{z}_i^\top)$: $$\text{LM}_{\text{BP-Gaussian}} = \frac{1}{2} \text{SSR}_{\text{aux}} \;\underset{H_0}{\xrightarrow{d}}\; \chi^2_m$$ - **Koenker–Bassett (1982) Studentized Breusch–Pagan Test**: Roger Koenker observed that the divisor $2\sigma^4$ relies critically on Gaussian kurtosis ($\kappa = 3$). If the true errors are heavy-tailed (leptokurtic, $\kappa > 3$), $\operatorname{Var}(\varepsilon_i^2) = (\kappa - 1)\sigma^4 > 2\sigma^4$, causing the original Breusch–Pagan test to reject homoscedasticity far too frequently merely because of non-normality! Replacing the theoretical Gaussian variance $2\sigma^4$ with the consistent empirical variance of the squared residuals $\hat{V} = \frac{1}{n}\sum_{i=1}^n (e_i^2 - \hat{\sigma}_{\text{MLE}}^2)^2$ yields the **Koenker–Bassett studentized LM statistic**, which equals the sample size $n$ multiplied by the unadjusted coefficient of determination $R_{\text{aux}}^2$ from the auxiliary OLS regression of $e_i^2$ on an intercept and $z_{i1}, \dots, z_{im}$: $$e_i^2 = \delta_0 + \delta_1 z_{i1} + \dots + \delta_m z_{im} + u_i \quad \Longrightarrow \quad \text{LM}_{\text{BP}} = n R_{\text{aux}}^2 \;\underset{H_0}{\xrightarrow{d}}\; \chi^2_m$$ This $n R_{\text{aux}}^2$ formulation is completely robust to non-normal kurtosis and is the standard version implemented in modern statistical software. ### 3. Park (1966) and Glejser (1969) Auxiliary Residual Tests Two classical precursor tests estimate functional forms directly: - **Park Test**: Regresses $\ln(e_i^2) = \alpha_0 + \gamma \ln(X_i) + v_i$ and performs a $t$-test on $\hat{\gamma}$. - **Glejser Test**: Regresses the absolute residuals $|e_i|$ on candidate functional forms of $X_i$, such as $|e_i| = \delta_0 + \delta_1 X_i + u_i$, $|e_i| = \delta_0 + \delta_1 \sqrt{X_i} + u_i$, or $|e_i| = \delta_0 + \delta_1 (1/X_i) + u_i$. This helps identify the specific power-law relationship needed to construct FGLS weights. ### 4. White's General Heteroscedasticity and Specification Test (1980) Halbert White (1980) recognized that heteroscedasticity compromises OLS inference **only** insofar as the middle matrix $\frac{1}{n}\mathbf{X}^\top \boldsymbol{\Sigma} \mathbf{X} = \frac{1}{n}\sum_{i=1}^n \sigma_i^2 \mathbf{x}_i \mathbf{x}_i^\top$ diverges from $\sigma^2 \left(\frac{1}{n}\sum_{i=1}^n \mathbf{x}_i \mathbf{x}_i^\top\right)$. For the difference between these two $p \times p$ matrices to vanish, $\sigma_i^2$ must be uncorrelated with every unique element of the outer product matrix $\mathbf{x}_i \mathbf{x}_i^\top$—namely, all linear terms $X_{ij}$, all pure quadratic terms $X_{ij}^2$, and all cross-product interaction terms $X_{ij} X_{i\ell}$ ($1 \le j < \ell \le k$). Accordingly, **White's Test** requires no prior knowledge of which variables drive the variance: 1. Estimate the primary regression by OLS and compute the squared residuals $e_i^2$. 2. Run an auxiliary OLS regression of $e_i^2$ on an intercept, all $k$ original regressors, all $k$ squared regressors, and all $\binom{k}{2} = \frac{k(k-1)}{2}$ pairwise cross-products (omitting any redundant columns, such as the square of a binary $0/1$ dummy variable since $D_i^2 = D_i$): $$e_i^2 = \delta_0 + \sum_{j=1}^k \delta_j X_{ij} + \sum_{j=1}^k \gamma_{jj} X_{ij}^2 + \sum_{j=1}^{k-1}\sum_{\ell=j+1}^k \gamma_{j\ell} X_{ij} X_{i\ell} + u_i$$ 3. With $q$ non-redundant auxiliary regressors (where $q = 2k + \frac{k(k-1)}{2} = \frac{k(k+3)}{2}$ if all $k$ predictors are continuous), compute the White test statistic: $$W = n R_{\text{White}}^2 \;\underset{H_0}{\xrightarrow{d}}\; \chi^2_q$$ Reject $H_0$ if $W > \chi^2_{\alpha, \, q}$. Note that because $e_i^2$ also captures omitted non-linear terms in $\mathbb{E}[Y_i \mid \mathbf{x}_i]$, a significant White statistic can signal either pure heteroscedasticity or functional form misspecification (omitted quadratic/interaction terms in the mean equation). ---

Suppose diagnostic tests reject homoscedasticity, or theoretical considerations dictate that $\sigma_i^2$ varies across observations, yet the true functional form of $\sigma_i^2$ is unknown. Rather than guessing a parametric skedastic model for FGLS (which risks misspecification), we can retain the unbiased and consistent OLS point estimator $\hat{\boldsymbol{\beta}}_{\text{OLS}}$ and construct a consistent estimator of its true sampling variance under arbitrary, unknown heteroscedasticity.

Recall the true variance of the scaled OLS estimation error $\sqrt{n}(\hat{\boldsymbol{\beta}}_{\text{OLS}} - \boldsymbol{\beta})$:

Under pure heteroscedasticity, $\boldsymbol{\Sigma} = \operatorname{diag}(\sigma_1^2, \sigma_2^2, \dots, \sigma_n^2)$, so the middle $p \times p$ meat matrix expands as a weighted sum of outer products of the regressor rows $\mathbf{x}_i \in \mathbb{R}^{p \times 1}$:

Wait: can we estimate the $n$ unknown individual variances $\sigma_1^2, \dots, \sigma_n^2$ consistently? No, because there is only one observation per $\sigma_i^2$. However, the fundamental insight of Friedhelm Eicker (1963), Peter Huber (1967), and Halbert White (1980) is that we do not need consistent estimates of each individual $\sigma_i^2$; we only need a consistent estimator of the $p \times p$ average matrix $\boldsymbol{\Psi}_n = \frac{1}{n}\sum_{i=1}^n \sigma_i^2 \mathbf{x}_i \mathbf{x}_i^\top$!

Replacing each unobserved $\sigma_i^2 = \mathbb{E}[\varepsilon_i^2 \mid \mathbf{X}]$ with its empirical counterpart—the squared OLS residual $e_i^2 = (Y_i - \mathbf{x}_i^\top \hat{\boldsymbol{\beta}}_{\text{OLS}})^2$—yields the empirical meat matrix $\hat{\boldsymbol{\Psi}}_n = \frac{1}{n}\sum_{i=1}^n e_i^2 \mathbf{x}_i \mathbf{x}_i^\top$. Because $e_i = \varepsilon_i - \mathbf{x}_i^\top (\hat{\boldsymbol{\beta}}_{\text{OLS}} - \boldsymbol{\beta})$ and $\hat{\boldsymbol{\beta}}_{\text{OLS}} - \boldsymbol{\beta} = O_p(n^{-1/2})$, expanding $e_i^2 = \varepsilon_i^2 - 2\varepsilon_i \mathbf{x}_i^\top (\hat{\boldsymbol{\beta}}_{\text{OLS}} - \boldsymbol{\beta}) + \left[\mathbf{x}_i^\top (\hat{\boldsymbol{\beta}}_{\text{OLS}} - \boldsymbol{\beta})\right]^2$ and applying Markov's Weak Law of Large Numbers for independent non-identically distributed random variables (under bounded fourth moments $\sup_i \mathbb{E}[\varepsilon_i^4] < \infty$ and $\sup_i \|\mathbf{x}_i\|^4 < \infty$) proves that the error averages out across the $n$ summands:

Finite-Sample Bias Corrections: The HC0, HC1, HC2, HC3, and HC4 Hierarchy

Although White's baseline estimator (HC0) is asymptotically consistent, in finite samples the squared OLS residual $e_i^2$ systematically underestimates $\sigma_i^2$ due to leverage shrinkage. Specifically, since $\mathbf{e} = (\mathbf{I}_n - \mathbf{H})\boldsymbol{\varepsilon}$, under homoscedasticity $\mathbb{E}[e_i^2 \mid \mathbf{X}] = (1 - h_{ii})\sigma^2 < \sigma^2$, where $h_{ii} = \mathbf{x}_i^\top (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{x}_i \in [1/n, 1]$ is the $i$-th diagonal of the hat matrix. When an observation has high leverage ($h_{ii}$ close to $1$), the OLS regression line is pulled toward $Y_i$, shrinking $e_i^2$ toward zero and biasing HC0 standard errors downward in small samples. James MacKinnon and Halbert White (1985), Russell Davidson and James MacKinnon (1993), and Francisco Cribari-Neto (2004) developed a hierarchy of leverage-adjusted sandwich estimators of the form:

where the diagonal elements $\hat{v}_i^{(k)}$ of $\hat{\boldsymbol{\Sigma}}_{\text{HC}k} = \operatorname{diag}(\hat{v}_1^{(k)}, \dots, \hat{v}_n^{(k)})$ are defined as follows:

1. HC0 (Baseline White, 1980): $\hat{v}_i^{(0)} = e_i^2$. Asymptotically valid, but downwardly biased in finite samples ($n < 250$), especially when high-leverage observations exist.

2. HC1 (Hinkley, 1977; Default in Stata `robust`): $\hat{v}_i^{(1)} = \frac{n}{n - p} e_i^2$. Applies a uniform degrees-of-freedom inflation factor $\frac{n}{n-p}$ (matching $\bar{h} = p/n$), but does not account for heterogeneous individual leverages $h_{ii}$.

3. HC2 (Horn, Horn, and Duncan, 1975): $\hat{v}_i^{(2)} = \frac{e_i^2}{1 - h_{ii}}$. Exactly unbiased for $\sigma^2$ under homoscedasticity since $\mathbb{E}\left[\frac{e_i^2}{1 - h_{ii}}\right] = \frac{(1 - h_{ii})\sigma^2}{1 - h_{ii}} = \sigma^2$.

4. HC3 (MacKinnon and White, 1985; Jackknife Approximation, Default in R `sandwich`): $\hat{v}_i^{(3)} = \frac{e_i^2}{(1 - h_{ii})^2} = e_{(i)}^2$, where $e_{(i)} = \frac{e_i}{1 - h_{ii}}$ is the leave-one-out deletion (PRESS) residual! By squaring the leverage penalty $(1 - h_{ii})^{-2}$, HC3 heavily penalizes high-leverage points and provides superior finite-sample Type I error control.

5. HC4 (Cribari-Neto, 2004): $\hat{v}_i^{(4)} = \frac{e_i^2}{(1 - h_{ii})^{\delta_i}}$, where $\delta_i = \min\left(4, \, \frac{n h_{ii}}{p}\right) = \min\left(4, \, \frac{h_{ii}}{\bar{h}}\right)$ scales the exponent up to $4$ for extreme leverage outliers ($h_{ii} > 2\bar{h}$).

---

When regression data are collected sequentially over time ($t = 1, 2, \dots, T$), the disturbance $\varepsilon_t$ represents the cumulative influence of omitted macroeconomic shocks, supply-chain bottlenecks, consumer habit persistence, and measurement conventions that persist across multiple periods. The canonical benchmark model for temporal correlation is the stationary First-Order Autoregressive Process, AR(1):

where $\rho \in (-1, 1)$ is the first-order autocorrelation coefficient and $u_t$ is a white-noise innovation sequence satisfying $\mathbb{E}[u_t \mid \mathbf{X}] = 0$, $\operatorname{Var}(u_t \mid \mathbf{X}) = \sigma_u^2$, and $\operatorname{Cov}(u_t, u_s \mid \mathbf{X}) = 0$ for $t \neq s$.

Derivation of the Stationary AR(1) Covariance Matrix $\boldsymbol{\Omega}$ and Its Analytical Inverse $\boldsymbol{\Omega}^{-1}$

By recursive backward substitution under the stationarity condition $|\rho| < 1$, we express $\varepsilon_t$ as an infinite moving-average representation of past white-noise shocks:

From this causal expansion, we compute the unconditional mean, variance, and lag-$s$ autocovariance ($s \ge 0$):

• Mean: $\mathbb{E}[\varepsilon_t] = \sum_{j=0}^\infty \rho^j \mathbb{E}[u_{t-j}] = 0$.
• Stationary Variance: Because the innovations $u_{t-j}$ are mutually uncorrelated, the variance is a convergent geometric series with ratio $\rho^2 < 1$:

• Lag-$s$ Autocovariance and Autocorrelation ($s \ge 1$): Multiplying $\varepsilon_t = \rho \varepsilon_{t-1} + u_t = \rho^s \varepsilon_{t-s} + \sum_{j=0}^{s-1} \rho^j u_{t-j}$ by $\varepsilon_{t-s}$ and taking expectations (noting that future innovations $u_{t-j}$ for $j < s$ are uncorrelated with $\varepsilon_{t-s}$) yields:

Setting the scale factor $\sigma^2 = \sigma_u^2$, the $T \times T$ covariance matrix $\operatorname{Var}(\boldsymbol{\varepsilon}) = \sigma_u^2 \boldsymbol{\Omega}$ has Toeplitz autocorrelation structure:

Remarkably, the inverse of this dense $T \times T$ matrix is an sparse, symmetric tridiagonal matrix $\boldsymbol{\Omega}^{-1} = \mathbf{P}^\top \mathbf{P}$, factored via the lower-bidiagonal Prais–Winsten transformation matrix $\mathbf{P} \in \mathbb{R}^{T \times T}$:

Let us verify why $\mathbf{P}\boldsymbol{\varepsilon}$ whitens the disturbance vector:

• For $t = 2, 3, \dots, T$, the $t$-th row of $\mathbf{P}\boldsymbol{\varepsilon}$ computes the quasi-difference $\varepsilon_t^* = \varepsilon_t - \rho \varepsilon_{t-1} = u_t$, which has variance $\sigma_u^2$ and is independent of all other rows!
• For the initial observation $t = 1$, there is no $\varepsilon_0$ in the sample. However, multiplying $\varepsilon_1$ by $\sqrt{1 - \rho^2}$ produces $\varepsilon_1^* = \sqrt{1 - \rho^2}\,\varepsilon_1$, whose variance is $(1 - \rho^2) \operatorname{Var}(\varepsilon_1) = (1 - \rho^2)\frac{\sigma_u^2}{1 - \rho^2} = \sigma_u^2$, and since $u_2, \dots, u_T$ are future innovations independent of $\varepsilon_1$, $\operatorname{Cov}(\varepsilon_1^*, \varepsilon_t^*) = 0$ for all $t \ge 2$!

Consequences of Positive Autocorrelation ($\rho > 0$) for Naive OLS Standard Errors

In economic time series, both the regressors $X_t$ and the disturbances $\varepsilon_t$ typically display positive autocorrelation ($\rho > 0$ and $\rho_X = \operatorname{Corr}(X_t, X_{t-1}) > 0$). In simple regression $Y_t = \beta_0 + \beta_1 x_t + \varepsilon_t$, the true asymptotic variance of $\hat{\beta}_{1,\text{OLS}}$ is inflated relative to the naive formula $\sigma_\varepsilon^2 / S_{xx}$ by the autocorrelation factor:

Simultaneously, positive autocorrelation causes successive disturbances to cluster on the same side of the true regression line, allowing the OLS fitted line to track the correlated error drift and artificially depressing the residual variance estimator:

Thus, when $\rho > 0$ and $\rho_X > 0$, naive OLS suffers a double downward bias: $s^2$ underestimates $\sigma_\varepsilon^2$, while $\sigma_\varepsilon^2 / S_{xx}$ omits the positive multiplier $\frac{1 + \rho \rho_X}{1 - \rho \rho_X} > 1$. For instance, if $\rho = 0.8$ and $\rho_X = 0.8$, the inflation factor is $\frac{1 + 0.64}{1 - 0.64} = \frac{1.64}{0.36} = 4.56$, meaning the true variance is more than $4.5$ times larger than the naive variance, and naive $t$-statistics are exaggerated by more than $\sqrt{4.56} \approx 2.13$-fold!

Diagnostic Tests for Serial Correlation

1. The Durbin–Watson $d$-Test (1950, 1951): To test $H_0: \rho = 0$ in an AR(1) model with strictly exogenous regressors and an intercept, James Durbin and Geoffrey Watson defined the ratio of the sum of squared first-differences of OLS residuals to the Residual Sum of Squares:

where $\hat{\rho} = \frac{\sum_{t=2}^T e_t e_{t-1}}{\sum_{t=1}^T e_t^2}$ is the sample first-order residual autocorrelation, and the approximation holds because the endpoint terms $(e_1^2 + e_T^2) / \sum_{t=1}^T e_t^2 \to 0$ as $T \to \infty$. Because $-1 \le \hat{\rho} \le 1$, the statistic satisfies $0 \le d \le 4$:

• Zero autocorrelation ($\hat{\rho} \approx 0$): $d \approx 2$.
• Strong positive autocorrelation ($\hat{\rho} \to +1$): $d \to 0$.
• Strong negative autocorrelation ($\hat{\rho} \to -1$): $d \to 4$.

Because $\mathbf{e} = (\mathbf{I}_T - \mathbf{H})\boldsymbol{\varepsilon}$, the exact finite-sample distribution of $d = \frac{\boldsymbol{\varepsilon}^\top (\mathbf{I}_T - \mathbf{H})\mathbf{A}(\mathbf{I}_T - \mathbf{H})\boldsymbol{\varepsilon}}{\boldsymbol{\varepsilon}^\top (\mathbf{I}_T - \mathbf{H})\boldsymbol{\varepsilon}}$ depends on the specific configuration of the design matrix $\mathbf{X}$. However, Durbin and Watson proved that for any design matrix $\mathbf{X}$ with $k$ regressors and an intercept, the exact critical value $d_\alpha$ is bounded between two tabulated universal values $d_L$ (lower bound) and $d_U$ (upper bound):

• If $d < d_L$: Reject $H_0: \rho = 0$ in favor of positive autocorrelation ($H_1: \rho > 0$).
• If $d_L \le d \le d_U$: The test falls in the Zone of Indecision (Inconclusive Region).
• If $d_U < d < 4 - d_U$: Fail to reject $H_0: \rho = 0$ (no evidence of AR(1) autocorrelation).
• If $4 - d_U \le d \le 4 - d_L$: Zone of Indecision for negative autocorrelation.
• If $d > 4 - d_L$: Reject $H_0$ in favor of negative autocorrelation ($H_1: \rho < 0$).

[!CAUTION] Invalidity of $d$ with Lagged Dependent Variables (Durbin's $h$-Test): When $\mathbf{X}$ contains a lagged dependent variable $Y_{t-1}$ (dynamic model $Y_t = \beta_0 + \gamma Y_{t-1} + \mathbf{x}_t^\top \boldsymbol{\beta} + \varepsilon_t$), $Y_{t-1}$ is correlated with $\varepsilon_{t-1}$, biasing $\hat{\rho}$ toward zero and biasing $d$ toward $2$. In this setting, Durbin (1970) proposed Durbin's $h$-statistic:

which is valid provided $T \cdot \widehat{\operatorname{Var}}(\hat{\gamma}) < 1$. Even more generally, the Breusch–Godfrey test remains valid regardless of whether $T \cdot \widehat{\operatorname{Var}}(\hat{\gamma}) < 1$.

2. The Breusch–Godfrey (BG) Lagrange Multiplier Test for $\text{AR}(q)$ or $\text{MA}(q)$ Errors (1978): To test $H_0: \rho_1 = \rho_2 = \dots = \rho_q = 0$ in the higher-order autoregressive process $\varepsilon_t = \rho_1 \varepsilon_{t-1} + \dots + \rho_q \varepsilon_{t-q} + u_t$ (even in the presence of lagged dependent variables $Y_{t-1}$), regress the OLS residuals $e_t$ on both the original regressors $\mathbf{x}_t$ and the $q$ lagged residuals $e_{t-1}, \dots, e_{t-q}$ (setting initial pre-sample residuals $e_0 = e_{-1} = \dots = 0$):

Including the original regressors $\mathbf{x}_t$ in this auxiliary regression accounts for the non-orthogonality between $\mathbf{x}_t$ and the lagged residuals (which is essential when $\mathbf{x}_t$ includes $Y_{t-1}$). The Breusch–Godfrey LM statistic is:

Corrective Estimation: Cochrane–Orcutt vs. Prais–Winsten FGLS and Newey–West HAC

When AR(1) serial correlation is confirmed, two classical FGLS procedures transform the model using $\hat{\rho} = \frac{\sum_{t=2}^T e_t e_{t-1}}{\sum_{t=1}^{T-1} e_t^2}$:

1. Cochrane–Orcutt (1949) Two-Step / Iterative Procedure: Drops the initial observation $t = 1$ and applies OLS to the $T - 1$ quasi-differenced observations ($t = 2, \dots, T$):

Note that the intercept estimated in the quasi-differenced regression is $\hat{\beta}_0^* = \widehat{\beta_0(1 - \rho)}$, so the original intercept is recovered via $\hat{\beta}_0 = \hat{\beta}_0^* / (1 - \hat{\rho})$. In iterative Cochrane–Orcutt, the updated parameter vector $\hat{\boldsymbol{\beta}}$ is used to recompute structural residuals $\tilde{e}_t = Y_t - \mathbf{x}_t^\top \hat{\boldsymbol{\beta}}$, update $\hat{\rho}$, and re-estimate until convergence.

2. Prais–Winsten (1954) Full-Sample FGLS: Retaining only $T - 1$ observations in Cochrane–Orcutt discards the initial state $t = 1$, which causes severe finite-sample efficiency loss when $T$ is small or when $X_t$ exhibits a strong trend. Prais and Winsten retain $t = 1$ by applying the first row of $\mathbf{P}$:

Combining this scaled $t = 1$ equation with the $T - 1$ quasi-differenced equations ($t = 2, \dots, T$) executes exact FGLS across all $T$ observations.

3. Newey–West (1987) Heteroscedasticity and Autocorrelation Consistent (HAC) Standard Errors: When the exact ARMA order is unknown, Whitney Newey and Kenneth West generalized White's sandwich estimator to time series by including Bartlett-weighted sample autocovariance matrices up to truncation lag $L$ (typically $L = \lfloor 4(T/100)^{2/9} \rfloor$):

The linear Bartlett kernel weights $w(\ell, L) = 1 - \frac{\ell}{L+1}$ guarantee that $\hat{\boldsymbol{\Psi}}_{\text{NW}}$ is strictly positive semi-definite in every finite sample while converging in probability to $\lim_{T \to \infty} \frac{1}{T}\mathbf{X}^\top \boldsymbol{\Sigma} \mathbf{X}$.

---