OLS Inference

Asymptotic Theory of the OLS Estimator

OLS Consistency

Theorem: Assume that $(x_i,y_i{)}_{i=1}^n$ i.i.d. , $\mathbb{E}[x_i{x}_i^{\prime }]=Q$ positive definite, $\mathbb{E}[x_i{x}_i^{\prime }]<\mathrm{\infty }$ and $\mathbb{E}[y_i^2]<\mathrm{\infty }$, then ${\hat{\beta }}_{OLS}$ is a consistent estimator of ${\beta }_0$, i.e. $\hat{\beta }={\mathbb{E}}_n[x_i{x}_i^{\prime }]{\mathbb{E}}_n[x_i{y}_i]\stackrel{p}{\to }{\beta }_0$.

Proof:
We consider 4 steps:

1. ${\mathbb{E}}_n[x_i{x}_i^{\prime }]\stackrel{p}{\to }\mathbb{E}[x_i{x}_i^{\prime }]$ by WLLN since $x_i{x}_i^{\prime }$ iid and $\mathbb{E}[x_i{x}_i^{\prime }]<\mathrm{\infty }$.
2. ${\mathbb{E}}_n[x_i{y}_i]\stackrel{p}{\to }\mathbb{E}[x_i{y}_i]$ by WLLN, due to $x_i{y}_i$ iid, Cauchy-Schwarz and finite second moments of $x_i$ and $y_i$ $$\mathbb{E}\left[x_i{y}_i\right]\le \sqrt{\mathbb{E}[x_i^2]\mathbb{E}[y_i^2]}<\mathrm{\infty }$$
3. ${\mathbb{E}}_n[x_i{x}_i^{\prime }{]}^{-1}\stackrel{p}{\to }\mathbb{E}[x_i{x}_i^{\prime }{]}^{-1}$ by CMT.
4. ${\mathbb{E}}_n[x_i{x}_i^{\prime }{]}^{-1}{\mathbb{E}}_n[x_i{y}_i]\stackrel{p}{\to }\mathbb{E}[x_i{x}_i^{\prime }{]}^{-1}\mathbb{E}[x_i{y}_i]=\beta$ by CMT. $$\begin{array}{}◼& \end{array}$$

Variance and Assumptions

Now we are going to investigate the variance of ${\hat{\beta }}_{OLS}$ progressively relaxing the underlying assumptions.

• Gaussian error term.
• Homoskedastic error term.
• Heteroskedastic error term.
• Heteroskedastic and autocorrelated error term.

Gaussian Error Term

Theorem: Under the GM assumption (1)-(5), $\hat{\beta }-\beta |X\sim N(0,{\sigma }^2(X^{\prime }X{)}^{-1})$

Proof:
We follow 2 steps:

1. We can rewrite $\hat{\beta }$ as $$\begin{array}{rl}\hat{\beta }& =(X^{\prime }X{)}^{-1}{X}^{\prime }y=(X^{\prime }X{)}^{-1}{X}^{\prime }(X\beta +\epsilon )\\ & =\beta +(X^{\prime }X{)}^{-1}{X}^{\prime }\epsilon =\\ & =\beta +{\mathbb{E}}_n[x_i{x}_i^{\prime }{]}^{-1}{\mathbb{E}}_n[x_i{\epsilon }_i]\end{array}$$
2. Therefore: $\hat{\beta }-\beta ={\mathbb{E}}_n[x_i{x}_i^{\prime }{]}^{-1}{\mathbb{E}}_n[x_i{\epsilon }_i]$. $$\begin{array}{rl}\hat{\beta }-\beta |X& \sim (X^{\prime }X{)}^{-1}{X}^{\prime }N(0,{\sigma }^2{I}_n)=\\ & =N(0,{\sigma }^2(X^{\prime }X{)}^{-1}{X}^{\prime }X(X^{\prime }X{)}^{-1})=\\ & =N(0,{\sigma }^2(X^{\prime }X{)}^{-1})\end{array}$$ $$\begin{array}{}◼& \end{array}$$

Does it make sense to assume that $\epsilon$ is gaussian? Not much. But does it make sense that $\hat{\beta }$ is gaussian? Yes, because it’s an average.

Homoskedastic Error Term

Theorem: Under the assumptions of the previous theorem, plus $\mathbb{E}[x^4]<\mathrm{\infty }$, the OLS estimate has an asymptotic normal distribution: $\hat{\beta }|X\stackrel{d}{\to }N(\beta ,{\sigma }^2(X^{\prime }X{)}^{-1})$.

