Mastering the OLS Method in Multiple Regression Estimation

If you’ve ever tried to find a pattern in the real world, you’ve engaged in a form of regression. We’re constantly trying to figure out how different factors influence an outcome. For instance, what *really* drives the price of a house? It’s not just its size. It’s also the location, the number of bedrooms, the age of the property, and maybe even the quality of the local school district. This is a classic multiple regression problem, where we have multiple ‘X’ variables (predictors) trying to explain a single ‘Y’ variable (the outcome). But with all these factors, how do we find the “best” formula? How do we assign the right weight to size versus location? This is where the workhorse of econometrics comes in: the Ordinary Least Squares (OLS) method. It’s a powerful and elegant technique for finding the best-fitting line through a cloud of data points. In this post, we’re going to pop the hood and look at the engine of OLS, using a bit of matrix algebra to see exactly how it works, why we trust it, and how we measure its precision.

Deriving the OLS estimator with matrix algebra 🔗

When you first learned about regression, you probably started with a simple model: $Y = \beta_0 + \beta_1X + \epsilon$. You might have even learned some messy-looking formulas for $\beta_0$ and $\beta_1$ full of summation signs. This works fine for one, maybe two, explanatory variables. But what about our housing example with four, five, or even twenty predictors? The math becomes a nightmare of subscripts and sums.

This is where matrix algebra comes to the rescue. It’s a powerful mathematical shorthand that lets us express this complex system in one clean, simple equation:

$Y = X\beta + \epsilon$

Let’s quickly break down what’s in these matrices. If we have ‘n’ observations (e.g., $n$ houses) and ‘k’ parameters (e.g., $k-1$ predictors plus one intercept):

• $Y$ is an $n \times 1$ vector (a single column) containing all our observed outcomes (e.g., the prices of all $n$ houses).
• $X$ is our $n \times k$ matrix of data. Each row is one observation (one house), and each column is one predictor variable (e.g., square footage, number of bedrooms). We also include a column of ones to represent the intercept.
• $\beta$ (beta) is the $k \times 1$ vector of coefficients we *want to find*. This is the “answer” we’re looking for (e.g., the specific price effect of one extra bedroom).
• $\epsilon$ (epsilon) is the $n \times 1$ vector of errors, or “residuals”-the parts of $Y$ that our model *can’t* explain.

The goal of OLS is to find the specific vector of coefficients (which we’ll call $\hat{\beta}$, or “beta-hat”) that makes the errors as small as possible. But what does “small as possible” mean? We can’t just make the simple sum of errors zero, because large positive errors and large negative errors would cancel each other out, telling us nothing. Instead, OLS seeks to minimize the Residual Sum of Squares (RSS). We square each error (making them all positive) and then add them all up. Our goal is to find the $\hat{\beta}$ that results in the *lowest possible RSS*.

In matrix terms, the vector of residuals is $e = Y – X\hat{\beta}$. The RSS is simply $e’e$ (e-transpose times e), which is a shorthand way of saying “square each element in the vector $e$ and add them all up.”

So, our objective is to minimize: $RSS = e’e = (Y – X\hat{\beta})'(Y – X\hat{\beta})$

How do we find the minimum of a function? Calculus! We take the derivative of the RSS with respect to our $\hat{\beta}$ vector, set that derivative to zero (finding the point where the “slope” is flat, i.e., the bottom of the valley), and solve for $\hat{\beta}$. We’ll spare you the full matrix calculus derivation, but this process of differentiation leads us to a famous intermediate step called the Normal Equations:

$(X’X)\hat{\beta} = X’Y$

This equation is the heart of OLS. Look at it closely. $(X’X)$ is just a $k \times k$ matrix we can calculate from our data. $X’Y$ is a $k \times 1$ vector we can also calculate. The only unknown is $\hat{\beta}$. This is now just a simple algebra problem, like $5x = 10$. How do you solve for $x$? You multiply by the inverse of 5 (i.e., $1/5$).

We do the exact same thing here. To isolate $\hat{\beta}$, we pre-multiply both sides by the inverse of the $(X’X)$ matrix, which we write as $(X’X)^{-1}$.

$(X’X)^{-1}(X’X)\hat{\beta} = (X’X)^{-1}X’Y$

A matrix multiplied by its own inverse becomes the Identity matrix (the matrix equivalent of the number 1), so it just “disappears.” This leaves us with the single most important formula in multiple regression:

$\hat{\beta} = (X’X)^{-1}X’Y$

This is the OLS estimator. It’s a complete “recipe” for finding the best-fitting coefficients. You just plug in your matrix of $X$ data and your vector of $Y$ data, perform these matrix operations, and the result is the vector of $\hat{\beta}$ coefficients that best explains $Y$.

Why we trust OLS: The mean and unbiasedness of estimators 🔗

