1. Statistical Foundations: From Samples to Regression

Question: If a fund’s estimated exposure or average excess return is nonzero, how much of that result could be sampling noise?

We build the tools for answering this question: descriptive statistics, sampling distributions, hypothesis tests, ordinary least squares, and covariance estimators that allow unequal variance and serial dependence. Small examples make the algebra checkable by hand; simulations show where asymptotic approximations remain imperfect.

By the end you should be able to derive $\hat\beta=(X^\top X)^{-1}X^\top y$, explain the difference between a coefficient and its standard error, and assemble the HC and Newey–West sandwich covariance estimators.

Background: random variables, expectation and variance, matrix multiplication and inversion. Python verifies the calculations, but the mathematical argument can be read independently.

First-pass route: follow Sections 1–3.2, then hypothesis testing and OLS in Sections 4–5, followed by robust inference in Sections 6–9. The generalised CLTs in Sections 3.3–3.6 explain how the assumptions change; revisit them when studying dependence in Chapter 2.

Book introduction · Next: Financial time series

Section 1: Descriptive Statistics — A Hand-Checkable Example¶

We start with a dataset small enough to compute everything on paper. Suppose we observe 5 monthly excess returns (in percent):

1.1 Sample Mean¶

1.2 Sample Variance (with Bessel’s correction)¶

The deviations from the mean are:

1.3 Why $n - 1$? (Bessel’s Correction)¶

If we knew the true population mean $\mu$, we would divide by $n$. But we estimated $\mu$ from the same data (using $\bar{x}$), which “uses up” one degree of freedom. Dividing by $n - 1$ corrects for this, making $s^2$ an unbiased estimator of the population variance:

With only $n = 5$ observations, the difference between dividing by 4 vs. 5 is 25% — far from negligible!

1.4 Covariance and Correlation (Two Variables)¶

Now suppose we have a second variable — say, the market excess return for the same 5 months:

Sample covariance:

Sample correlation:

We already know $s_x = \sqrt{2.5}$. For $y$: the squared deviations are $1.8^2 + 3.2^2 + 2.8^2 + 0.2^2 + 1.2^2 = 3.24 + 10.24 + 7.84 + 0.04 + 1.44 = 22.80$, so $s_y^2 = 22.80/4 = 5.70$ and $s_y = \sqrt{5.70} \approx 2.3875$.

This strong positive correlation ($\approx 0.93$) means $x$ and $y$ move together closely.

Section 2: Random Variables, Expectation, and the Population–Sample Distinction¶

2.1 Population vs. Sample¶

A crucial distinction underlies all of statistics:

• The population is the entire (usually infinite or very large) set of possible outcomes. For example, “all possible monthly returns that SPY could ever produce.”

• A sample is a finite collection of observations drawn from the population. For example, “the 168 monthly returns we actually observed from 2010–2023.”

We use sample statistics ($\bar{x}$, $s^2$, etc.) to estimate population parameters ($\mu$, $\sigma^2$, etc.). The sample mean $\bar{x}$ is our best guess of the true mean $\mu$, but it comes with uncertainty because different samples would give different $\bar{x}$ values.

2.2 The Standard Error of the Mean¶

If $X_1, X_2, \ldots, X_n$ are independent draws from a population with mean $\mu$ and variance $\sigma^2$, then the sample mean

has:

Why $n$ and not $n - 1$ here? In Section 1.3 we used Bessel’s correction ($n - 1$) when estimating the population variance $\sigma^2$ from a sample — that corrects for the bias introduced by using $\bar{x}$ in place of the unknown $\mu$. Here the situation is different: we are not estimating anything. The formula $\text{Var}(\bar{X}) = \sigma^2 / n$ is an exact algebraic result that follows directly from the properties of variance. Since the $X_i$ are independent, $\text{Var}\!\left(\sum X_i\right) = \sum \text{Var}(X_i) = n\sigma^2$, and dividing a random variable by $n$ scales its variance by $1/n^2$, giving $\text{Var}(\bar{X}) = n\sigma^2 / n^2 = \sigma^2 / n$. The $n$ here simply counts how many observations are being averaged — it has nothing to do with degrees of freedom or bias correction.

The standard error of the mean is therefore:

In practice $\sigma$ is unknown, so we replace it with the sample standard deviation $s$ (computed with the Bessel-corrected $n - 1$ denominator from Section 1.2). Note that the $n - 1$ inside $s$ and the $n$ inside $\sqrt{n}$ play completely different roles: $n - 1$ corrects the bias when estimating the variance, while $n$ reflects the number of observations being averaged. The estimated standard error is:

Hand calculation with our data ($\bar{x} = 1$, $s = \sqrt{2.5}$, $n = 5$):

This tells us: “If we repeatedly drew samples of size 5 from this population and computed the mean each time, those means would have a standard deviation of about 0.71.”

Section 3: The Central Limit Theorem — The Engine of Inference¶

The CLT is the single most important result in statistics. It is the reason we can do hypothesis testing, build confidence intervals, and compute p-values. In this section we state the classical version first, then progressively relax its assumptions toward versions that actually apply to financial time series. By the end, you will understand exactly which CLT variant underwrites the Newey-West standard errors used throughout the Fama-French tutorials.

Reading guide. Sections 3.1–3.2 (classical CLT) and Section 3.7 (visual demo) are essential for everyone. Sections 3.3–3.5 develop the theoretical machinery that justifies Newey-West for time-series data — they can be skimmed on a first reading and revisited after working through Section 7.

What are Newey-West standard errors? In short, they are a method for computing standard errors that remain valid even when the data are serially correlated and/or heteroscedastic — two features that are ubiquitous in financial time series. Where ordinary standard errors assume each observation is independent, Newey-West standard errors account for the fact that nearby observations may be correlated, allowing intervals to become wider or narrower depending on the covariance structure. We develop the full theory and a hand-checkable example in Section 7.

3.1 Version 1: The Classical CLT (Lindeberg-Lévy)¶

Theorem. Let $X_1, X_2, \ldots, X_n$ be independent and identically distributed (i.i.d.) random variables with mean $\mu$ and finite variance $\sigma^2 > 0$. Then as $n \to \infty$:

In words: the standardized sample mean converges in distribution to a standard normal, regardless of the shape of the underlying distribution.

Assumptions: (1) identical distributions, (2) independence, (3) finite variance.

What this buys us: We can compute the distribution of $\bar{X}_n$ without knowing the shape of the $X_i$. This is extraordinary — we need no parametric model for the data.

Relevance for finance: If monthly returns were truly i.i.d. (same distribution every month, no dependence across time), this version would be all we need. But they are not — so we need to relax the assumptions.

3.2 Why This Is Revolutionary¶

The CLT tells us something extraordinary:

• The $X_i$ might follow any distribution — uniform, exponential, chi-squared, a bizarre bimodal shape, anything with finite variance.

• Yet the average of many such observations will be approximately normal.

• The approximation improves as $n$ grows.

This is why the normal distribution appears everywhere in statistics: it governs the behavior of averages and sums, not necessarily of individual observations.

3.3 Version 2: Dropping “Identical” — The Lindeberg-Feller CLT¶

The problem with “identically distributed.” In a regression $y_i = \alpha + \beta x_i + \epsilon_i$, the coefficient estimator $\hat{\beta} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$ is a weighted average of the $y_i$ — where the weights depend on the $x_i$ values and differ across observations. Even if the errors $\epsilon_i$ are i.i.d., the weighted terms are not identically distributed. The classical CLT does not directly apply.

The solution. The Lindeberg-Feller CLT drops the “identically distributed” requirement:

Theorem (Lindeberg-Feller). Let $X_1, \ldots, X_n$ be independent (but not necessarily identically distributed) random variables with $E[X_i] = \mu_i$ and $\text{Var}(X_i) = \sigma_i^2$. If the Lindeberg tail condition holds — for every $\eta>0$, $s_n^{-2}\sum_i E[(X_i-\mu_i)^2\mathbf{1}\{|X_i-\mu_i|>\eta s_n\}]\to0$, where $s_n^2=\sum_i\sigma_i^2$ — then:

Assumptions: (1) independence, (2) finite variances, (3) the Lindeberg tail condition (variance non-dominance alone is insufficient).

What this buys us: Regression coefficients $\hat{\beta}$ are asymptotically normal even when the regressors $x_i$ take different values across observations. This is why OLS inference works even for regression (not just sample means).

Relevance for finance: This handles the “weighted average” structure of $\hat{\beta}$, but still requires independence — which financial time series violate.

3.4 Version 3: Dropping “Independent” (Partially) — The Martingale Difference CLT¶