Proof: $$\sqrt{n}(\hat{\beta }-\beta )=\underset{\stackrel{p}{\to }{Q}^{-1}}{\underbrace{{\mathbb{E}}_n[x_i{x}_i^{\prime }{]}^{-1}}}\underset{\stackrel{d}{\to }N(0,\mathrm{\Omega })}{\underbrace{\sqrt{n}{\mathbb{E}}_n[x_i{\epsilon }_i]}}\to N(0,\mathrm{\Sigma })$$ where in general $\mathrm{\Omega }=Var(x_i{\epsilon }_i)=\mathbb{E}[(x_i{\epsilon }_i{)}^2]$ and $\mathrm{\Sigma }=Q^{-1}\mathrm{\Omega }{Q}^{-1}$. $$\begin{array}{}◼& \end{array}$$

Given that $Q=\mathbb{E}[x_i{x}_i^{\prime }]$ is unobserved, we estimate it with $\hat{Q}={\mathbb{E}}_n[x_i{x}_i^{\prime }]$. Since we have assumed homoskedastic error term, we have $\mathrm{\Omega }={\sigma }^2(X^{\prime }X{)}^{-1}$. Since we do not observe ${\sigma }^2$ we estimate it as ${\hat{\sigma }}^2={\mathbb{E}}_n[{\hat{\epsilon }}_i^2]$.

The terms $x_i{\epsilon }_i$ are called scores and we can already see their central importance for inference.

Heteroskedastic Error Term

Assumption: $\mathbb{E}[{\epsilon }_i{x}_i{\epsilon }_j^{\prime }{x}_j^{\prime }]=0$, for all $j\ne i$ and $\mathbb{E}[{\epsilon }_i^4]\le \mathrm{\infty }$, $\mathbb{E}[||x_i|{|}^4]\le C<\mathrm{\infty }$ a.s.

Theorem: Under GM assumptions (1)-(4) plus heteroskedastic error term, the following estimators are consistent, i.e. $\hat{\mathrm{\Sigma }}\stackrel{p}{\to }\mathrm{\Sigma }$.

Note that we are only looking at $\mathrm{\Omega }$ of the $\mathrm{\Sigma }=Q^{-1}\mathrm{\Omega }{Q}^{-1}$ matrix.

• HC0: use the observed residual ${\hat{\epsilon }}_i$ $${\mathrm{\Omega }}_{HC0}={\mathbb{E}}_n[x_i{x}_i^{\prime }{\hat{\epsilon }}_i^2]$$ When $k$ is too big relative to $n$ – i.e., $k/n\to c>0$ – ${\hat{\epsilon }}_i^2$ are too small (${\mathrm{\Omega }}_{HC0}$ biased towards zero). ${\mathrm{\Omega }}_{HC1}$, ${\mathrm{\Omega }}_{HC2}$ and ${\mathrm{\Omega }}_{HC3}$ try to correct this small sample bias.
• HC1: degree of freedom correction (default robust in Stata) $${\mathrm{\Omega }}_{HC1}=\frac{1}{n-k}{\mathbb{E}}_n[x_i{x}_i^{\prime }{\hat{\epsilon }}_i^2]$$
• HC2: use standardized residuals $${\mathrm{\Omega }}_{HC2}={\mathbb{E}}_n[x_i{x}_i^{\prime }{\hat{\epsilon }}_i^2(1-h_{ii}{)}^{-1}]$$ where $h_{ii}=[X(X^{\prime }X{)}^{-1}{X}^{\prime }{]}_{ii}$ is the leverage of the $i^{th}$ observation. A large $h_{ii}$ means that observation $i$ is unusual in the sense that the regressor $x_i$ is far from its sample mean.
• HC3: use prediction error, equivalent to Jack-knife estimator, i.e., ${\mathbb{E}}_n[x_i{x}_i^{\prime }{\hat{\epsilon }}_{(-i)}^2]$ $${\mathrm{\Omega }}_{HC3}={\mathbb{E}}_n[x_i{x}_i^{\prime }{\hat{\epsilon }}_i^2(1-h_{ii}{)}^{-2}]$$ This estimator does not overfit when $k$ is relatively big with respect to $n$. Idea: you exclude the corresponding observation when estimating a particular ${\epsilon }_i$: ${\hat{\epsilon }}_i=y_i-x_i^{\prime }{\hat{\beta }}_{-i}$.

HC0 Consistency

Theorem