That’s a very impressive-looking formula. But how do we know it’s *any good*? Why should we trust the coefficients it gives us? This brings us to the properties of the OLS estimator, and the most fundamental one is unbiasedness.

What does it mean for an estimator to be unbiased? Let’s use an analogy. Imagine an archer shooting at a target. [Image: An archery target with arrows clustered around the bullseye] An unbiased archer is one whose arrows, on average, hit the dead-center bullseye. Any *single* shot might be slightly high, low, left, or right. But if you averaged out 1,000 of their shots, the average position would be the exact center. A *biased* archer, in contrast, might be very precise (all their arrows land close together) but consistently hits high and to the left. Their average shot is not on the bullseye.

In econometrics, the “bullseye” is the *true, unknown population parameter*, $\beta$. Our “shot” is the $\hat{\beta}$ we calculate from our one sample of data. Because our sample is just one random draw from the population, our $\hat{\beta}$ will almost certainly be slightly different from the true $\beta$. Unbiasedness means that if we could (hypothetically) draw thousands of different samples from the population and calculate $\hat{\beta}$ for each one, the *average* of all those $\hat{\beta}$’s would be equal to the true $\beta$.

Formally, we want to prove that the Expected Value of our estimator, $E[\hat{\beta}]$, is equal to the true $\beta$. Let’s see if it is. The proof is surprisingly elegant.

1. We start with our estimator: $\hat{\beta} = (X’X)^{-1}X’Y$
2. We don’t know $Y$, but we know the *true model* for $Y$: $Y = X\beta + \epsilon$. Let’s substitute that in.
3. $\hat{\beta} = (X’X)^{-1}X'(X\beta + \epsilon)$
4. Now, let’s distribute the $(X’X)^{-1}X’$ term across the parenthesis.
5. $\hat{\beta} = (X’X)^{-1}X’X\beta + (X’X)^{-1}X’\epsilon$
6. Look at that first term. We have $(X’X)^{-1}$ right next to $(X’X)$. As we said before, a matrix times its inverse is just the Identity matrix $I$. So the first term simplifies to $I\beta$, which is just $\beta$.
7. $\hat{\beta} = \beta + (X’X)^{-1}X’\epsilon$

This equation is itself very insightful. It says our *estimated* coefficient ($\hat{\beta}$) is equal to the *true* coefficient ($\beta$) plus some “error term” that depends on our $X$ data and the unobservable errors $\epsilon$.

Now, let’s take the expected value (the “average”) of both sides to check for bias:

$E[\hat{\beta}] = E[\beta + (X’X)^{-1}X’\epsilon]$

We can split the expectation: $E[\hat{\beta}] = E[\beta] + E[(X’X)^{-1}X’\epsilon]$

The true $\beta$ is a fixed, unknown constant, not a random variable. So, the expected value of $\beta$ is just $\beta$. And in OLS, we treat our $X$ data as fixed (non-stochastic). The only random part left is the error term $\epsilon$. So, we can pull the fixed parts out of the expectation:

$E[\hat{\beta}] = \beta + (X’X)^{-1}X’ \cdot E[\epsilon]$

What is $E[\epsilon]$, the expected value of the error term? A fundamental assumption of OLS is that the error term has a mean of zero. This means our model is correctly specified, and the errors are just random “noise,” not a systematic, predictable bias. If $E[\epsilon] = 0$, then that entire second term becomes $(X’X)^{-1}X’ \cdot 0$, which is just a big vector of zeros.

And so, we are left with:

$E[\hat{\beta}] = \beta$

This proves it! The OLS estimator is unbiased. On average, it hits the bullseye. This is a crucial property and a primary reason why it’s so trusted. It’s the “U” in BLUE-the “Best Linear Unbiased Estimator,” a title given to OLS by the famous Gauss-Markov Theorem.

How precise are our estimates? The variance of OLS estimators 🔗

Unbiasedness is great. We know our archer is aiming at the bullseye. But that’s not the whole story. Is the archer a “precise” shooter? Are all their (unbiased) shots clustered tightly around the center, or are they spread all over the target? We want an estimator with *low variance*-a precise archer.

We need a way to measure this “spread” for our $\hat{\beta}$ coefficients. This is the variance-covariance matrix of $\hat{\beta}$. It’s a $k \times k$ matrix that tells us two things:

• The diagonal elements: These are the variances of each coefficient estimate (e.g., $\text{Var}(\hat{\beta}_1)$, $\text{Var}(\hat{\beta}_2)$). This is a measure of precision for each *individual* variable. A small variance means a precise estimate.
• The off-diagonal elements: These are the *covariances* between estimates (e.g., $\text{Cov}(\hat{\beta}_1, \hat{\beta}_2)$). They tell us how our estimates “move together.” For example, if we overestimate the effect of square footage, do we also tend to overestimate the effect of location?

