Introduction
Generalized least squares (GLS) is useful when regression errors are correlated across observations or have unequal variances. The central task in an applied problem is not to transform the data mechanically, but to specify a defensible covariance model and then use it consistently to estimate coefficients and quantify their uncertainty. Unequal variances, serial correlation, and shared group-level shocks imply different covariance structures and therefore different GLS procedures.
This tutorial works through those three settings. In each, the regression mean is linear, but the error covariance changes how observations should contribute to estimation. The calculations assume the stated covariance model is appropriate; in practice, its form and estimated parameters need to be assessed rather than treated as automatic.
From the covariance model to the estimator
The matrix \(\Omega^{-1}\) determines how the observations are weighted and how their dependence is accounted for. For diagonal \(\Omega\), its diagonal entries assign inverse-variance weights. For correlated errors, off-diagonal entries matter as well: the information in similar, dependent observations is not counted as if they were independent.
A useful justification comes from whitening the errors. Since \(\Omega\) is positive definite, it has a factorization \(\Omega=CC'\) with nonsingular \(C\). Premultiplying the model by \(C^{-1}\) gives \(C^{-1}y=C^{-1}X\beta+C^{-1}u\), and the transformed error has covariance \(I\). Ordinary least squares on this transformed model yields the GLS formula:
This derivation also explains the main assumption behind the procedure: the chosen transformation must represent the error covariance well enough for the transformed errors to have the intended covariance. When covariance parameters are estimated from the data, substituting their estimates gives feasible GLS (FGLS). The resulting covariance estimate is a plug-in approximation; it does not generally reproduce exact finite-sample inference that would be available if the covariance parameters were known.
Worked Example: Unequal Variances and Inverse-Variance Weighting
Worked Example: Unequal Variances and Inverse-Variance Weighting
Suppose four hypothetical observations follow \(y_i=\beta_0+\beta_1x_i+u_i\), with independent normal errors having known variances \(1,1,4,4\), respectively. The data are \((x_i,y_i)=(0,1),(1,2),(2,2),(3,4)\). Thus \(\Omega=\operatorname{diag}(1,1,4,4)\), and the inverse-variance weights are \(1,1,1/4,1/4\). The model has a linear conditional mean, full-rank design, zero conditional error means, and the stated diagonal covariance; the normality assumption supports the exact normal interval below.
The weighted normal equations use \(X'\Omega^{-1}X\) and \(X'\Omega^{-1}y\). The weighted sums are \(\sum w_i=2.5\), \(\sum w_ix_i=2.25\), \(\sum w_ix_i^2=4.25\), \(\sum w_iy_i=4.5\), and \(\sum w_ix_iy_i=6\). Therefore,
The determinant is \(2.5(4.25)-2.25^2=5.5625\). Solving the two normal equations gives
Because the error covariance is specified in absolute units, the slope variance is the lower-right entry of \((X'\Omega^{-1}X)^{-1}\), which is \(2.5/5.5625=0.4494\). Hence \(SE(\hat{\beta}_1)=\sqrt{0.4494}=0.6704\). Under the normal model with known \(\Omega\), a 95% normal interval is
The estimated slope is about 0.876 outcome units per predictor unit. The interval is a model-based 95% confidence interval for the population slope; it is not a probability statement about the fixed parameter after the interval has been computed. The observations with variance 4 receive one quarter the weight of observations with variance 1 because their errors are less precise.
Worked Example: AR(1) Errors and the Prais–Winsten Transformation
Worked Example: AR(1) Errors and the Prais–Winsten Transformation
Consider a hypothetical time-ordered regression \(y_t=\beta_0+\beta_1x_t+u_t\), for \(t=1,\ldots,4\), where \(x=(0,1,2,3)\) and \(y=(1,3,3,8)\). Suppose the errors follow a stationary Gaussian AR(1) model: \(\operatorname{Var}(u_t)=\sigma^2\) and \(\operatorname{Cov}(u_t,u_s)=\sigma^2\rho^{|t-s|}\). For this calculation, take \(\rho=0.5\) as known and assume the model is correctly ordered in time, has conditional mean zero, and has this covariance structure.
The Prais–Winsten transformation retains the first observation and transforms the later ones. The first transformed row is multiplied by \(\sqrt{1-\rho^2}=\sqrt{0.75}\); for \(t\ge2\), subtract \(\rho\) times the preceding value from both the outcome and each regressor. The transformed outcome is \((0.8660,2.5,1.5,6.5)\). The transformed design rows for intercept and slope are \((0.8660,0),(0.5,1),(0.5,1.5),(0.5,2.5)\).
For these transformed data, the residual vector at \((\beta_0,\beta_1)=(1,2)\) is \((0,0,-2,1)\). Its inner product with the intercept column is \(0.5(0+0-2+1)=-0.5\), and with the slope column is \(1(0)+1.5(-2)+2.5(1)=-0.5\). Thus it is not orthogonal to both columns, so \((1,2)\) is not the transformed least-squares estimate. Solving the transformed normal equations gives \((\hat{\beta}_0,\hat{\beta}_1)=(0.5625,2.0625)\), with residual sum of squares \(77/16=4.8125\).
The transformed design cross-product is
There are \(n-k=4-2=2\) residual degrees of freedom. The unbiased estimate of the transformed error variance is \(\hat{\sigma}^2=(77/16)/2=77/32=2.40625\), and the slope variance estimate is \((77/32)(1.5/8)=231/512\approx0.4512\). Therefore \(SE(\hat{\beta}_1)=\sqrt{231/512}=0.6717\). Conditional on the known \(\rho\), Gaussian errors give the usual \(t\) interval with \(df=2\):
The estimated slope is 2.0625 outcome units per predictor unit, but the interval is wide and includes zero. It indicates substantial uncertainty in this small sample, not evidence that the slope must be zero. If \(\rho\) were estimated rather than known, replacing it by an estimate would produce FGLS; the displayed \(t\) interval would not in general be exact because it ignores covariance-parameter estimation.
Worked Example: Random-Intercept Covariance
Worked Example: Random-Intercept Covariance
Suppose two hypothetical groups each contribute two observations. The model is \(y_{ij}=\beta_0+\beta_1x_{ij}+a_j+e_{ij}\), where \(a_j\) is a group-specific random intercept with variance 1 and \(e_{ij}\) are independent errors with variance 1. Assume the group intercepts and observation errors are mutually independent, Gaussian, and mean zero, and that \(E[a_j+e_{ij}\mid X]=0\). The observations are \((x,y)=(0,1),(1,3)\) in group A and \((0,2),(1,5)\) in group B.
Within each group, the shared intercept produces covariance 1 between the two errors, while each marginal variance is \(1+1=2\). Thus each block of \(\Omega\) is
The two groups are independent, so the full covariance is block diagonal. For each group, the contribution to \(X'\Omega^{-1}X\) is \(\frac{1}{3}\begin{pmatrix}2&1\\1&2\end{pmatrix}\), giving total \(\begin{pmatrix}4/3&2/3\\2/3&4/3\end{pmatrix}\). The groups contribute \((4/3,5/3)'\) and \((7/3,8/3)'\) to \(X'\Omega^{-1}y\), respectively, so the total right-hand side is \((11/3,13/3)'\). Solving gives
With the stated covariance treated as known, \((X'\Omega^{-1}X)^{-1}=\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix}\). The slope standard error is therefore 1 outcome unit per predictor unit. Under the Gaussian model, a 95% normal interval is \(2.5\pm1.96(1)\), or \((0.540,4.460)\). It describes uncertainty in the population slope under the random-intercept covariance model. With variance components estimated from data, their estimation would also need to be considered in the inferential procedure.
Choosing and Applying GLS in Practice
A covariance structure should follow the sampling or data-generating arrangement, not the appearance of a convenient answer. Unequal observation precision suggests a diagonal covariance model; time-ordered dependence may motivate an AR(1) model; repeated observations sharing a group-level shock may motivate a random intercept. These are not interchangeable descriptions. A useful workflow is to identify the observational ordering or grouping, state the covariance assumptions, and then construct the corresponding transformation or inverse covariance.
Write \(y=X\beta+u\), specify \(E[u\mid X]=0\), and state \(\operatorname{Var}(u\mid X)=\Omega\), including the source of any covariance parameters.
If \(\Omega\) is specified, use GLS. If it depends on parameters estimated from the sample, substitute those estimates to obtain FGLS and acknowledge plug-in uncertainty.
Use \(\Omega^{-1}\) in the estimator, or a factorization that whitens the errors. Apply the same transformation to the outcome and every regressor, including the intercept.
Use the coefficient covariance implied by the fitted model, and use finite-sample reference distributions only when their conditions hold. Report estimates, uncertainty, and coefficient units.
Common Pitfalls and Why the Distinction Matters
A frequent mistake is to treat a variance as a weight rather than using its inverse. If one observation has error variance four times as large as another, its GLS weight is one quarter as large, not four times as large. Another mistake is to transform only \(y\) while leaving the regressors unchanged; this changes the regression model and does not implement GLS. In AR(1) transformations, the first observation also requires the Prais–Winsten factor if it is to be retained.
A further pitfall is to interpret a fitted GLS coefficient as causal merely because the covariance model is sophisticated. GLS addresses the covariance of the regression errors; it does not by itself establish that regressors are exogenous or remove confounding. Nor does a small standard error certify that the covariance model is correct. Diagnostics and sensitivity to plausible covariance specifications matter, particularly for FGLS where estimated covariance parameters enter the calculation.
| Setting | Covariance feature | GLS implication |
|---|---|---|
| Unequal variances | Diagonal entries differ | Weight observations by inverse variance |
| AR(1) errors | Covariance decays with time separation | Use a serial-correlation transformation or full covariance matrix |
| Random intercept | Observations in a group share covariance | Use block covariance matrices or the equivalent random-effects GLS form |
Check Your Understanding
Use the covariance assumptions, not just the regression formula, to answer each question.
- If two observations have error variances 2 and 8 and are independent, what is the ratio of their inverse-variance weights?
- In a Prais–Winsten transformation with \(\rho=0.4\), what multiplier is applied to the first outcome and first regressor row?
- Why does a random intercept create positive covariance between two observations from the same group?
- What changes when a covariance parameter is estimated and substituted into the GLS estimator?
- Does using GLS alone justify a causal interpretation of the slope? Explain briefly.