OLS Algebra

The Gauss Markov Model

Definition

A statistical model for regression data is the Gauss Markov Model if each of its distributions satisfies the conditions

1. Linearity: a statistical model $\mathcal{F}$ over data $\mathcal{D}$ satisfies linearity if for each element of $\mathcal{F}$, the data can be decomposed in $$\begin{array}{rl}{y}_i& ={\beta }_1{x}_{i1}+\cdots +{\beta }_k{x}_{ik}+{\epsilon }_i=x_i^{\prime }\beta +{\epsilon }_i\\ \underset{n\times 1}{\phantom{{\beta }_{\beta }}y}& =\underset{n\times k}{\phantom{\beta }X}\cdot \underset{k\times 1}{\beta }+\underset{n\times 1}{\phantom{\beta }\epsilon }\end{array}$$

2. Strict Exogeneity: $\mathbb{E}[{\epsilon }_i|x_1,\dots ,x_n]=0,\mathrm{\forall }i$.

3. No Multicollinerity: ${\mathbb{E}}_n[x_i{x}_i^{\prime }]$ is strictly positive definite almost surely. Equivalent to require $rank(X)=k$ with probability $p\to 1$. Intuition: no regressor is a linear combination of other regressors.

4. Spherical Error Variance: -$\mathbb{E}[{\epsilon }_i^2|x]={\sigma }^2>0,\mathrm{\forall }i$ -$\mathbb{E}[{\epsilon }_i{\epsilon }_j|x]=0,\mathrm{\forall }$ $1\le i<j\le n$

The Extended Gauss Markov Model also satisfies assumption

1. Normal error term: $\epsilon |X\sim N(0,{\sigma }^2{I}_n)$ and $\epsilon \perp X$.

Implications

• Note that by (2) and (4) you get homoskedasticity:

$$Var({\epsilon }_i|x)=\mathbb{E}[{\epsilon }_i^2|x]-\mathbb{E}[{\epsilon }_i|x{]}^2={\sigma }^{2}I\mathrm{\forall }i$$

• Strict exogeneity is not restrictive since it is sufficient to include a constant in the regression to enforce it $$y_i=\alpha +x_i^{\prime }\beta +({\epsilon }_i-\alpha )\Rightarrow \mathbb{E}[{\epsilon }_i]={\mathbb{E}}_x[\mathbb{E}[{\epsilon }_i|x]]=0$$
• This implies $\mathbb{E}[x_{jk}{\epsilon }_i]=0$ by the LIE.
• These two conditions together imply $Cov(x_{jk}{\epsilon }_i)=0$.

Projection

A map $\mathrm{\Pi }:V\to V$ is a projection if $\mathrm{\Pi }\circ \mathrm{\Pi }=\mathrm{\Pi }$.

The Gauss Markov Model assumes that the conditional expectation function (CEF) $f(X)=\mathbb{E}[Y|X]$ and the linear projection $g(X)=X\beta$ coincide.

Code - DGP

This code draws 100 observations from the model $y=2x_1-x_2+\epsilon$ where $x_1,x_2\sim U[0,1]$ and $\epsilon \sim N(0,1)$.

# Set seed
Random.seed!(123);

# Set the number of observations
n = 100;

# Set the dimension of X
k = 2;

# Draw a sample of explanatory variables
X = rand(Uniform(0,1), n, k);

# Draw the error term
σ = 1;
ε = rand(Normal(0,1), n, 1) * sqrt(σ);

# Set the parameters
β = [2; -1];

# Calculate the dependent variable
y = X*β + ε;

The OLS estimator

Definition

The sum of squared residuals (SSR) is given by $$Q_n(\beta )\equiv \frac{1}{n}\sum _{i=1}^n{(y_i-x_i^{\prime }\beta )}^2=\frac{1}{n}(y-X\beta {)}^{\prime }(y-X\beta )$$

Consider a dataset $\mathcal{D}$ and define $Q_n(\beta )={\mathbb{E}}_n[(y_i-x_i^{\prime }\beta {)}^2]$. Then the ordinary least squares (OLS) estimator ${\hat{\beta }}_{OLS}$ is the value of $\beta$ that minimizes $Q_n(\beta )$.

When we can write $D=(y,X)$ in matrix form, then $${\hat{\beta }}_{OLS}=\arg \underset{\beta }{min}\frac{1}{n}(y-X\beta {)}^{\prime }(y-X\beta )$$

Derivation

Theorem

Under the assumption that $X$ has full rank, the OLS estimator is unique and it is determined by the normal equations. More explicitly, $\hat{\beta }$ is the OLS estimate precisely when $X^{\prime }X\hat{\beta }=X^{\prime }y$.