To find this matrix, we start from our equation $\hat{\beta} = \beta + (X’X)^{-1}X’\epsilon$. The variance of $\hat{\beta}$ is the variance of this whole expression. Since $\beta$ is a constant, it has no variance, so we just need the variance of the error part: $\text{Var}(\hat{\beta}) = \text{Var}((X’X)^{-1}X’\epsilon)$.

We’ll skip the matrix algebra rules, but the derivation relies on a key assumption: homoscedasticity and no autocorrelation. This means we assume that all error terms $\epsilon_i$ have the same variance (call it $\sigma^2$, sigma-squared) and are not correlated with each other. In matrix form, this is written as $\text{Var}(\epsilon) = \sigma^2I$, where $I$ is the identity matrix.

When you plug this assumption into the variance formula and simplify the matrix algebra, you get another beautifully clean result:

$\text{Var}(\hat{\beta}) = \sigma^2(X’X)^{-1}$

This formula is incredibly important. It tells us exactly what determines the precision of our OLS estimates.

1. $\sigma^2$ (The Error Variance): This is the “noise” in the model, the inherent randomness in $Y$ that cannot be explained by $X$. If this number is large (the data points are very spread out far from the true regression line), then $\text{Var}(\hat{\beta})$ will also be large. This makes perfect sense: in a “noisy” dataset, it’s harder to be precise about our estimates.
2. $(X’X)^{-1}$ (The Data Structure): This part is all about our $X$ variables. It generally tells us that the more “spread out” (more variation) our $X$ variables are, the *smaller* this matrix becomes, and thus the *more precise* our estimates will be. It also contains the impact of multicollinearity. If two of your $X$ variables are very highly correlated (e.g., you include both “house size in sq. ft.” and “house size in sq. meters”), the $(X’X)$ matrix becomes “ill-conditioned” or “singular,” and its inverse $(X’X)^{-1}$ “blows up,” leading to massive, unstable variances.

Finding the missing piece: Estimating the error variance (σ²) 🔗

We have a problem. Our wonderful variance formula, $\sigma^2(X’X)^{-1}$, contains $\sigma^2$-the *true population variance of the errors*. We can never, ever know this true value. All we have is our sample.

So, we must *estimate* it. Our best guess for the unobservable errors ($\epsilon$) is our vector of *residuals* ($e$). Remember, $e = Y – X\hat{\beta}$. These are the “leftovers” from our regression, the differences between the actual $Y$ values and our model’s predicted $\hat{Y}$ values. [Image: A scatter plot with a regression line, showing the vertical residual lines from each point to the line] We can use these residuals to calculate our estimate for the error variance, which we’ll call $\hat{\sigma}^2$ (“sigma-hat-squared”).

A simple first guess might be to just take the average of the squared residuals: $\frac{RSS}{n}$, where $RSS = e’e$. But it turns out this estimate is *biased*. It will, on average, be slightly too small. Why? Because we *used* our $n$ data points to estimate our $k$ coefficients. Our fitted line is *by definition* as close as possible to those $n$ points. We’ve “used up” some of the randomness in the data just by fitting the line. We’ve lost some degrees of freedom.

To get an *unbiased* estimator for the error variance, we must divide the RSS not by $n$, but by the *degrees of freedom* remaining: $n – k$. (Where $k$ is the number of parameters we estimated, including the intercept).

This gives us the final, crucial formula for our estimated error variance, also known as the Mean Square Error (MSE):

$\hat{\sigma}^2 = \frac{RSS}{n-k}$

This $\hat{\sigma}^2$ is our unbiased, data-driven estimate for the “noise” in the model. Now we can finally compute a usable variance-covariance matrix by plugging this estimate back into our variance formula:

$\widehat{\text{Var}(\hat{\beta})} = \hat{\sigma}^2(X’X)^{-1} = \left(\frac{RSS}{n-k}\right)(X’X)^{-1}$

This is it! This is the matrix your statistical software (like R, Stata, or Python) calculates. The numbers you see in your regression output table under “Standard Error” are simply the square roots of the diagonal elements of this very matrix. These standard errors are what allow us to build confidence intervals and perform t-tests, giving us a complete picture of not only our “best guess” ($\hat{\beta}$) but also our “level of precision” (its standard error).

What do you think? Having seen the matrix formula $\hat{\beta} = (X’X)^{-1}X’Y$, does it make more sense why high correlation between your $X$ variables (multicollinearity) is such a big problem for OLS? (Hint: It makes the $(X’X)$ matrix difficult or impossible to invert.) And now that you know standard errors depend on $n-k$, what do you think happens to the precision of your estimates if you add 20 new variables to a model but only have 30 data points?

References

1. https://online.stat.psu.edu/stat501/lesson/9/9.3
2. https://nptel.ac.in/courses/110104075
3. https://stats.oarc.ucla.edu/other/annotated-output/regression-analysis/