Ordinary Least Squares

\[\mathbf{Y = X\beta + \varepsilon}\]

The model is $Y = X\beta + \varepsilon$ with $X \in \mathbb{R}^{n \times p}$, and OLS minimizes $|Y - X\beta|^2$. When $X^\top X$ is invertible the solution is

\[\hat\beta = (X^\top X)^{-1}X^\top Y\]

Proof: Expanding $|Y - X\beta|^2 = Y^\top Y - 2\beta^\top X^\top Y + \beta^\top X^\top X\beta$ and setting the gradient to zero gives $-2X^\top Y + 2X^\top X\beta = 0$, the normal equations $X^\top X\beta = X^\top Y$, so $\hat\beta = (X^\top X)^{-1}X^\top Y$ when $X^\top X$ is invertible.

Moore–Penrose Pseudoinverse

When $X^\top X$ is singular the solution is written through the Moore–Penrose pseudoinverse $X^+$, using the SVD $X = UDV^\top$ ($U, V$ orthogonal, $D$ diagonal with singular values $\sigma_i \geq 0$):

\[\hat\beta = X^+ Y \qquad X^+ = VD^+U^\top \qquad D^+ = \operatorname{diag}\Big(\tfrac{1}{\sigma_i} \text{ for } \sigma_i > 0,\ 0 \text{ otherwise}\Big)\]

This $\hat\beta$ is always the minimum-norm least-squares solution: among all $\beta$ minimizing $|Y - X\beta|^2$ it has the smallest $|\beta|$. Its concrete form depends on the shape of $X$.

Tall, full column rank ($n \geq p$): $X^\top X$ is invertible and $X^+ = (X^\top X)^{-1}X^\top$
Wide, full row rank ($n \leq p$): $XX^\top$ is invertible and $X^+ = X^\top(XX^\top)^{-1}$. There are infinitely many $\beta$ that solve $X\beta = Y$, so this form picks the smallest-norm one.
Rank deficient: Neither inverse exists, only the SVD form $X^+ = VD^+U^\top$ applies, and it still returns the minimum-norm least-squares solution.

Assumptions

Linearity: The model is correctly specified as $Y = X\beta + \varepsilon$.
Strict exogeneity, $\mathbb{E}[\varepsilon \mid X] = 0$: Implies $\mathbb{E}[\varepsilon] = 0$ and $\operatorname{Cov}(X, \varepsilon) = \mathbb{E}[X^\top\varepsilon] = 0$. If it holds the regressors are exogenous, if violated they are endogenous, OLS is biased, and instrumental variables are needed.
No perfect multicollinearity, $\Pr[\operatorname{rank}(X) = p] = 1$: The regressors are linearly independent. If violated $\beta$ is not identified as $(X^\top X)^{-1}$ cannot be computed.
Spherical errors, $\operatorname{Var}[\varepsilon \mid X] = \sigma^2 I_n$: Combines homoscedasticity, $\mathbb{E}[\varepsilon_i^2 \mid X] = \sigma^2$, and no autocorrelation, $\mathbb{E}[\varepsilon_i\varepsilon_j \mid X] = 0$ for $i \neq j$. If violated OLS is unbiased but inefficient, use GLS or robust estimation.
Normality, $\varepsilon \mid X \sim \mathcal{N}(0, \sigma^2 I_n)$: Needed only for exact finite-sample inference, not for unbiasedness.

The Estimator and its Distribution

Substituting $Y = X\beta + \varepsilon$,

\[\hat\beta = (X^\top X)^{-1}X^\top Y = \beta + (X^\top X)^{-1}X^\top\varepsilon\]

so, conditional on $X$,

\[\begin{align*} \mathbb{E}[\hat\beta \mid X] &= \beta + (X^\top X)^{-1}X^\top\,\mathbb{E}[\varepsilon \mid X] = \beta \\ \operatorname{Var}(\hat\beta \mid X) &= (X^\top X)^{-1}X^\top\,\operatorname{Var}(\varepsilon)\,X(X^\top X)^{-1} = \sigma^2(X^\top X)^{-1} \end{align*}\]