The problem with “independent.” Monthly stock returns exhibit volatility clustering — periods of high volatility (like 2008–2009 or March 2020) alternate with calm periods. Even if returns themselves are nearly unpredictable (weak autocorrelation), they are not independent, because knowing yesterday’s return magnitude tells you something about today’s volatility.

A useful middle ground. Many financial models assume that returns are a martingale difference sequence (MDS). This means:

In words: the conditional mean of the next return is zero (given past data), even though other aspects of the conditional distribution (like the variance) may depend on the past. This is weaker than independence — it says returns are unpredictable in conditional mean (including nonlinear functions of past information) but not necessarily independent.

Theorem (MDS-CLT). Let $\{X_t\}$ be a stationary, ergodic martingale difference sequence with $E[X_t^2] = \sigma^2 < \infty$. Then:

Assumptions: (1) MDS property (unpredictability), (2) stationarity and ergodicity, (3) finite variance. Independence is not required.

What this buys us: We can handle volatility clustering. Even though $\text{Var}(X_t \mid \text{past})$ changes over time (GARCH effects), as long as returns themselves are unpredictable, $\bar{X}_n$ is asymptotically normal with variance $\sigma^2/n$ where $\sigma^2 = E[X_t^2]$ is the unconditional variance.

Relevance for finance: An MDS can model demeaned innovations. Market efficiency by itself does not imply zero or constant conditional expected raw returns; risk premia can vary. This CLT version is what makes standard hypothesis tests on mean returns valid — even with GARCH-type volatility clustering — as long as the standard errors are computed correctly (i.e., using the unconditional variance, or better yet, a HAC estimator).

3.5 Version 4: Allowing Mild Dependence — The CLT for Mixing Processes¶

The problem with MDS. In reality, monthly returns may exhibit mild autocorrelation (short-run momentum or reversal effects), and regression residuals may be serially correlated. If $E[X_t \mid X_{t-1}] \neq 0$, we can’t use the MDS-CLT.

The solution. One useful class of CLTs for finance allows weak dependence — observations can be correlated, but observations far apart in time must be nearly independent. This is formalized by mixing conditions (such as $\alpha$-mixing or strong mixing):

Theorem (CLT for stationary mixing sequences). One sufficient version assumes strict stationarity, $E|X_t-\mu|^{2+\delta}<\infty$ for some $\delta>0$, and strong-mixing coefficients satisfying $\sum_{k\ge1}\alpha(k)^{\delta/(2+\delta)}<\infty$, with positive long-run variance. Write $\gamma_k=\operatorname{Cov}(X_t,X_{t+k})$. Then:

The asymptotic variance is not simply $\gamma_0 = \sigma^2$. It is the long-run variance:

This equals $2\pi$ times the spectral density at frequency zero — the “DC component” of the time series.

Assumptions: stationarity, a sufficiently fast mixing rate and suitable moments. Summable autocovariances alone do not imply mixing or a CLT. In regressions, these conditions concern the score process $x_t\epsilon_t$ and the regressor moments, not only the residuals.

What this buys us: A useful sufficient framework for weakly dependent time series. The sample mean $\bar{X}_n$ is still asymptotically normal, but its variance is $\sigma_{\text{LR}}^2 / n$, not $\sigma^2 / n$. If autocovariances are positive, $\sigma_{\text{LR}}^2 > \sigma^2$ and the “naive” standard error $s/\sqrt{n}$ understates the true uncertainty.

Relevance for finance: This is the CLT that Newey-West standard errors are built on. Newey-West estimates $\sigma_{\text{LR}}^2$ from the data using a weighted sum of sample autocovariances. We will derive this in full detail in Section 7.

3.6 Summary: Which CLT Do We Need?¶

For empirical asset pricing with financial time series, we need the Mixing CLT (or at minimum the MDS-CLT). The classical i.i.d. version is a pedagogical starting point — useful for intuition — but it does not match the properties of real financial data.

The key takeaway: When its assumptions hold, a CLT delivers asymptotic normality of the sample mean (or regression coefficients). What changes across versions is the variance expression. Getting the variance right is the entire game — and it is exactly what robust standard errors (Newey-West) are designed to do.

3.7 A Visual Demonstration¶

Let us verify the CLT by repeatedly drawing samples from a decidedly non-normal distribution — the exponential distribution (which is right-skewed with skewness = 2) — and examining the distribution of sample means.

3.8 Why the Long-Run Variance Matters — A Quick Demonstration¶