Under regularity conditions HC0 is consistent, i.e. ${\hat{\mathrm{\Sigma }}}_{HC0}\stackrel{p}{\to }\mathrm{\Sigma }$. $$\hat{\mathrm{\Sigma }}={\hat{Q}}^{-1}\hat{\mathrm{\Omega }}{\hat{Q}}^{-1}\stackrel{p}{\to }\mathrm{\Sigma }\text{ with }\hat{\mathrm{\Omega }}={\mathbb{E}}_n[x_i{x}_i^{\prime }{\hat{\epsilon }}_i^2]\text{ and }\hat{Q}={\mathbb{E}}_n[x_i{x}_i^{\prime }{]}^{-1}$$

Why is the proof relevant? You cannot directly apply the WLLN to $\hat{\mathrm{\Sigma }}$.

Proof

For the case $\dim (x_i)=1$.

1. ${\hat{Q}}^{-1}\stackrel{p}{\to }{Q}^{-1}$ by WLLN since $x_i$ is iid, $\mathbb{E}[x_i^4]<\mathrm{\infty }$
2. $\overline{\mathrm{\Omega }}={\mathbb{E}}_n[{\epsilon }_i^2{x}_i{x}_i^{\prime }]\stackrel{p}{\to }\mathrm{\Omega }$ by WLLN since ${\mathbb{E}}_n[{\epsilon }_i^4]<c$ and $x_i$ bounded.
3. By the triangle inequality, $$|\hat{\mathrm{\Omega }}-\hat{\mathrm{\Omega }}|\le \underset{\stackrel{p}{\to }0}{\underbrace{|\mathrm{\Omega }-\overline{\mathrm{\Omega }}|}}+\underset{\text{WTS:}\stackrel{p}{\to }0}{\underbrace{|\overline{\mathrm{\Omega }}-\hat{\mathrm{\Omega }}|}}$$
4. We want to show $|\overline{\mathrm{\Omega }}-\hat{\mathrm{\Omega }}|\stackrel{p}{\to }0$ $$\begin{array}{rl}|\overline{\mathrm{\Omega }}-\hat{\mathrm{\Omega }}|& ={\mathbb{E}}_n[{\epsilon }_i^2{x}_i^2]-{\mathbb{E}}_n[{\hat{\epsilon }}_i^2{x}_i^2]=\\ & ={\mathbb{E}}_n[({\epsilon }_i^2-{\hat{\epsilon }}_i^2)x_i^2]\le \\ & \le {\mathbb{E}}_n{\left[{({\epsilon }_i^2-{\hat{\epsilon }}_i^2)}^2\right]}^{\frac{1}{2}}{\mathbb{E}}_n[x_i^4{]}^{\frac{1}{2}}\end{array}$$ where ${\mathbb{E}}_n[x_i^4{]}^{\frac{1}{2}}\stackrel{p}{\to }\mathbb{E}[x_i^4{]}^{\frac{1}{2}}$ by $x_i$ bounded, iid and CMT.
5. We want to show that ${\mathbb{E}}_n\left[{({\epsilon }_i^2-{\hat{\epsilon }}_i^2)}^2\right]\le \eta$ with $\eta \to 0$. Let $L=\underset{i}{max}|{\hat{\epsilon }}_i-{\epsilon }_i|$ (RV depending on $n$), with $L\stackrel{p}{\to }0$ since $$|{\hat{\epsilon }}_i-{\epsilon }_i|=|x_i\hat{\beta }-x_i\beta |\le |x_i||\hat{\beta }-\beta |\stackrel{p}{\to }c\cdot 0$$ We can depompose $$\begin{array}{rl}{({\epsilon }_i^2-{\hat{\epsilon }}_i^2)}^2& ={({\epsilon }_i-{\hat{\epsilon }}_i)}^2{({\epsilon }_i+{\hat{\epsilon }}_i)}^2\le \\ & \le {({\epsilon }_i+{\hat{\epsilon }}_i)}^2{L}^2=\\ & ={(2{\epsilon }_i-{\epsilon }_i+{\hat{\epsilon }}_i)}^2{L}^2\le \\ & \le {(2(2{\epsilon }_i{)}^2+2({\hat{\epsilon }}_i-{\epsilon }_i{)}^2)}^2{L}^2\le \\ & \le (8{\epsilon }_i^2+2L^2)L^2\end{array}$$ Hence $$\mathbb{E}\left[{({\epsilon }_i^2-{\hat{\epsilon }}_i^2)}^2\right]\le L^2(8{\mathbb{E}}_n[{\epsilon }_i^2]+2{\mathbb{E}}_n[L^2])\stackrel{p}{\to }0$$ $$\begin{array}{}◼& \end{array}$$