Proof

Taking the FOC: $$\frac{\partial Q_n(\beta )}{\partial \beta }=-\frac{2}{n}{X}^{\prime }y+\frac{2}{n}{X}^{\prime }X\beta =0\iff X^{\prime }X\beta =X^{\prime }y$$ Since $(X^{\prime }X{)}^{-1}$ exists by assumption,

Finally, $\frac{{\partial }^2{Q}_n(\beta )}{\partial \beta \partial {\beta }^{\prime }}=X^{\prime }X/n$ is positive definite since $X^{\prime }X$ is positive semi-definite and $(X^{\prime }X{)}^{-1}$ exists because $X$ is full rank. Therefore, $Q_n(\beta )$ minimized at ${\hat{\beta }}_n$. $$\begin{array}{}◼& \end{array}$$

The $k$ equations $X^{\prime }X\hat{\beta }=X^{\prime }y$ are called normal equations.

Futher Objects

• Fitted coefficient: ${\hat{\beta }}_{OLS}=(X^{\prime }X{)}^{-1}{X}^{\prime }y={\mathbb{E}}_n[x_i{x}_i^{\prime }]{\mathbb{E}}_n[x_i{y}_i]$
• Fitted residual: ${\hat{\epsilon }}_i=y_i-x_i^{\prime }\hat{\beta }$
• Fitted value: ${\hat{y}}_i=x_i^{\prime }\hat{\beta }$
• Predicted coefficient: ${\hat{\beta }}_{-i}={\mathbb{E}}_n[x_{-i}{x}_{-i}^{\prime }]{\mathbb{E}}_n[x_{-i}{y}_{-i}]$
• Prediction error: ${\hat{\epsilon }}_{-i}=y_i-x_i^{\prime }{\hat{\beta }}_{-i}$
• Predicted value: ${\hat{y}}_i=x_i^{\prime }{\hat{\beta }}_{-i}$

Notes on Orthogonality Conditions

• The normal equations are equivalent to the moment condition ${\mathbb{E}}_n[x_i{\epsilon }_i]=0$.
• The algebraic result ${\mathbb{E}}_n[x_i{\hat{\epsilon }}_i]=0$ is called ortogonality property of the OLS residual ${\hat{\epsilon }}_i$.
• If we have included a constant in the regression, ${\mathbb{E}}_n[{\hat{\epsilon }}_i]=0$.
• $\mathbb{E}[{\mathbb{E}}_n[x_i{\epsilon }_i]]=0$ by strict exogeneity (assumed in GM), but ${\mathbb{E}}_n[x_i{\epsilon }_i]\ne \mathbb{E}[x_i{\epsilon }_i]=0$. This is why ${\hat{\beta }}_{OLS}$ is just an estimate of ${\beta }_0$.
• Calculating OLS is like replacing the $j$ equations $\mathbb{E}[x_{ij}{\epsilon }_i]=0$ $\mathrm{\forall }j$ with ${\mathbb{E}}_n[x_{ij}{\epsilon }_i]=0$ $\mathrm{\forall }j$ and forcing them to hold (remindful of GMM).

The Projection Matrix

The projection matrix is given by $P=X(X^{\prime }X{)}^{-1}{X}^{\prime }$. It has the following properties: - $PX=X$ - $P\hat{\epsilon }=0$ ($P$, $\epsilon$ orthogonal) - $Py=X(X^{\prime }X{)}^{-1}{X}^{\prime }y=X\hat{\beta }=\hat{y}$ - Symmetric: $P=P^{\prime }$, Idempotent: $PP=P$ - $tr(P)=tr(X(X^{\prime }X{)}^{-1}{X}^{\prime })=tr(X^{\prime }X(X^{\prime }X{)}^{-1})=tr(I_k)=k$ - Its diagonal elements are $h_{ii}=x_i(X^{\prime }X{)}^{-1}{x}_i^{\prime }$ and are called leverage.

$h_{ii}\in [0,1]$ is a normalized length of the observed regressor vector $x_i$. In the OLS regression framework it captures the relative influence of observation $i$ on the estimated coefficient. Note that $\sum _n{h}_{ii}=k$.

The Annihilator Matrix

The annihilator matrix is given by $M=I_n-P$. It has the following properties: - $MX=0$ ($M$, $X$ orthogonal) - $M\hat{\epsilon }=\hat{\epsilon }$ - $My=\hat{\epsilon }$ - Symmetric: $M=M^{\prime }$, idempotent: $MM=M$ - $tr(M)=n-k$ - Its diagonal elements are $1-h_{ii}\in [0,1]$