Under normality, $\varepsilon \sim \mathcal{N}(0, \sigma^2 I_n)$ so $\hat\beta \sim \mathcal{N}\big(\beta,\ \sigma^2(X^\top X)^{-1}\big)$.

Forecast variance: For a query point $x_0 \in \mathbb{R}^p$ the fitted value is $\hat y_0 = x_0^\top\hat\beta$, with

\[\operatorname{Var}(\hat y_0 \mid X) = x_0^\top\operatorname{Var}(\hat\beta)\,x_0 = \sigma^2\,x_0^\top(X^\top X)^{-1}x_0\]

A prediction of the actual $y_0 = x_0^\top\beta + \varepsilon_0$ adds the new noise,

\[\operatorname{Var}(y_0 - \hat y_0) = \sigma^2\big(1 + x_0^\top(X^\top X)^{-1}x_0\big)\]

Both the estimator and the forecast variance carry $(X^\top X)^{-1}$, so ill-conditioning of $X^\top X$ (multicollinearity) inflates both.

Hat Matrix and Leverage

The fitted values are a linear projection of $Y$, $\hat Y = X\hat\beta = HY$, through the hat matrix

\[H = X(X^\top X)^{-1}X^\top\]

which is symmetric and idempotent ($H^2 = H$), satisfies $HX = X$, and has $\operatorname{tr}(H) = p$. The residuals are $e = (I - H)Y = (I - H)\varepsilon$, so $\operatorname{Var}(e) = \sigma^2(I - H)$ and

\[\operatorname{Var}(e_i) = \sigma^2(1 - h_{ii})\]

where the leverage $h_{ii} \in [0, 1]$ measures how outlying row $i$ is in $X$-space, with $\sum_i h_{ii} = p$. Rescaling each residual by its own standard deviation gives the internally studentized residual $r_i = e_i / \big(\hat\sigma\sqrt{1 - h_{ii}}\big)$. The same $1 - h_{ii}$ governs the leave-one-out identity: the residual from refitting without observation $i$ is the full residual inflated by the leverage,

\[e_{(i)} = \frac{e_i}{1 - h_{ii}}\]

so leave-one-out cross-validation costs only a single full fit,

\[\mathrm{CV} = \frac{1}{n}\sum_{i=1}^n\Big(\frac{e_i}{1 - h_{ii}}\Big)^2\]

Normalizing the deleted residual by a variance estimated without $i$ gives the externally studentized (studentized deleted) residual $t_i = e_i / \big(\hat\sigma_{(i)}\sqrt{1 - h_{ii}}\big)$, used for outlier detection: it is the same leave-one-out correction as the CV residual, applied for testing rather than prediction.

Inference

Write $\mathrm{RSS} = e^\top e = \sum_i(y_i - \hat y_i)^2$ and estimate $\sigma^2$ by

\[\hat\sigma^2 = \frac{\mathrm{RSS}}{n - p} \qquad \mathbb{E}[\hat\sigma^2] = \sigma^2 \qquad \mathrm{RSE} = \sqrt{\hat\sigma^2} = \sqrt{\tfrac{\mathrm{RSS}}{n - p}}\]

where $n - p$ makes it unbiased.

Testing $\beta_j = 0$: The $t$-statistic and confidence interval are

\[t = \frac{\hat\beta_j}{\mathrm{SE}(\hat\beta_j)} \sim t_{n-p} \qquad \mathrm{SE}(\hat\beta_j) = \hat\sigma\sqrt{[(X^\top X)^{-1}]_{jj}} \qquad \hat\beta_j \pm t_{n-p,\,1-\alpha/2}\,\mathrm{SE}(\hat\beta_j)\]

with $t_{n-p} \to \mathcal{N}(0, 1)$ as $n \to \infty$.

Comparing nested models: A restricted model with $p_0$ parameters against a full model with $p_1 > p_0$:

\[F = \frac{(\mathrm{RSS}_0 - \mathrm{RSS}_1)/(p_1 - p_0)}{\mathrm{RSS}_1/(n - p_1)} \sim F_{\,p_1 - p_0,\ n - p_1}\]

Metrics:

\[R^2 = 1 - \frac{\mathrm{RSS}}{\mathrm{TSS}}\ (= \rho_{xy}^2 \text{ in simple regression}) \qquad R^2_{\text{adj}} = 1 - (1 - R^2)\frac{n - 1}{n - p - 1}\] \[\mathrm{AIC} = 2k - 2\ln L \qquad \mathrm{BIC} = k\ln n - 2\ln L\]

The two differ only in the complexity penalty, $2k$ versus $k\ln n$. BIC penalizes harder once $n > e^2 \approx 7.4$, and its penalty grows with the sample size, making it consistent: it recovers the true model as $n \to \infty$. AIC’s fixed penalty is not consistent and tends to over-select, but is asymptotically efficient, minimizing prediction error. So BIC targets the true sparse model, AIC targets prediction.

Gauss–Markov Theorem

Theorem (Gauss–Markov): Under linearity $Y = X\beta + \varepsilon$, $\mathbb{E}[\varepsilon] = 0$, $\operatorname{Var}(\varepsilon) = \sigma^2 I$, and full column rank $X$, the OLS estimator is the best linear unbiased estimator (BLUE): it has the smallest variance among all linear unbiased estimators.

Proof: Let $\tilde\beta = CY$ be any linear estimator, with $C = (X^\top X)^{-1}X^\top + D$. Then

\[\mathbb{E}[\tilde\beta] = \big((X^\top X)^{-1}X^\top + D\big)X\beta = (I + DX)\beta\]

so $\tilde\beta$ is unbiased iff $DX = 0$. Its variance is

\[\begin{align*} \operatorname{Var}(\tilde\beta) = \sigma^2 CC^\top &= \sigma^2\big((X^\top X)^{-1} + (X^\top X)^{-1}X^\top D^\top + DX(X^\top X)^{-1} + DD^\top\big) \\ &= \sigma^2(X^\top X)^{-1} + \sigma^2 DD^\top \end{align*}\]

using $DX = 0$. Since $DD^\top \succeq 0$,

\[\operatorname{Var}(\tilde\beta) \succeq \sigma^2(X^\top X)^{-1} = \operatorname{Var}(\hat\beta)\]

so $\hat\beta$ has the smallest variance and is BLUE.

Frisch–Waugh–Lovell Theorem

For $Y = X_1\beta_1 + X_2\beta_2 + \varepsilon$, the coefficient $\hat\beta_2$ is unchanged if one first partials out $X_1$ from both $Y$ and $X_2$ and regresses the residuals. Let $M_1$ be the residual-maker (annihilator), projecting onto the orthogonal complement of the column space of $X_1$,

\[M_1 = I - X_1(X_1^\top X_1)^{-1}X_1^\top\]

Then $\hat\beta_2$ from the full regression equals $\tilde\beta_2$ from $M_1Y = \beta_2\,M_1X_2 + M_1\varepsilon$. The procedure:

Residualize $Y$ on $X_1$: $\tilde Y = M_1Y = Y - X_1(X_1^\top X_1)^{-1}X_1^\top Y$
Residualize $X_2$ on $X_1$: $\tilde X_2 = M_1X_2 = X_2 - X_1(X_1^\top X_1)^{-1}X_1^\top X_2$
Regress $\tilde Y$ on $\tilde X_2$: $\tilde\beta_2 = (\tilde X_2^\top\tilde X_2)^{-1}\tilde X_2^\top\tilde Y = \hat\beta_2$

Violated Assumptions

Model Misspecification / Functional Form

Problem: Wrong functional form, omitted variables, or unmodelled nonlinearity, which breaks $\mathbb{E}[\varepsilon \mid X] = 0$ when the true relationship is not the fitted one.