Heteroskedastic and Autocorrelated Error Term

Assumption

There esists a $\overline{d}$ such that:

• $\mathbb{E}[{\epsilon }_i{x}_i{\epsilon }_{i-d}^{\prime }{x}_{i-d}^{\prime }]\ne 0$ for $d\le \overline{d}$
• $\mathbb{E}[{\epsilon }_i{x}_i{\epsilon }_{i-d}^{\prime }{x}_{i-d}^{\prime }]=0$ for $d>\overline{d}$

Intuition: observations far enough from each other are not correlated.

We can express the variance of the score as $$\begin{array}{rl}{\mathrm{\Omega }}_n& =Var(\sqrt{n}{\mathbb{E}}_n[x_i{\epsilon }_i])=\\ & =\mathbb{E}\left[\left(\frac{1}{n}\sum _{i=1}^n{x}_i{\epsilon }_i\right)\left(\frac{1}{n}\sum _{j=1}^n{x}_j{\epsilon }_j\right)\right]=\\ & =\frac{1}{n}\sum _{i=1}^n\sum _{j=1}^n\mathbb{E}[x_i{\epsilon }_i{x}_j^{\prime }{\epsilon }_j^{\prime }]=\\ & =\frac{1}{n}\sum _{i=1}^n\sum _{j:|i-j|\le \overline{d}}\mathbb{E}[x_i{\epsilon }_i{x}_j^{\prime }{\epsilon }_j^{\prime }]=\\ & =\frac{1}{n}\sum _{d=0}^{\overline{d}}\sum _{i=d}^n\mathbb{E}[x_i{\epsilon }_i{x}_{i-d}^{\prime }{\epsilon }_{i-d}^{\prime }]\end{array}$$

We estimate ${\mathrm{\Omega }}_n$ by $${\hat{\mathrm{\Omega }}}_n=\frac{1}{n}\sum _{d=0}^{\overline{d}}\sum _{i=d}^n{x}_i{\hat{\epsilon }}_i{x}_{i-d}^{\prime }{\hat{\epsilon }}_{i-d}^{\prime }$$

Theorem

If $\overline{d}$ is a fixed integer, then $${\hat{\mathrm{\Omega }}}_n-{\mathrm{\Omega }}_n\stackrel{p}{\to }0$$

What if $\overline{d}$ does not exist (all $x_i,x_j$ are correlated)? $${\hat{\mathrm{\Omega }}}_n=\frac{1}{n}\sum _{d=0}^n\sum _{i=d}^n{x}_i{\hat{\epsilon }}_i{x}_{i-d}^{\prime }{\hat{\epsilon }}_{i-d}^{\prime }=n{\mathbb{E}}_n[x_i{\hat{\epsilon }}_i{]}^2=0$$ By the orthogonality property of the OLS residual.

HAC with Uniform Kernel $${\hat{\mathrm{\Omega }}}_h=\frac{1}{n}\sum _{i,j}{x}_i{\hat{\epsilon }}_i{x}_j^{\prime }{\hat{\epsilon }}_j^{\prime }\mathbb{I}\{|i-j|\le h\}$$ where $h$ is the bandwidth of the kernel. The bandwidth is chosen such that $\mathbb{E}[x_i{\epsilon }_i{x}_{i-d}^{\prime }{\epsilon }_{i-d}^{\prime }]$ is small for $d>h$. How small? Small enough for the estimates to be consistent.

HAC with General Kernel $${\hat{\mathrm{\Omega }}}_{k,h}^{HAC}=\frac{1}{n}\sum _{i,j}{x}_i{\hat{\epsilon }}_i{x}_j^{\prime }{\hat{\epsilon }}_j^{\prime }k\left(\frac{|i-j|}{n}\right)$$

HAC Consistency

Theorem If the joint distribution is stationary and $\alpha$-mixing with $\sum _{k=1}^{\mathrm{\infty }}{k}^2\alpha (k)<\mathrm{\infty }$ and

• $\mathbb{E}[|x_{ij}{\epsilon }_i{|}^{\nu }]<\mathrm{\infty }$ $\mathrm{\forall }\nu$
• ${\hat{\epsilon }}_i=y_i-x_i^{\prime }\hat{\beta }$ for some $\hat{\beta }\stackrel{p}{\to }{\beta }_0$
• $k$ smooth, symmetric, $k(0)\to \mathrm{\infty }$ as $z\to \mathrm{\infty }$, $\int k^2<\mathrm{\infty }$
• $\frac{h}{n}\to 0$
• $h\to \mathrm{\infty }$