The classical CLT says $\text{Var}(\bar{X}_n) = \sigma^2 / n$. But when data are autocorrelated, the true variance of the sample mean is $\sigma_{\text{LR}}^2 / n$, where $\sigma_{\text{LR}}^2 = \sum_{k=-\infty}^{\infty} \gamma_k$ can be much larger than $\sigma^2 = \gamma_0$.

The following simulation makes this concrete. We draw from an AR(1) process with autocorrelation $\rho = 0.5$ and compare the actual dispersion of sample means against what the naive formula $\sigma / \sqrt{n}$ predicts.

Section 4: Hypothesis Testing, t-Statistics, and p-Values — and the CLT Connection¶

4.1 The Framework¶

Hypothesis testing asks: “Could this result have arisen by chance?”

Suppose we observe a sample mean $\bar{x} = 1.0$ for a set of monthly excess returns. That sounds like a positive return — but with only 5 noisy observations, maybe the true mean is actually zero and we just got an unlucky draw. Hypothesis testing gives us a principled way to weigh the evidence.

We set up two competing hypotheses:

• $H_0$ (null hypothesis): The effect is zero. For example, $\mu = 0$ (“the true mean excess return is zero — the asset earns no premium”).

• $H_1$ (alternative hypothesis): The effect is nonzero. For example, $\mu \neq 0$.

The logic compares the statistic with its null distribution, allowing a controlled probability of false rejection: we temporarily assume $H_0$ is true and ask how surprising our observed data would be under that assumption. If the data would be very surprising (low probability), we reject $H_0$. If the data is not particularly unusual under $H_0$, we have no grounds to reject it.

This is analogous to a courtroom trial: the null hypothesis is “innocent until proven guilty.” We need sufficient evidence (data) to reject this default presumption. Importantly, failing to reject $H_0$ does not prove $H_0$ is true — it merely means we lack enough evidence to overturn it.

4.2 The t-Distribution — What It Is and Where It Comes From¶

Before constructing our test statistic, we need to understand a distribution that plays a central role in small-sample inference: the Student’s $t$-distribution.

The problem. If $X_1, \ldots, X_n$ are i.i.d. $\mathcal{N}(\mu, \sigma^2)$ and we knew the true $\sigma$, then

But we almost never know $\sigma$. When we replace it with the sample standard deviation $s$, we introduce extra randomness in the denominator — $s$ itself varies from sample to sample. The resulting ratio is no longer standard normal.

The solution (Student, 1908). William Sealy Gosset (writing under the pseudonym “Student”) showed that if the underlying data are normally distributed, then

follows a $t$-distribution with $n - 1$ degrees of freedom. The parameter $n - 1$ (degrees of freedom) matches the denominator in Bessel’s correction, for the same reason: we used up one degree of freedom by estimating $\mu$ with $\bar{x}$.

Key properties of the $t$-distribution:

The heavier tails reflect the extra uncertainty from estimating $\sigma$ with $s$. With small samples, $s$ can differ substantially from $\sigma$, so the $t$-distribution is noticeably wider than the normal. With large samples, $s \approx \sigma$ and the $t$-distribution is nearly indistinguishable from $\mathcal{N}(0,1)$.

For our data with $n = 5$, we have $\nu = 4$ degrees of freedom. The $t_4$ distribution has variance $4/(4-2) = 2$ — double that of the standard normal. This matters: using the normal instead of the $t$ would make us overconfident.

The plot below visualizes exactly this — notice how the $t_4$ distribution spreads more probability into the tails compared to $\mathcal{N}(0,1)$, and how larger $\nu$ brings the $t$-distribution closer to the normal.

4.3 The t-Statistic¶

Under $H_0: \mu = 0$, the test statistic is:

This measures how many standard errors the sample mean is away from zero. The intuition is simple:

• If $\bar{x}$ is far from 0 (the null value) relative to its own uncertainty, that is strong evidence against $H_0$.

• If $\bar{x}$ is close to 0 relative to the noise, the data are consistent with $H_0$.

The $t$-statistic standardizes the sample mean into a “signal-to-noise ratio”: the numerator is the signal ($\bar{x} - 0$) and the denominator is the noise ($s/\sqrt{n}$).

Hand calculation with our data:

So our observed mean is about 1.41 standard errors above zero — present, but not overwhelmingly so.

4.4 The p-Value¶

The p-value is the probability of observing a test statistic at least as extreme as the one we got, assuming $H_0$ is true. “At least as extreme” means at least as far from zero in either direction (for a two-sided test).