Why it matters: $\hat\beta$ is biased and inconsistent. With $Y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon$ but only $x_1$ regressed, the omitted-variable bias is

\[\mathbb{E}[\hat\beta_1] = \beta_1 + \beta_2\,\frac{\operatorname{Cov}(x_1, x_2)}{\operatorname{Var}(x_1)}\]

vanishing only if $x_1 \perp x_2$ or $\beta_2 = 0$.

Diagnostics: Residuals versus fitted values or versus each $x_j$ should not show any patterns, so any systematic shapes (U-shapes, trends) signal misspecification. Partial residual plots, $e + \hat\beta_j x_j$ against $x_j$, isolate nonlinearity in $x_j$ controlling for the rest.

Remedies: Nonlinear transforms of predictors ($\log x$, $\sqrt{x}$, $x^2$, $1/x$), interaction terms $x_1 x_2$, and the Box–Cox transform of the response $y^{(\lambda)} = (y^\lambda - 1)/\lambda$ (or $\log y$ at $\lambda = 0$), with $\lambda$ chosen by maximum likelihood.

Multicollinearity

Problem: High correlation among regressors makes $X^\top X$ ill-conditioned.

Why it matters: $\operatorname{Var}(\hat\beta) = \sigma^2(X^\top X)^{-1}$ inflates, deflating $t$-statistics and raising the Type II rate, and the individual effects become impossible to disentangle. The forecast variance $\sigma^2 x_0^\top(X^\top X)^{-1}x_0$ inflates for the same reason.

Diagnostics: With $R_j^2$ from regressing $x_j$ on the other predictors, the variance inflation factor multiplies the variance relative to the orthogonal case:

\[\mathrm{VIF}_j = \frac{1}{1 - R_j^2} \qquad \operatorname{Var}(\hat\beta_j) = \operatorname{Var}(\hat\beta_j)^{\text{orth}} \times \mathrm{VIF}_j\]

The condition number $\kappa = \sqrt{\lambda_{\max}/\lambda_{\min}}$ of $X^\top X$ above $30$ also signals a problem.

Remedies: Shrinkage (Ridge, Lasso, Elastic-Net), dimensionality reduction (principal component regression), FWL, drop or combine redundant predictors.

Heteroskedasticity

Problem: Non-constant error variance, $\operatorname{Var}(\varepsilon_i) = \sigma_i^2$. OLS stays unbiased but is no longer BLUE, and $\sigma^2(X^\top X)^{-1}$ is the wrong covariance, invalidating inference.

Diagnostics: Plot residuals and squared residuals against predictors or fitted values (a flat band around $\sigma^2$ is expected). Formal tests: Breusch–Pagan (variance depends on predictors), White (general), Goldfeld–Quandt (variance differs across two subsets).

Remedies: Robust Eicker–Huber–White standard errors, replacing $\sigma^2(X^\top X)^{-1}$ by the sandwich

\[\widehat V_{\mathrm{EHW}} = (X^\top X)^{-1}\big(X^\top\hat\Omega X\big)(X^\top X)^{-1} \qquad \hat\Omega = \operatorname{diag}(\hat\varepsilon_1^2, \ldots, \hat\varepsilon_n^2)\]

after which $\hat\beta \overset{a}{\sim} \mathcal{N}(\beta, \widehat V_{\mathrm{EHW}})$. Weighted least squares with $w_i = 1/\sigma_i^2$ if the variances are known.

Autocorrelated Errors

Problem: Non-zero covariance between errors, $\operatorname{Cov}(\varepsilon_i, \varepsilon_j) \neq 0$, very common in time series.

Why it matters: Consecutive observations are not independent, cutting the effective sample size. For AR(1) the variance is inflated by roughly $(1+\rho)/(1-\rho)$. OLS is unbiased but inefficient, and conventional standard errors are inconsistent: for $\rho > 0$ they are underestimated (narrow intervals, Type I errors), for $\rho < 0$ overestimated (Type II).

Diagnostics: Residuals against time and ACF / PACF plots. Durbin–Watson

\[d = \frac{\sum_{t \geq 2}(e_t - e_{t-1})^2}{\sum_t e_t^2} \approx 2(1 - \hat\rho)\]

with $d \approx 2$ no autocorrelation, $d \ll 2$ positive, $d \gg 2$ negative (AR(1) only, inconclusive near $2$). Box–Pierce $Q_{BP} = n\sum_{k=1}^h\hat\rho_k^2 \sim \chi^2_h$ and the small-sample Ljung–Box $Q_{LB} = n(n+2)\sum_{k=1}^h\hat\rho_k^2/(n-k) \sim \chi^2_h$.

Remedies: Generalized least squares: for AR(1), transform $y_t^* = y_t - \rho y_{t-1}$, $X_t^* = X_t - \rho X_{t-1}$, with Cochrane–Orcutt estimating $\rho$ from residuals iteratively. Including a lagged dependent variable can absorb the autocorrelation into model dynamics. Or Newey–West HAC standard errors, below.

Newey–West HAC: The heteroskedasticity- and autocorrelation-consistent covariance keeps the sandwich form,

\[\widehat V_{\mathrm{NW}} = (X^\top X)^{-1}\big(X^\top\hat\Omega X\big)(X^\top X)^{-1}\]

but the meat now carries residual autocovariances out to lag $L$, down-weighted by Bartlett weights $w_k = 1 - \tfrac{k}{L+1}$. The residual covariance $\hat\Omega$ is banded, with $e_i^2$ on the diagonal and $w_k e_i e_{i\pm k}$ on the $k$-th band, zero beyond $L$:

\[\hat\Omega = \begin{pmatrix} e_1^2 & w_1 e_1 e_2 & w_2 e_1 e_3 & \cdots & 0\\ w_1 e_2 e_1 & e_2^2 & w_1 e_2 e_3 & \ddots & \vdots\\ w_2 e_3 e_1 & w_1 e_3 e_2 & e_3^2 & \ddots & w_2 e_{n-2}e_n\\ \vdots & \ddots & \ddots & \ddots & w_1 e_{n-1}e_n\\ 0 & \cdots & w_2 e_n e_{n-2} & w_1 e_n e_{n-1} & e_n^2 \end{pmatrix}\]

The weights decay linearly, the nearest lag counting most: for $L = 5$ the lags $1, \dots, 5$ get $\tfrac{5}{6}, \tfrac{4}{6}, \tfrac{3}{6}, \tfrac{2}{6}, \tfrac{1}{6}$. This linear taper is what guarantees a positive-semidefinite estimate.

Heavy-tailed / Non-normal Errors

Problem: Errors deviate from normality, especially heavy tails (high kurtosis). Normality is not needed for unbiasedness, but it matters for inference and efficiency.

Why it matters: $\hat\sigma^2$ is sensitive to extremes and is typically overestimated under heavy tails, inflating standard errors and cutting power. Exact $t$- and $F$-tests need normality, otherwise they rely on the CLT, which fails for infinite-variance distributions, so inference is invalid even asymptotically.

Diagnostics: Q–Q plot (tails curve away at the extremes), excess kurtosis, Shapiro–Wilk $W = \frac{(\sum_i a_i x_{(i)})^2}{\sum_i(x_i - \bar x)^2}$ (most powerful for $n < 2000$), and the goodness-of-fit tests Kolmogorov–Smirnov and Anderson–Darling (the latter weighting the tails more).

Remedies: Robust standard errors (HC, HAC) and bootstrap standard errors (residual or pairs), which stay valid asymptotically without a distributional assumption, report studentized residuals $r_i = e_i/(\hat\sigma\sqrt{1 - h_{ii}})$.

References

The Elements of Statistical Learning: Data Mining, Inference, and Prediction. (2009)
Trevor Hastie, Robert Tibshirani, Jerome Friedman
Book

The Truth about Linear Regression (2015)
Cosma Rohilla Shalizi
Book

Linear Model and Extensions (2024)
Peng Ding
Book