Then the HAC estimator is consistent. $${\hat{\mathrm{\Omega }}}_{k,h}^{HAC}-{\mathrm{\Omega }}_n\stackrel{p}{\to }0$$

Comments

We want to choose $h$ small relative to $n$ in order to avoid estimation problems. But we also want to choose $h$ large so that the remainder is small: $$\begin{array}{rl}{\mathrm{\Omega }}_n& =Var(\sqrt{n}{\mathbb{E}}_n[x_i{\epsilon }_i])=\\ & =\underset{{\mathrm{\Omega }}_n^h}{\underbrace{\frac{1}{n}\sum _{i,j:|i-j|\le h}\mathbb{E}[x_i{\epsilon }_i{x}_j^{\prime }{\epsilon }_j^{\prime }]}}+\underset{\text{remainder: }{R}_n}{\underbrace{\frac{1}{n}\sum _{i,j:|i-j|>h}\mathbb{E}[x_i{\epsilon }_i{x}_j^{\prime }{\epsilon }_j^{\prime }]}}=\\ & ={\mathrm{\Omega }}_n^h+R_n\end{array}$$

In particular, HAC theory requires: $${\hat{\mathrm{\Omega }}}^{HAC}\stackrel{p}{\to }\mathrm{\Omega }\text{ if }\{\begin{array}{ll}& \frac{h}{n}\to 0\\ & h\to \mathrm{\infty }\end{array}$$

But in practice, long-run estimation implies $\frac{h}{n}\simeq 0$ which is not ``safe” in the sense that it does not imply $R_n\simeq 0$. On the other hand, if $h\simeq n$, ${\hat{\mathrm{\Omega }}}^{HAC}$ does not converge in probability because it’s too noisy.

Choice of h

How to choose $h$? Look at the score autocorrelation function (ACF).

It looks like after 10 periods the empirical autocorrelation is quite small but still not zero.

Fixed b Asymptotics

[Neave, 1970]: “When proving results on the asymptotic behavior of estimates of the spectrum of a stationary time series, it is invariably assumed that as the sample size $n$ tends to infinity, so does the truncation point $h$, but at a slower rate, so that $\frac{h}{n}$ tends to zero. This is a convenient assumption mathematically in that, in particular, it ensures consistency of the estimates, but it is unrealistic when such results are used as approximations to the finite case where the value of $\frac{h}{n}$ cannot be zero.””

Fixed b Theorem

Theorem

Under regularity conditions, $$\sqrt{n}(V_{k,h}^{HAC})(\hat{\beta }-{\beta }_0)\stackrel{d}{\to }F$$

The asymptotic critical values of the $F$ statistic depend on the choice of the kernel. In order to do hypothesis testing, Kiefer and Vogelsang(2005) provide critical value functions for the t-statistic for each kernel-confidence level combination using a cubic equation: $$cv(b)=a_0+a_{1}b+a_2{b}^2+a_3{b}^3$$

Example

Example for the Bartlett kernel:

Fixed G Asymptotics

[Bester, 2013]: “Cluster covariance estimators are routinely used with data that has a group structure with independence assumed across groups. Typically, inference is conducted in such settings under the assumption that there are a large number of these independent groups.””

“However, with enough weakly dependent data, we show that groups can be chosen by the researcher so that group-level averages are approximately independent. Intuitively, if groups are large enough and well shaped (e.g. do not have gaps), the majority of points in a group will be far from other groups, and hence approximately independent of observations from other groups provided the data are weakly dependent. The key prerequisite for our methods is the researcher’s ability to construct groups whose averages are approximately independent. As we show later, this often requires that the number of groups be kept relatively small, which is why our main results explicitly consider a fixed (small) number of groups.””

Assumption

Assumption Suppose you have data $D=(y_{it},x_{it}{)}_{i=1,t=1}^{N,T}$ where $y_{it}=x_{it}^{\prime }\beta +{\alpha }_i+{\epsilon }_{it}$ where $i$ indexes the observational unit and $t$ indexes time (could also be space).

Let $$\begin{array}{rl}& {\tilde{y}}_{it}=y_{it}-\frac{1}{T}\sum _{t=1}^T{y}_{it}\\ & {\tilde{x}}_{it}=x_{it}-\frac{1}{T}\sum _{t=1}^T{x}_{it}\\ & {\tilde{\epsilon }}_{it}={\epsilon }_{it}-\frac{1}{T}\sum _{t=1}^T{\epsilon }_{it}\end{array}$$ Then $${\tilde{y}}_{it}={\tilde{x}}_{it}^{\prime }\beta +{\tilde{\epsilon }}_{it}$$

The ${\tilde{\epsilon }}_{it}$ are by construction correlated between each other even if the original $\epsilon$ was iid. The cluster score variance estimator is given by: $${\hat{\mathrm{\Omega }}}^{CL}=\frac{1}{T-1}\sum _{i=1}^n\sum _{t=1}^T\sum _{s=1}^T{\tilde{x}}_{it}{\hat{\tilde{\epsilon }}}_{it}{\tilde{x}}_{is}{\hat{\tilde{\epsilon }}}_{is}$$

It’s very similar too the HAC estimator since we have dependent cross-products here as well. However, here we do not consider the $i\times j$ cross-products. We only have time-dependency (state).

Comments (1)

On $T$ and $n$:

• If $T$ is fixed and $n\to \mathrm{\infty }$, then the number of cross-products considered is much smaller than the total number of cross-products.
• If $T»n$ issues arise since the number of cross products considered is close to the total number of cross products. As in HAC estimation, this is a problem because it implies that the algebraic estimate of the cluster score variance gets close to zero because of the orthogonality property of the residuals.
• The panel assumption is that observations across individuals are not correlated.

Strategy: as in HAC, we want to limit the correlation across clusters (individuals). We hope that observations are negligibly dependent between cluster sufficiently distant from each other.

Comments (2)

Classical cluster robust estimator: $${\hat{\mathrm{\Omega }}}^{CL}=\frac{1}{n}\sum _{i=1}^n{x}_i{\epsilon }_i{x}_j^{\prime }{\epsilon }_j^{\prime }\mathbb{I}\{i,j\text{ in the same cluster}\}$$

On clusters:

• If the number of observations near a boundary is small relative to the sample size, ignoring the dependence should not affect inference too adversely.
• The higher the dimension of the data, the easier it is to have observations near boundaries (curse of dimensionality).
• We would like to have few clusters in order to make less independence assumptions. However, few clusters means bigger blocks and hence a larger number of cross-products to estimate. If the number of cross-products is too large (relative to the sample size), ${\hat{\mathrm{\Omega }}}^{CL}$ does not converge

Theorem: Under regularity conditions: $$\hat{t}\stackrel{d}{\to }\sqrt{\frac{G}{G-1}}{t}_{G-1}$$

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*β + ε;

Ideal Estimate

# OLS estimator
β_hat = (X'*X)\(X'*y);

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

# Homoskedastic standard errors
std_h = var(ε_hat) * inv(X'*X);

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

# Leverage
h = diag(P);

HC Estimates

# HC0 variance and standard errors
Ω_hc0 = X' * (I(n) .* ε_hat.^2) * X;
std_hc0 = sqrt.(diag(inv(X'*X) * Ω_hc0 * inv(X'*X)))

## 2-element Array{Float64,1}:
##  0.24691300271914793
##  0.28044707935951835

# HC1 variance and standard errors
Ω_hc1 = n/(n-k) * X' * (I(n) .* ε_hat.^2) * X;
std_hc1 = sqrt.(diag(inv(X'*X) * Ω_hc1 * inv(X'*X)))

## 2-element Array{Float64,1}:
##  0.24941979797977423
##  0.2832943308272532

# HC2 variance and standard errors
Ω_hc2 = X' * (I(n) .* ε_hat.^2 ./ (1 .- h)) * X;
std_hc2 = sqrt.(diag(inv(X'*X) * Ω_hc2 * inv(X'*X)))

## 2-element Array{Float64,1}:
##  0.2506509902982869
##  0.2850878737103963

# HC3 variance and standard errors
Ω_hc3 = X' * (I(n) .* ε_hat.^2 ./ (1 .- h).^2) * X;
std_hc3 = sqrt.(diag(inv(X'*X) * Ω_hc3 * inv(X'*X)))

## 2-element Array{Float64,1}:
##  0.25446321015850176
##  0.2898264779289438

# Note what happens if you allow for full autocorrelation
omega_full = X'*ε_hat*ε_hat'*X;

Inference

Hypothesis Testing

In order to do inference on $\hat{\beta }$ we need to know its distribution. We have two options: (i) assume gaussian error term (extended GM) or (ii) rely on asymptotic approximations (CLT).

A statistical hypothesis is a subset of a statistical model, $\mathcal{K}\subset \mathcal{F}$. A hypothesis test is a map $\mathcal{D}\to \{0,1\}$, $D\mapsto T$. If $\mathcal{F}$ is the statistical model and $\mathcal{K}$ is the statistical hypothesis, we use the notation $H_0:Pr\in \mathcal{K}$.

Generally, we are interested in understanding whether it is likely that data $D$ are drawn from $\mathcal{K}$ or not.

A hypothesis test, $T$ is our tool for deciding whether the hypothesis is consistent with the data. $T(D)=0$ implies fail to reject $H_0$ and test inconclusive $T(D)=1$ $⟹$ reject $H_0$ and $D$ is inconsistent with any $Pr\in \mathcal{K}$.

Let $\mathcal{K}\subseteq \mathcal{F}$ be a statistical hypothesis and $T$ a hypothesis test.

1. Suppose $Pr\in \mathcal{K}$. A Type I error (relative to $Pr$) is an event $T(D)=1$ under $Pr$.
2. Suppose $Pr\in {\mathcal{K}}^c$. A Type II error (relative to $Pr$) is an event $T(D)=0$ under $Pr$.

The corresponding probability of a type I error is called size. The corresponding probability of a type II error is called power (against the alternative $Pr$).

In this section, we are interested in testing three hypotheses, under the assumptions of linearity, strict exogeneity, no multicollinearity, normality on the error term. They are:

1. $H_0:{\beta }_{0k}={\overline{\beta }}_{0k}$ (single coefficient, ${\overline{\beta }}_{0k}\in \mathbb{R}$, $k\le K$)
2. $a^{\prime }{\beta }_0=c$ (linear combination, $a\in {\mathbb{R}}^K,c\in \mathbb{R}$)
3. $R{\beta }_0=r$ (linear restrictions, $R\in {\mathbb{R}}^{p\times K}$, full rank, $r\in {\mathbb{R}}^p$)

Testing Problem

Consider the testing problem $H_0:{\beta }_{0k}={\overline{\beta }}_{0k}$ where ${\overline{\beta }}_{0k}$ is a pre-specified value under the null. The t-statistic for this problem is defined by $$t_k:=\frac{b_k-{\overline{\beta }}_{0k}}{SE(b_k)},SE(b_k):=\sqrt{s^2[(X^{\prime }X{)}^{-1}{]}_{kk}}$$

Theorem: In the testing procedure above, the sampling distribution under the null $H_0$ is given by $$t_k|X\sim t_{n-k}\text{and so}{t}_k\sim t_{n-k}$$

$t_{(n-K)}$ denotes the t-distribution with $(n-k)$ degress of freedom. The test can be one sided or two sided. The above sampling distribution can be used to construct a confidence interval.

Example

We want to asses whether or not the ``true” coefficient ${\beta }_0$ equals a specific value $\hat{\beta }$. Specifically, we are interested in testing $H_0$ against $H_1$, where:

• Null Hypothesis: $H_0:{\beta }_0=\hat{\beta }$
• Alternative Hypothesis: $H_1:{\beta }_0\ne \hat{\beta }$.

Hence, we are interested in a statistic informative about $H_1$, which is the Wald test statistic $$|T^{\ast }|=|\frac{\hat{\beta }-{\beta }_0}{\sigma (\hat{\beta })}|\sim N(0,1)$$

However, the true variance ${\sigma }^2(\hat{\beta })$ is not known and has to be estimated. Therefore we plug in the sample variance ${\hat{\sigma }}^2(\hat{\beta })=\frac{n}{n-1}{\mathbb{E}}_n[{\hat{e}}_i^2]$ and we use $$|T|=|\frac{\hat{\beta }-{\beta }_0}{\hat{\sigma }(\hat{\beta })}|\sim t_{(n-k)}$$

Comments

Hypothesis testing is like proof by contradiction. Imagine the sampling distribution was generated by $\beta$. If it is highly improbable to observe $\hat{\beta }$ given ${\beta }_0=\beta$ then we reject the hypothesis that the sampling distribution was generated by $\beta$.

Then, given a realized value of the statistic $|T|$, we take the following decision:

• Do not reject $H_0$: it is consistent with random variation under true $H_0$—i.e., $|T|$ small as it has an exact student t distribution with $(n-k)$ degree of freedom in the normal regression model.
• Reject $H_0$ in favor of $H_1$: $|T|>c$, with $c$ being the critical values selected to control for false rejections: $Pr(|t_{n-k}|\ge c)=\alpha$. Moreover, you can also reject $H_0$ if the p-value $p$ is such that: $p<\alpha$.

Comments (2)

The probability of false rejection is decreasing in $c$, i.e. the critical value for a given significant level. $$\begin{array}{rl}Pr(\text{Reject }{H}_0|H_0)& =Pr(|T|>c|H_0)=\\ & =Pr(T>c|H_0)+Pr(T<-c|H_0)=\\ & =1-F(c)+F(-c)=2(1-F(c))\end{array}$$

Example: Consider the testing problem $H_0:a^{\prime }{\beta }_0=c$ where $a$ is a pre-specified linear combination under study. The t-statistic for this problem is defined by: $$t_k:=\frac{a^{\prime }b-c}{SE(a^{\prime }b)},SE(a^{\prime }b):=\sqrt{s^2{a}^{\prime }(X^{\prime }X{)}^{-1}a}$$

t Stat

Theorem

In the testing procedure above, the sampling distribution under the null $H_0$ is given by $$t_a|X\sim t_{n-K}\text{and so}{t}_a\sim t_{n-K}$$

Like in the previous test, $t_{(n-K)}$ denotes the t-distribution with $(n-K)$ degress of freedom. The test can again be one sided or two sided. The above sampling distribution can be used to construct a confidence interval

F Stat

Example

Consider the testing problem $$H_0:R{\beta }_0=r$$ where $R\in {\mathbb{R}}^{p\times k}$ is a presepecified set of linear combinations and $r\in {\mathbb{R}}^p$ is a restriction vector.

The F-statistic for this problem is given by $$F:=\frac{(Rb-r{)}^{\prime }[R(X^{\prime }X)R^{\prime }{]}^{-1}(Rb-r)/p}{s^2}$$

Theorem

For the problem, the sampling distribution of the F-statistic under the null $H_0:$ $$F|X\sim F_{p,n-K}\text{and so}F\sim F_{p,n-K}$$

The test is intrinsically two-sided. The above sampling distribution can be used to construct a confidence interval.

Equivalence

Theorem

Consider the testing problem $H_0:R{\beta }_0=r$ where $R\in {\mathbb{R}}^{p\times K}$ is a presepecified set of linear combinations and $r\in {\mathbb{R}}^p$ is a restriction vector.

Consider the restricted least squares estimator, denoted ${\hat{\beta }}_R$: ${\hat{\beta }}_R:=\text{arg}\underset{\beta :R\beta =r}{min}Q(\beta )$. Let $SS{R}_U=Q(b),SS{R}_R=Q({\hat{\beta }}_R)$. Then the $F$ statistic is numerically equivalent to the following expression: $F=\frac{(SS{R}_R-SS{R}_U)/p}{SS{R}_U/(n-K)}$.

Confidence Intervals

A confidence interval at $(1-\alpha )$ is a random set $C$ such that $$Pr({\beta }_0\in C)\ge 1-\alpha$$ i.e. the probability that $C$ covers the true value $\beta$ is fixed at $(1-\alpha )$.

Since $C$ is not known, it has to be estimated ($\hat{C}$). We construct confidence intervals such that:

• they are symmetric around $\hat{\beta }$;
• their length is proportional to $\sigma (\hat{\beta })=\sqrt{Var(\hat{\beta })}$.

A CI is equivalent to the set of parameter values such that the t-statistic is less than $c$, i.e., $$\hat{C}=\{\beta :|T(\beta )|\le c\}=\{\beta :-c\le \frac{\beta -\hat{\beta }}{\sigma (\hat{\beta })}\le c\}$$

In practice, to construct a 95% confidence interval for a single coefficient estimate ${\hat{\beta }}_j$, we use the fact that $$Pr(\frac{|{\hat{\beta }}_j-{\beta }_{0,j}|}{\sqrt{{\sigma }^2[(X^{\prime }X{)}^{-1}{]}_{jj}}}>1.96)=0.05$$

Code

# t-test for beta=0
t = abs.(β_hat ./ (std_hc1));

# p-value
p_val = 1 .- cdf.(Normal(0,1), t);

# F statistic of joint significance
SSR_u = ε_hat'*ε_hat;
SSR_r = y'*y;
F = (SSR_r - SSR_u)/k / (SSR_u/(n-k));

# 95# confidente intervals
conf_int = [β_hat - 1.96*std_hc1, β_hat + 1.96*std_hc1];