Then we can equivalently write $\hat{y}$ (defined by stacking ${\hat{y}}_i$ into a vector) as $\hat{y}=Py$.

Estimating Beta

# Estimate beta
β_hat = inv(X'*X)*(X'*y)

## 2×1 Array{Float64,2}:
##   1.8821600407711814
##  -0.9429354944506099

# Equivalent but faster formulation
β_hat = (X'*X)\(X'*y)

## 2×1 Array{Float64,2}:
##   1.8821600407711816
##  -0.9429354944506098

# Even faster (but less intuitive) formulation
β_hat = X\y

## 2×1 Array{Float64,2}:
##   1.8821600407711807
##  -0.9429354944506088

Equivalent Formulation?

Generally it’s not true that $${\hat{\beta }}_{OLS}=\frac{Var(X)}{Cov(X,y)}$$

# Wrong formulation
β_wrong = inv(cov(X)) * cov(X, y)

## 2×1 Array{Float64,2}:
##   1.8490257777704475
##  -0.9709213554007003

Equivalent Formulation (correct)

But it’s true if you include a constant, $\alpha$ $$y=\alpha +X\beta +\epsilon$$

# Correct, with constant
α = 3;
y1 = α .+ X*β + ε;
β_hat1 = [ones(n,1) X] \ y1

## 3×1 Array{Float64,2}:
##   3.0362313477745615
##   1.8490257777704477
##  -0.9709213554007007

β_correct1 = inv(cov(X)) * cov(X, y1)

## 2×1 Array{Float64,2}:
##   1.8490257777704477
##  -0.9709213554007006

Some More Objects

# Predicted y
y_hat = X*β_hat;

# Residuals
ε_hat = y - X*β_hat;

# Projection matrix
P = X * inv(X'*X) * X';

# Annihilator matrix
M = I - P;

# Leverage
h = diag(P);

OLS Residuals

Homoskedasticity

The error is homoskedastic if $\mathbb{E}[{\epsilon }^2|x]={\sigma }^2$ does not depend on $x$. $$Var(\epsilon )=I{\sigma }^2=\left[\begin{array}{ccc}{\sigma }^2& \dots & 0\\ \\ ⋮& \ddots & ⋮\\ 0& \dots & {\sigma }^2\end{array}\right]$$

The error is heteroskedastic if $\mathbb{E}[{\epsilon }^2|x]={\sigma }^2(x)$ does depend on $x$. $$Var(\epsilon )=I{\sigma }_i^2=\left[\begin{array}{ccc}{\sigma }_1^2& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & {\sigma }_n^2\end{array}\right]$$

Residual Variance

The OLS residual variance can be an object of interest even in a heteroskedastic regression. Its method of moments estimator is given by $${\hat{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n{\hat{\epsilon }}_i^2$$

Note that ${\hat{\sigma }}^2$ can be rewritten as $${\hat{\sigma }}^2=\frac{1}{n}{\epsilon }^{\prime }{M}^{\prime }M\epsilon =\frac{1}{n}tr({\epsilon }^{\prime }M\epsilon )=\frac{1}{n}tr(M{\epsilon }^{\prime }\epsilon )$$

However, the method of moments estimator is a biesed estimator. In fact $$\mathbb{E}[{\hat{\sigma }}^2|X]=\frac{1}{n}\mathbb{E}[tr(M{\epsilon }^{\prime }\epsilon )|X]=\frac{1}{n}tr(M\mathbb{E}[{\epsilon }^{\prime }\epsilon |X])=\frac{1}{n}\sum _{i=1}^n(1-h_{ii}){\sigma }_i^2$$

Under conditional homoskedasticity, the above expression simplifies to $$\mathbb{E}[{\hat{\sigma }}^2|X]=\frac{1}{n}tr(M){\sigma }^2=\frac{n-k}{n}{\sigma }^2$$

Sample Variance

The OLS residual sample variance is denoted by $s^2$ and is given by $$s^2=\frac{SSR}{n-k}=\frac{{\hat{\epsilon }}^{\prime }\hat{\epsilon }}{n-k}=\frac{1}{n-k}\sum _{i=1}^n{\hat{\epsilon }}_i^2$$ Furthermore, the square root of $s^2$, denoted $s$, is called the standard error of the regression (SER) or the standard error of the equation (SEE). Not to be confused with other notions of standard error to be defined later in the course.

The sum of squared residuals can be rewritten as: $SSR={\hat{\epsilon }}^{\prime }\hat{\epsilon }={\epsilon }^{\prime }M\epsilon$.

The OLS residual sample variance is an unbiased estimator of the error variance ${\sigma }^2$.

Another unbiased estimator of ${\sigma }^2$ is given by $${\overline{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n(1-h_{ii}{)}^{-1}{\hat{\epsilon }}_i^2$$

Uncentered R^2

One measure of the variability of the dependent variable $y_i$ is the sum of squares $\sum _{i=1}^n{y}_i^2=y^{\prime }y$. There is a decomposition: $$\begin{array}{rl}{y}^{\prime }y& =(\hat{y}+e{)}^{\prime }(\hat{y}+\hat{\epsilon })\\ & ={\hat{y}}^{\prime }\hat{y}+2{\hat{y}}^{\prime }\hat{\epsilon }+{\hat{\epsilon }}^{\prime }\hat{\epsilon }e\\ & ={\hat{y}}^{\prime }\hat{y}+2b^{\prime }{X}^{\prime }\hat{\epsilon }+{\hat{\epsilon }}^{\prime }\hat{\epsilon }(\text{since}\hat{y}=Xb)\\ & ={\hat{y}}^{\prime }\hat{y}+{\hat{\epsilon }}^{\prime }\hat{\epsilon }(\text{since}{X}^{\prime }\hat{\epsilon }=0)\end{array}$$

The uncentered ${\mathbf{R}}^{\mathbf{2}}$ is defined as: $$R_{uc}^2\equiv 1-\frac{{\hat{\epsilon }}^{\prime }\hat{\epsilon }}{y^{\prime }y}=1-\frac{{\mathbb{E}}_n[{\hat{\epsilon }}_i^2]}{{\mathbb{E}}_n[y_i^2]}=\frac{\mathbb{E}[{\hat{y}}_i^2]}{\mathbb{E}[y_i^2]}$$

Centered R^2

A more natural measure of variability is the sum of centered squares $\sum _{i=1}^n(y_i-\overline{y}{)}^2,$ where $\overline{y}:=\frac{1}{n}\sum _{i=1}^n{y}_i$. If the regressors include a constant, it can be decomposed as $$\sum _{i=1}^n(y_i-\overline{y}{)}^2=\sum _{i=1}^n({\hat{y}}_i-\overline{y}{)}^2+\sum _{i=1}^n{\hat{\epsilon }}_i^2$$

The coefficient of determination, ${\mathbf{R}}^{\mathbf{2}}$, is defined as $$R^2\equiv 1-\frac{\sum _{i=1}^n{\hat{\epsilon }}_i^2}{\sum _{i=1}^n(y_i-\overline{y}{)}^2}=\frac{\sum _{i=1}^n({\hat{y}}_i-\overline{y}{)}^2}{\sum _{i=1}^n(y_i-\overline{y}{)}^2}=\frac{{\mathbb{E}}_n[({\hat{y}}_i-\overline{y}{)}^2]}{{\mathbb{E}}_n[(y_i-\overline{y}{)}^2]}$$

Always use the centered $R^2$ unless you really know what you are doing.

Code - Variance

# Biased variance estimator
σ_hat = ε_hat'*ε_hat / n;

# Unbiased estimator 1
σ_hat_2 = ε_hat'*ε_hat / (n-k);

# Unbiased estimator 2
σ_hat_3 = mean( ε_hat.^2 ./ (1 .- h) );

Code - R^2

# R squared - uncentered
R2_uc = (y_hat'*y_hat)/ (y'*y);

# R squared
y_bar = mean(y);
R2 = ((y_hat .- y_bar)'*(y_hat .- y_bar))/ ((y .- y_bar)'*(y .- y_bar));

Finite Sample Properties of OLS

Conditional Unbiasedness

Theorem

Under the GM assumptions (1)-(3), the OLS estimator is conditionally unbiased, i.e. the distribution of ${\hat{\beta }}_{OLS}$ is centered at ${\beta }_0$: $\mathbb{E}[\hat{\beta }|X]={\beta }_0$.

Proof $$\begin{array}{rl}\mathbb{E}[\hat{\beta }|X]& =\mathbb{E}[(X^{\prime }X{)}^{-1}{X}^{\prime }y|X]=\\ & =(X^{\prime }X{)}^{-1}{X}^{\prime }\mathbb{E}[y|X]=\\ & =(X^{\prime }X{)}^{-1}{X}^{\prime }\mathbb{E}[X\beta +\epsilon |X]=\\ & =(X^{\prime }X{)}^{-1}{X}^{\prime }X\beta +(X^{\prime }X{)}^{-1}{X}^{\prime }\mathbb{E}[\epsilon |X]=\\ & =\beta \end{array}$$ $$\begin{array}{}◼& \end{array}$$

OLS Variance

Theorem

Under the GM assumptions (1)-(3), $Var(\hat{\beta }|X)={\sigma }^2(X^{\prime }X{)}^{-1}$.

Proof: $$\begin{array}{rl}Var(\hat{\beta }|X)& =Var((X^{\prime }X{)}^{-1}{X}^{\prime }y|X)=\\ & =((X^{\prime }X{)}^{-1}{X}^{\prime })Var(y|X)((X^{\prime }X{)}^{-1}{X}^{\prime }{)}^{\prime }=\\ & =((X^{\prime }X{)}^{-1}{X}^{\prime })Var(X\beta +\epsilon |X)((X^{\prime }X{)}^{-1}{X}^{\prime }{)}^{\prime }=\\ & =((X^{\prime }X{)}^{-1}{X}^{\prime })Var(\epsilon |X)((X^{\prime }X{)}^{-1}{X}^{\prime }{)}^{\prime }=\\ & =((X^{\prime }X{)}^{-1}{X}^{\prime }){\sigma }^{2}I((X^{\prime }X{)}^{-1}{X}^{\prime }{)}^{\prime }=\\ & ={\sigma }^2(X^{\prime }X{)}^{-1}\end{array}$$ $$\begin{array}{}◼& \end{array}$$

Higher correlation of the $X$ implies higher variance of the OLS estimator.

Intuition: individual observations carry less information. You are exploring a smaller region of the $X$ space.

BLUE

Theorem

Under the GM assumptions (1)-(3), $Cov(\hat{\beta },\hat{\epsilon })=0$.

Theorem

Under the GM assumptions (1)-(3), ${\hat{\beta }}_{OLS}$ is the best (most efficient) linear, unbiased estimator (BLUE), i.e., for any unbiased linear estimator $b$: $Var(b|X)\ge Var(\hat{\beta }|X)$.

BLUE Proof

Consider four steps:

1. Define three objects: (i) $b=Cy$, (ii) $A=(X^{\prime }X{)}^{-1}{X}^{\prime }$ such that $\hat{\beta }=Ay$, and (iii) $D=C-A$.
2. Decompose $b$ as $$\begin{array}{rl}b& =(D+A)y=\\ & =Dy+Ay=\\ & =D(X\beta +\epsilon )+\hat{\beta }=\\ & =DX\beta +D\epsilon +\hat{\beta }\end{array}$$
3. By assumption, $b$ must be unbiased: $$\begin{array}{rl}\mathbb{E}[b|X]& =\mathbb{E}[D(X\beta +\epsilon )+Ay|X]=\\ & =\mathbb{E}[DX\beta |X]+\mathbb{E}[D\epsilon |X]+\mathbb{E}[\hat{\beta }|X]=\\ & =DX\beta +D\mathbb{E}[\epsilon |X]+\beta \\ & =DX\beta +\beta \end{array}$$ Hence, it must be that $DX=0$

BLUE Proof (2)

1. We know by (2)-(3) that $b=D\epsilon +\hat{\beta }$. We can now calculate its variance. $$\begin{array}{rl}Var(b|X)& =Var(\hat{\beta }+D\epsilon |X)=\\ & =Var(Ay+D\epsilon |X)=\\ & =Var(AX\beta +(D+A)\epsilon |X)=\\ & =Var((D+A)\epsilon |X)=\\ & =(D+A){\sigma }^{2}I(D+A{)}^{\prime }=\\ & ={\sigma }^{2}I(D{D}^{\prime }+A{A}^{\prime }+D{A}^{\prime }+A{D}^{\prime })=\\ & ={\sigma }^{2}I(D{D}^{\prime }+A{A}^{\prime })\ge \\ & \ge {\sigma }^{2}A{A}^{\prime }=\\ & ={\sigma }^2(X^{\prime }X{)}^{-1}=\\ & =Var(\hat{\beta }|X)\end{array}$$ since $D{A}^{\prime }=A{D}^{\prime }=0$, $DX=0$ and $A{A}^{\prime }=(X^{\prime }X{)}^{-1}$. $$\begin{array}{}◼& \end{array}$$

$Var(b|X)\ge Var(\hat{\beta }|X)$ is meant in a positive definite sense.

Code - Variance

# Ideal variance of the OLS estimator
var_β = σ * inv(X'*X)

## 2×2 Array{Float64,2}:
##   0.0609402  -0.0467732
##  -0.0467732   0.0656808

# Standard errors
std_β = sqrt.(diag(var_β))

## 2-element Array{Float64,1}:
##  0.24686077212177054
##  0.25628257446345265