Formally, for a two-sided test with $n - 1 = 4$ degrees of freedom:

where $T_4$ denotes a random variable following a $t$-distribution with 4 degrees of freedom.

How to interpret the p-value:

• A small p-value (e.g., $p < 0.05$) means: “Data this extreme would be very unlikely if $H_0$ were true.” This is evidence against $H_0$.

• A large p-value means: “Data this extreme is not unusual under $H_0$.” We have no strong reason to reject $H_0$.

• The p-value is not the probability that $H_0$ is true. It is the probability of seeing data this extreme given that $H_0$ is true — a crucial distinction.

The conventional significance level is $\alpha = 0.05$, meaning we reject $H_0$ if $p < 0.05$. This corresponds to a $|t|$ threshold of about 2 (more precisely, $t_{0.025, \nu}$ for the relevant degrees of freedom). The “rule of thumb” you often hear — “a $t$-statistic above 2 is significant” — comes from this.

4.5 Why a t-Distribution Instead of the Standard Normal?¶

When the sample size $n$ is small and we estimate $\sigma$ with $s$, the ratio $\bar{x} / (s/\sqrt{n})$ does not exactly follow $\mathcal{N}(0,1)$. As derived in Section 4.2, it follows a $t$-distribution with $n - 1$ degrees of freedom under i.i.d. normal sampling. Small sample size alone does not justify that distribution.

Using the standard normal when you should use the $t$-distribution would understate the p-value — you’d think data is more significant than it actually is. This is because the $t$-distribution allocates more probability to the tails.

For example, $P(|Z| \geq 2) = 0.046$ under the standard normal, but $P(|T_4| \geq 2) = 0.116$ under $t_4$. That is a factor of 2.5 difference. With 4 degrees of freedom, the normal would mislead you substantially.

As $n \to \infty$, the $t_{n-1}$ distribution converges to $\mathcal{N}(0,1)$. For $n \geq 30$, the difference is already small. This convergence is related to the CLT: with enough data, estimating $\sigma$ with $s$ introduces negligible additional uncertainty, and the extra randomness in the denominator vanishes.

4.6 Confidence Intervals¶

A confidence interval inverts the hypothesis test. Instead of asking “Is $\mu = 0$?”, it asks: “What range of $\mu$ values are consistent with our data?”

A 95% confidence interval for $\mu$ is:

where $t_{0.025,\, 4} \approx 2.776$ is the 97.5th percentile of the $t_4$ distribution. We use the $t$-distribution (not the normal) for the same reason as above: we are estimating $\sigma$ from a small sample.

The critical value 2.776 is notably larger than the 1.96 we would use from the standard normal. This makes the confidence interval wider, appropriately reflecting our greater uncertainty with $n = 5$.

Hand calculation:

Since this interval contains zero, we cannot reject $H_0: \mu = 0$ at the 5% level — consistent with our p-value being above 0.05.

Interpretation. “If we repeated this experiment many times and computed a 95% CI each time, about 95% of those intervals would contain the true $\mu$.” It does not mean there is a 95% probability that $\mu$ lies in this particular interval — $\mu$ is a fixed (unknown) number, not a random variable.

4.7 Why the CLT Justifies Normal p-Values¶

Let us pause and appreciate the deep connection between the CLT and hypothesis testing. It is the intellectual heart of this entire notebook and fundamental to understanding how Newey-West works.

Above we used the t-distribution to compute the p-value. This was exact — it required that the $X_i$ are normally distributed (or that we use an exact finite-sample result). But in practice, financial returns are not normally distributed: they have fat tails, skewness, and time-varying volatility.

So why do our p-values still work?

The answer is the CLT. Here is the argument, spelled out carefully:

4.8 The Argument in Five Steps¶

Step 1. We want to test $H_0: \mu = 0$. Our test statistic is:

Step 2. By the CLT, for large $n$:

Under $H_0$ ($\mu = 0$): $\frac{\bar{X}_n}{\sigma / \sqrt{n}} \xrightarrow{d} \mathcal{N}(0,1)$.

Step 3. By the Law of Large Numbers, $s^2 \xrightarrow{p} \sigma^2$, i.e. the sample variance converges in probability to the true variance.

Step 4. By Slutsky’s theorem, replacing $\sigma$ with $s$ in the denominator does not change the limiting distribution.

Slutsky’s Theorem. If $X_n \xrightarrow{d} X$ and $Y_n \xrightarrow{p} c$ (a constant), then $X_n / Y_n \xrightarrow{d} X / c$, and more generally $g(X_n, Y_n) \xrightarrow{d} g(X, c)$ for any continuous $g$.

Applying this with $X_n = \bar{X}_n / (\sigma/\sqrt{n}) \xrightarrow{d} \mathcal{N}(0,1)$ and $Y_n = s/\sigma \xrightarrow{p} 1$:

Step 5. Therefore, for large $n$, we can compute p-values using the standard normal $\mathcal{N}(0,1)$ table (or the $t_{n-1}$ distribution, which is nearly identical for large $n$).

4.9 What “Nonparametric” Means Here¶

Notice what we did not assume in Steps 1–5:

• We did not assume the $X_i$ are normally distributed.

• We did not assume the $X_i$ have any particular distributional shape.

• We only assumed they have a finite mean and variance and are (approximately) independent.

This is what makes the CLT-based inference nonparametric (more precisely, semiparametric; it is not an exact finite-sample distribution-free test): it works regardless of the population distribution, provided $n$ is large enough.

4.10 How Large is “Large Enough”?¶

The speed of CLT convergence depends on the distribution’s skewness and kurtosis:

The table is only heuristic: no fixed observation count guarantees a good approximation. Heavy tails, persistent scores and high leverage can make inference inaccurate even with hundreds of observations. This is why finance researchers can use standard-normal-based p-values with Newey-West standard errors even though returns are not normally distributed.

4.11 A Subtle but Critical Point¶

The CLT tells us that $\bar{X}_n$ is approximately normal. This does not mean:

• ❌ Individual observations $X_i$ are normal

• ❌ Residuals from a regression are normal

• ❌ The error term $\epsilon$ is normal

It does mean:

• ✅ This works even when the underlying data is decidedly non-normal

• ✅ Functions of averages (like regression coefficients $\hat{\beta}$, which are weighted averages of the data) are approximately normal for large $n$

• ✅ We can use the normal distribution to compute p-values for these statistics

Section 5: Linear Regression and OLS¶

We now shift from testing a single mean to analyzing the relationship between variables. This is the core of empirical finance: the Fama-French model asks whether exposure to market, size, and value factors explains the cross-section of expected returns — a question that is fundamentally about regression.

5.1 What Is Regression?¶

Regression answers the question: how does one variable change, on average, when another variable changes?

Concretely, given paired observations $(x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n)$, we model:

where:

• $\alpha$ (the intercept) is the expected value of $y$ when $x = 0$,

• $\beta$ (the slope) is the expected change in $y$ per one-unit change in $x$,

• $\epsilon_i$ (the error or residual) captures everything about $y_i$ that is not explained by $x_i$.

In the Fama-French context, $y_i$ is the excess return of a portfolio, $x_i$ is the excess return of the market (or a factor), and $\beta$ measures the portfolio’s sensitivity (or loading) on that factor. A positive $\beta$ on the market factor means the portfolio tends to go up when the market goes up, and vice versa. The intercept $\alpha$ — often called Jensen’s alpha — measures the average return that is not explained by factor exposure. If $\alpha > 0$, the portfolio earns a premium beyond what its risk exposure would predict.

5.2 The Idea: Minimizing Squared Errors¶

How do we choose $\alpha$ and $\beta$? The Ordinary Least Squares (OLS) principle says: choose the line that makes the residuals $\hat{\epsilon}_i = y_i - (\hat{\alpha} + \hat{\beta} x_i)$ as small as possible, in the sense of minimizing the sum of squared residuals:

Why squared errors rather than absolute errors? Three reasons:

1. Differentiability — Squared errors are smooth, so we can use calculus to find the minimum.

2. Unique solution — The squared-error problem always has exactly one solution (for non-degenerate data).

3. Connection to the CLT — Under standard assumptions, the OLS estimator is the best linear unbiased estimator (BLUE), and its sampling distribution is governed by the CLT.

5.3 Deriving the OLS Formulas (Scalar Case)¶

To minimize $\text{RSS}$, we take partial derivatives with respect to $\alpha$ and $\beta$ and set them to zero.

Derivative with respect to $\alpha$:

This yields the first normal equation:

Derivative with respect to $\beta$:

Substituting $\hat{\alpha} = \bar{y} - \hat{\beta}\bar{x}$ and simplifying gives the second normal equation:

This is an elegant result: the slope equals the sample covariance of $x$ and $y$ divided by the sample variance of $x$. This makes intuitive sense — $\hat{\beta}$ measures how much $y$ moves per unit of $x$-variability.

5.4 A Geometric Intuition¶

Think of the $n$ observed $y$-values as a point in $n$-dimensional space. The set of all vectors of the form $\hat{\alpha}\mathbf{1} + \hat{\beta}\mathbf{x}$ (where $\mathbf{1}$ is the column of ones) traces out a 2-dimensional plane within that $n$-dimensional space.

OLS finds the point on that plane that is closest to $\mathbf{y}$ — that is, the orthogonal projection of $\mathbf{y}$ onto the column space of $\mathbf{X}$. The residual vector $\hat{\boldsymbol{\epsilon}} = \mathbf{y} - \hat{\mathbf{y}}$ is perpendicular to the column space, which is precisely the condition $\mathbf{X}^T\hat{\boldsymbol{\epsilon}} = \mathbf{0}$ — the normal equations in matrix form.

This geometric picture explains several properties of OLS that can otherwise seem mysterious:

• The residuals sum to zero (because $\hat{\boldsymbol{\epsilon}}$ is perpendicular to the constant column $\mathbf{1}$).

• The fitted values and residuals are uncorrelated (because $\hat{\mathbf{y}}$ lies in the column space while $\hat{\boldsymbol{\epsilon}}$ is perpendicular to it).

• With an intercept, $R^2$ measures the cosine squared of the angle between the centered vectors $\mathbf{y}-\bar y\mathbf1$ and $\hat{\mathbf{y}}-\bar y\mathbf1$.

5.5 Hand Calculation with Our Data¶

Using the dataset from Section 1: $x = [2, -1, 3, 0, 1]$ and $y = [3, -2, 4, 1, 0]$.

From Section 1 we already computed: $\bar{x} = 1.0$, $\bar{y} = 1.2$.

Step 1: Compute the deviations.

Step 2: Compute $\hat{\beta}$ and $\hat{\alpha}$.

So the fitted line is $\hat{y} = -0.2 + 1.4x$.

Interpretation: Each additional unit of $x$ is associated with a 1.4-unit increase in $y$. The intercept -0.2 is the predicted value of $y$ when $x = 0$.

5.6 Visualizing the Fit¶

The code below plots the data and the fitted regression line, and marks the residuals as vertical distances from each point to the line — making concrete the idea that OLS minimizes the sum of these squared vertical distances.

5.7 The Matrix Formulation¶

The scalar formulas from Section 5.3 work for simple regression (one regressor). But the Fama-French model has three regressors (market, size, value) — and in general we need to handle $k$ regressors simultaneously. Matrix algebra provides a compact, general framework that works for any number of regressors. It also sets up the variance formulas we will need for robust standard errors in subsequent sections.

The model in matrix form. Stack all $n$ observations into a single equation:

where:

• $\mathbf{y}$ is the $n \times 1$ vector of outcomes,

• $\mathbf{X}$ is the $n \times k$ design matrix — its first column is a column of ones (the intercept), and the remaining columns are the regressors,

• $\boldsymbol{\beta}$ is the $k \times 1$ vector of parameters $(\alpha, \beta_1, \ldots, \beta_{k-1})^T$,

• $\boldsymbol{\epsilon}$ is the $n \times 1$ vector of errors.

For our 5-observation simple regression, $k = 2$ (intercept + one slope):

Deriving the OLS estimator. The residual sum of squares in matrix form is:

Expanding:

To minimize, we take the matrix derivative (gradient) and set it to zero. The key matrix calculus identities are:

• $\frac{\partial}{\partial \boldsymbol{\beta}}(\boldsymbol{\beta}^T\mathbf{a}) = \mathbf{a}$, for any constant vector $\mathbf{a}$;

• $\frac{\partial}{\partial \boldsymbol{\beta}}(\boldsymbol{\beta}^T\mathbf{A}\boldsymbol{\beta}) = 2\mathbf{A}\boldsymbol{\beta}$, for any symmetric matrix $\mathbf{A}$.

Applying these:

Rearranging gives the normal equations:

This is a system of $k$ linear equations in $k$ unknowns. If $\mathbf{X}^T\mathbf{X}$ is invertible (which requires that no regressor is a perfect linear combination of the others), we can solve directly:

Interpreting the components. Each piece of this formula has a clear meaning:

The formula says: $\hat{\boldsymbol{\beta}}$ equals the covariation of regressors with outcomes, scaled by the inverse of the regressors’ own variation. This is the multivariate generalization of $\hat{\beta} = \text{Cov}(x,y) / \text{Var}(x)$.

Connection to projection (Section 5.4 revisited). The normal equations can be rewritten as:

This says the residual vector is orthogonal to every column of $\mathbf{X}$ — exactly the geometric projection condition from Section 5.4. The fitted values are the orthogonal projection of $\mathbf{y}$ onto the column space of $\mathbf{X}$:

where $\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$ is the hat matrix (or projection matrix). It “puts the hat on $\mathbf{y}$.” The hat matrix is symmetric ($\mathbf{H}^T = \mathbf{H}$) and idempotent ($\mathbf{H}^2 = \mathbf{H}$) — properties that projections always have.

5.8 From Coefficients to Standard Errors: Deriving $\text{Var}(\hat{\boldsymbol{\beta}})$¶

Computing $\hat{\boldsymbol{\beta}}$ is only half the job. To do inference (test hypotheses, build confidence intervals), we need to know how uncertain our estimates are. This requires deriving the variance of the estimator.

Step 1 — Express $\hat{\boldsymbol{\beta}}$ in terms of the true errors. Substituting $\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}$ into the OLS formula:

This is a fundamental decomposition: $\hat{\boldsymbol{\beta}}$ equals the truth $\boldsymbol{\beta}$ plus a noise term $(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\boldsymbol{\epsilon}$ that depends on the errors. Two immediate consequences:

1. Unbiasedness: If $E[\boldsymbol{\epsilon} \mid \mathbf{X}] = \mathbf{0}$, then $E[\hat{\boldsymbol{\beta}} \mid \mathbf{X}] = \boldsymbol{\beta}$.

2. Consistency: As $n \to \infty$, the noise term vanishes (under regularity conditions), so $\hat{\boldsymbol{\beta}} \xrightarrow{p} \boldsymbol{\beta}$.

Step 2 — The variance-covariance matrix: what it is and why we need it. For a scalar random variable $Z$, its variance is $\text{Var}(Z) = E[(Z - E[Z])^2]$ — a single number. But $\hat{\boldsymbol{\beta}}$ is a vector of $k$ random variables $(\hat{\beta}_0, \hat{\beta}_1, \ldots, \hat{\beta}_{k-1})^T$. We need a matrix that captures the variance of each component and the covariances between components. This is the variance-covariance matrix, defined as:

This is a $k \times k$ matrix. Its structure is:

• The diagonal entries $[\text{Var}(\hat{\boldsymbol{\beta}})]_{jj} = \text{Var}(\hat{\beta}_j)$ are the variances — the standard errors come directly from these: $\text{SE}(\hat{\beta}_j) = \sqrt{\text{Var}(\hat{\beta}_j)}$.

• The off-diagonal entries are covariances between different coefficient estimates — these matter when testing joint hypotheses (e.g., are multiple coefficients simultaneously zero?).

“Conditional on $\mathbf{X}$” means we treat the regressor matrix as fixed (non-random) and ask: over repeated samples of $\boldsymbol{\epsilon}$ drawn from the same error distribution, how much would $\hat{\boldsymbol{\beta}}$ vary? This is the standard frequentist perspective — the randomness comes only from $\boldsymbol{\epsilon}$, not from $\mathbf{X}$.

Step 3 — Derive the formula. From Step 1, we know $\hat{\boldsymbol{\beta}} - \boldsymbol{\beta} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\boldsymbol{\epsilon}$. Since $E[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta}$ (unbiasedness), we can substitute directly into the definition:

Using the transpose rule $(\mathbf{A}\mathbf{b})^T = \mathbf{b}^T\mathbf{A}^T$:

Since $\mathbf{X}$ is treated as fixed, only $\boldsymbol{\epsilon}\boldsymbol{\epsilon}^T$ is random, so the expectation passes through:

Here $\boldsymbol{\Omega} = E[\boldsymbol{\epsilon}\boldsymbol{\epsilon}^T]$ is the $n \times n$ variance-covariance matrix of the errors — its $(i,j)$ entry is $E[\epsilon_i \epsilon_j]$. The diagonal entries are the individual error variances $E[\epsilon_i^2]$, and the off-diagonals capture any correlation between errors at different observations.

This gives us the sandwich formula — so named for its bread-meat-bread structure: