Deriving the OLS Estimators

Setting up the problem

You have \(n\) pairs \((x_i, y_i)\) and want to fit a line.

  • A line is just two numbers: an intercept \(b_0\) and a slope \(b_1\).
  • At point \(i\) the line predicts \(\hat{y}_i = b_0 + b_1 x_i\), while the real value is \(y_i\).
  • The gap between them is \(e_i = y_i - b_0 - b_1 x_i\).
  • Change \(b_0\) or \(b_1\) and every gap moves. Those two numbers are all you get to choose.

Residuals are not errors — keep them apart:

  • A residual \(e_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i\) is measured from the fitted line, so you can compute it.
  • An error \(u_i = y_i - \beta_0 - \beta_1 x_i\) is measured from the true, unknown \(\beta_0, \beta_1\), so you can never observe it.
  • The difference becomes essential the moment you reach unbiasedness and estimating the error variance.

What we minimize

You cannot minimize \(n\) gaps at once, so collapse them into one score. Not \(\sum_i e_i\): positive and negative misses cancel, and you can push it to \(-\infty\) by sliding the intercept. You want something that treats \(+5\) and \(-5\) as equally bad, so square them:

\[S(b_0, b_1) = \sum_{i=1}^{n} (y_i - b_0 - b_1 x_i)^2.\]

Absolute values would also do the job (that is least absolute deviations), but squares can be differentiated, which gives a clean formula, and the resulting estimator has good statistical properties.

First-order conditions

To find the minimum, set both partial derivatives to zero. The chain rule gives an outer derivative of \(2(\cdot)\) and an inner derivative of \(-1\) for \(b_0\) and \(-x_i\) for \(b_1\):

\[\frac{\partial S}{\partial b_0} = -2\sum_i (y_i - b_0 - b_1 x_i) = 0, \qquad \frac{\partial S}{\partial b_1} = -2\sum_i x_i (y_i - b_0 - b_1 x_i) = 0.\]

Normal equations

Divide each by \(-2\) and write the solutions as \(\hat{\beta}_0, \hat{\beta}_1\):

\[\sum_i (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i) = 0 \quad (1), \qquad \sum_i x_i (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i) = 0 \quad (2).\]

In residual form these are \(\sum_i e_i = 0\) and \(\sum_i x_i e_i = 0\). Both hold in every sample, whether or not the model is correct: they are a consequence of how we chose \(\hat{\beta}_0\) and \(\hat{\beta}_1\), not an assumption. Equation (1) holds because the model has an intercept (it is the derivative with respect to \(b_0\)). Fit a line through the origin instead and the guarantee \(\sum_i e_i = 0\) disappears. These two restrictions are also exactly why estimating the error variance later costs two degrees of freedom: \(\hat{\sigma}^2 = SSR/(n-2)\).

Why “normal” equations? Nothing to do with the normal distribution, and it doesn’t mean “ordinary.” Normal is the old geometry word for perpendicular (a line “normal” to a surface meets it at a right angle). Equation (2), \(\sum_i x_i e_i = 0\), says the residuals are perpendicular to \(x\): OLS fits the line by dropping a straight-down projection from the data onto it, and the leftover residual points at a right angle to the predictors. That right angle is exactly what makes \(SSR\) as small as it can be, since the shortest distance from a point to a line is the perpendicular one.

Solve (1) for the intercept

Expand and divide by \(n\), using \(\bar{y} = \tfrac{1}{n}\sum_i y_i\) and \(\bar{x} = \tfrac{1}{n}\sum_i x_i\):

\[\bar{y} - \hat{\beta}_0 - \hat{\beta}_1 \bar{x} = 0 \;\Longrightarrow\; \boxed{\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}.}\]

The fitted line passes through the point of means \((\bar{x}, \bar{y})\). This too depends on there being an intercept.

Substitute into (2) for the slope

Replace \(\hat{\beta}_0\) with \(\bar{y} - \hat{\beta}_1 \bar{x}\) and collect the \(\hat{\beta}_1\) terms:

\[\sum_i x_i\big[(y_i - \bar{y}) - \hat{\beta}_1 (x_i - \bar{x})\big] = 0 \;\Longrightarrow\; \hat{\beta}_1 = \frac{\sum_i x_i (y_i - \bar{y})}{\sum_i x_i (x_i - \bar{x})}.\]

Subtracting \(\bar{x}\) from the leading \(x_i\) changes nothing, because \(\sum_i \bar{x}(y_i - \bar{y}) = \bar{x}\sum_i (y_i - \bar{y}) = 0\) and likewise in the denominator. So

\[\boxed{\hat{\beta}_1 = \frac{\sum_i (x_i - \bar{x})(y_i - \bar{y})}{\sum_i (x_i - \bar{x})^2} = \frac{S_{xy}}{S_{xx}} = \frac{\widehat{\mathrm{Cov}}(x, y)}{\widehat{\mathrm{Var}}(x)}.}\]

We just divided, so the denominator must be non-zero. It equals \(\sum_i (x_i - \bar{x})^2\), which is zero only if every \(x_i\) is the same number. So \(x\) must vary across observations. That makes sense: if all your \(x\) values are identical, the data sit in a vertical strip and there is no slope to estimate.

Is it really a minimum?

Setting derivatives to zero finds a flat point, not necessarily the lowest one. Here it clearly is the lowest. \(S\) is a sum of squares, so it can never go below zero, and as a function of \(b_0\) and \(b_1\) it is a quadratic whose squared terms have positive coefficients (\(n\) and \(\sum_i x_i^2\)). It is a bowl that opens upward, so its one flat point is the bottom, provided \(x\) varies.

Find the minimum yourself

Drag the two sliders to move your line (dark) around the scatter. The red segments are your residuals, and the title reports your \(SSR = \sum_i e_i^2\). The green dashed line is the least-squares fit. No setting you pick can beat it: the right panel shows \(SSR\) as a bowl in the slope, and your line sits at its bottom only when \(b_1 = \hat{\beta}_1\).

#| standalone: true
#| viewerHeight: 560

library(shiny)

set.seed(12)
n  <- 40
x  <- round(runif(n, 1, 9), 2)
y  <- 2 + 1.3 * x + rnorm(n, sd = 2.2)
fit <- lm(y ~ x)
b0hat <- unname(coef(fit)[1]); b1hat <- unname(coef(fit)[2])
ssr_min <- sum(resid(fit)^2)

ui <- fluidPage(
  tags$head(tags$style(HTML("
    .eq-box { background:#f0f4f8; border-radius:6px; padding:14px;
              margin-top:14px; font-size:14px; line-height:1.8; }
    .eq-box b { color:#2c3e50; }
    .win  { color:#27ae60; font-weight:bold; }
    .miss { color:#e67e22; font-weight:bold; }
  "))),
  sidebarLayout(
    sidebarPanel(
      width = 3,
      sliderInput("b0", "Intercept  b0:", min = -4, max = 8, value = 0.0, step = 0.1),
      sliderInput("b1", "Slope  b1:",     min = -1, max = 3.5, value = 0.4, step = 0.05),
      uiOutput("box")
    ),
    mainPanel(
      width = 9,
      fluidRow(
        column(6, plotOutput("scatter", height = "380px")),
        column(6, plotOutput("bowl",    height = "380px"))
      )
    )
  )
)

server <- function(input, output, session) {

  ssr <- reactive(sum((y - input$b0 - input$b1 * x)^2))

  output$scatter <- renderPlot({
    par(mar = c(4.5, 4.5, 3, 1))
    plot(x, y, pch = 19, col = "#95a5a6",
         main = sprintf("Your line:  SSR = %.0f", ssr()),
         xlab = "x", ylab = "y")
    yhat <- input$b0 + input$b1 * x
    segments(x, y, x, yhat, col = "#e74c3c", lwd = 1)
    abline(a = b0hat, b = b1hat, col = "#27ae60", lwd = 2, lty = 2)
    abline(a = input$b0, b = input$b1, col = "#2c3e50", lwd = 2.5)
    legend("topleft", bty = "n", cex = 0.9,
           legend = c("Your line", "OLS (least squares)"),
           col = c("#2c3e50", "#27ae60"), lwd = c(2.5, 2), lty = c(1, 2))
  })

  output$bowl <- renderPlot({
    par(mar = c(4.5, 4.5, 3, 1))
    b1seq <- seq(-1, 3.5, length.out = 200)
    prof  <- sapply(b1seq, function(b) sum((y - input$b0 - b * x)^2))
    plot(b1seq, prof, type = "l", lwd = 2, col = "#2980b9",
         main = "SSR as you change the slope",
         xlab = "slope  b1", ylab = "SSR (at your intercept)")
    points(input$b1, ssr(), pch = 19, col = "#2c3e50", cex = 1.5)
    abline(v = b1hat, col = "#27ae60", lwd = 2, lty = 2)
    legend("top", bty = "n", cex = 0.85,
           legend = c("your slope", "OLS slope"),
           col = c("#2c3e50", "#27ae60"), pch = c(19, NA), lwd = c(NA, 2), lty = c(NA, 2))
  })

  output$box <- renderUI({
    gap <- ssr() - ssr_min
    cls <- if (gap < 1) "win" else "miss"
    tags$div(class = "eq-box",
      HTML(sprintf(
        "<b>Your SSR:</b> %.0f<br><b>Best possible:</b> %.0f<br><b>Gap above OLS:</b> <span class='%s'>%.0f</span><br><br><b>OLS estimates</b><br>&beta;&#770;<sub>0</sub> = %.2f<br>&beta;&#770;<sub>1</sub> = %.2f",
        ssr(), ssr_min, cls, gap, b0hat, b1hat)))
  })
}

shinyApp(ui, server)

A word on “MSE”

Two different things share this name. Don’t mix them up.

  • The regression-output kind (Root MSE). Take the residuals you just minimized, square them, and average with the right degrees of freedom: \(\text{MSE} = SSR/(n-2) = s^2\), your estimate of the error variance \(\sigma^2\). Its square root, Root MSE \(= s\), is the number Stata prints top-right of regress — the typical size of a residual, in the units of \(y\). (Where the \(n-2\) comes from is on the sampling distribution page.)
  • The estimator kind. For any estimate \(\hat{\theta}\) of a target \(\theta\), \(\text{MSE}(\hat{\theta}) = E[(\hat{\theta} - \theta)^2] = \mathrm{Var}(\hat{\theta}) + \text{Bias}(\hat{\theta})^2\). It measures how far the estimate lands from the truth, on average. If the estimator is unbiased, MSE is just its variance.

They measure different things: one is how far the data sit from the line (noise), the other how far an estimator sits from the truth (accuracy). But they connect, and the slope is where you see it:

  • Start with the estimator MSE of the slope: \(\text{MSE}(\hat{\beta}_1) = \mathrm{Var}(\hat{\beta}_1) + \text{Bias}(\hat{\beta}_1)^2\).
  • Under the OLS assumptions \(\hat{\beta}_1\) is unbiased, so its bias is zero and its MSE is just its variance, \(\mathrm{Var}(\hat{\beta}_1) = \sigma^2 / SST_x\).
  • That formula needs \(\sigma^2\), the error variance you can’t observe. So you plug in your estimate of it — the regression MSE, \(s^2 = \hat{\sigma}^2 = SSR/(n-2)\).
  • The result is the estimated variance \(s^2 / SST_x\), and its square root \(s / \sqrt{SST_x}\) is the slope’s standard error.

So the “how noisy is the data” MSE is exactly the ingredient that sets the “how uncertain is my slope” MSE: the scatter of the points around the line feeds straight through into the uncertainty of \(\hat{\beta}_1\).

Why the matrix version

Everything above is for one regressor. Real models have many: \(y = \beta_0 + \beta_1 x_1 + \dots + \beta_k x_k + u\). You could grind out \(k+1\) first-order conditions by hand, but the algebra explodes, and \(\hat{\beta}_1 = S_{xy}/S_{xx}\) no longer holds once the regressors are correlated with each other. Stacking the data into a matrix \(X\) solves all of it at once:

\[\hat{\beta} = (X'X)^{-1} X'y.\]

It is worth the extra notation for four reasons.

  • One formula, any number of X’s. With a single regressor you got \(\hat{\beta}_1 = S_{xy}/S_{xx}\). Add a second, third, or tenth regressor and you’d need a fresh hand-derivation each time. The matrix formula returns all the coefficients in one shot, and each one is already computed holding the others fixed — that is the “controlling for everything else” (ceteris paribus) reading you actually want in a real regression.
  • Every standard error comes out together. To run a \(t\)-test or build a confidence interval you need the standard error of each coefficient. The matrix \(\sigma^2 (X'X)^{-1}\) is a little table whose diagonal lists the variance of every coefficient at once. Take square roots and you have all the standard errors, from one object, with no extra work.
  • It reveals when two X’s are “too alike” (multicollinearity). If two regressors carry nearly the same information (say height in inches and height in centimeters), the data can’t separate their individual effects. In matrix terms \(X'X\) becomes almost impossible to invert, and the standard errors blow up, so the coefficients look wildly uncertain. The matrix form makes this problem visible instead of hiding it.
  • It is what the software actually runs. When you type regress in Stata or lm in R, it is not grinding through by-hand algebra — it is solving these normal equations in matrix form, \(\hat{\beta} = (X'X)^{-1}X'y\). So the matrix version is not optional theory; it is what happens under the hood every time you fit a regression.

So: the scalar derivation is where you understand least squares; the matrix form is the engine you use once there is more than one \(x\). The full matrix treatment is on The Algebra Behind OLS.

Where the standard errors come from (matrix form)

In matrix form, all the variances and covariances of the estimates live in one object, the variance–covariance matrix:

\[\widehat{\mathrm{Var}}(\hat{\beta}) = \hat{\sigma}^2 (X'X)^{-1}.\]

The standard error of coefficient \(j\) is \(\hat{\sigma}\) times the square root of the \(j\)-th diagonal entry of \((X'X)^{-1}\):

\[\mathrm{SE}(\hat{\beta}_j) = \hat{\sigma}\,\sqrt{\big[(X'X)^{-1}\big]_{jj}}.\]

This is where “blow up” becomes concrete. Inverting a matrix divides by its determinant, and when two regressors are nearly collinear the columns of \(X\) almost line up, so \(\det(X'X)\) sits close to zero. Dividing by a near-zero number makes the entries of \((X'X)^{-1}\) huge — the diagonal ones included — so the variances, and with them the standard errors, explode. “\(X'X\) is impossible to invert” and “the standard errors blow up” are the same sentence.

What makes a standard error big or small

Once there is more than one regressor, the standard error of coefficient \(j\) is

\[\mathrm{SE}(\hat{\beta}_j) = \frac{\sigma}{\sqrt{SST_j \, (1 - R_j^2)}},\]

with three moving parts:

  • \(\sigma\) — the noise around the line (the error SD).
  • \(SST_j = \sum_i (x_{ji} - \bar{x}_j)^2\) — how much regressor \(j\) varies.
  • \(R_j^2\) — the \(R^2\) from regressing \(x_j\) on all the other regressors: how much of \(x_j\) the others already explain.

Reading off the levers:

  • The SE goes up when the noise \(\sigma\) is larger, when \(x_j\) barely varies (small \(SST_j\) — too few points or too narrow a range of \(x_j\)), or when \(x_j\) is nearly explained by the other regressors (\(R_j^2 \to 1\)). That last case is multicollinearity: the factor \(1/(1 - R_j^2)\) is the variance inflation factor, and it runs off to infinity as \(R_j^2 \to 1\) — the same thing as \(X'X\) becoming impossible to invert.
  • The SE goes down with less noise, more spread in \(x_j\) (or simply more data, which grows \(SST_j\)), and regressors that are uncorrelated with each other (\(R_j^2 \to 0\)).

With exactly two regressors, \(R_j^2\) is just \(\rho^2\), the square of their correlation, so the inflation factor is \(1/(1-\rho^2)\) — the multiplier you can watch in the demo below.

See it happen. Drag the correlation between the two regressors toward 1. The left panel shows the coefficient estimates across many repeated samples: as the X’s line up, the cloud stretches into a diagonal — the data can pin down the two effects’ combination but not each one on its own — and both standard errors blow up. The right panel is the SE multiplier \(1/\sqrt{1-\rho^2}\), which runs off to infinity as \(\rho \to 1\) (that is \(X'X\) becoming impossible to invert).

#| standalone: true
#| viewerHeight: 520

library(shiny)

ui <- fluidPage(
  tags$head(tags$style(HTML("
    .eq-box { background:#f0f4f8; border-radius:6px; padding:14px;
              margin-top:14px; font-size:14px; line-height:1.8; }
    .eq-box b { color:#2c3e50; }
    .hot { color:#e74c3c; font-weight:bold; }
  "))),
  sidebarLayout(
    sidebarPanel(
      width = 3,
      sliderInput("rho", "Correlation between X1 and X2:",
                  min = 0, max = 0.99, value = 0.5, step = 0.03),
      helpText("Push it toward 1 and the two regressors become nearly the same variable."),
      uiOutput("box")
    ),
    mainPanel(
      width = 9,
      fluidRow(
        column(6, plotOutput("cloud",   height = "380px")),
        column(6, plotOutput("securve", height = "380px"))
      )
    )
  )
)

server <- function(input, output, session) {

  sims <- reactive({
    rho <- input$rho; n <- 100; ns <- 400
    b1 <- numeric(ns); b2 <- numeric(ns)
    for (i in seq_len(ns)) {
      z1 <- rnorm(n); z2 <- rnorm(n)
      x1 <- z1
      x2 <- rho * z1 + sqrt(1 - rho^2) * z2
      y  <- 1 + 2 * x1 + 1 * x2 + rnorm(n, sd = 1)
      cf <- coef(lm(y ~ x1 + x2))
      b1[i] <- cf["x1"]; b2[i] <- cf["x2"]
    }
    list(b1 = b1, b2 = b2)
  })

  output$cloud <- renderPlot({
    s <- sims(); par(mar = c(4.5, 4.5, 3, 1))
    plot(s$b1, s$b2, pch = 19, col = adjustcolor("#2980b9", 0.30),
         xlim = c(-3, 7), ylim = c(-4, 6),
         main = "Estimates across repeated samples",
         xlab = "b1-hat", ylab = "b2-hat")
    points(2, 1, pch = 4, col = "#e74c3c", cex = 2.2, lwd = 3)
    legend("topright", bty = "n", legend = "truth (2, 1)",
           pch = 4, col = "#e74c3c", pt.cex = 1.6, pt.lwd = 2)
  })

  output$securve <- renderPlot({
    par(mar = c(4.5, 4.5, 3, 1))
    rseq <- seq(0, 0.99, length.out = 120)
    plot(rseq, 1 / sqrt(1 - rseq^2), type = "l", lwd = 2, col = "#8e44ad",
         main = "How much the standard error inflates",
         xlab = "correlation between X1 and X2", ylab = "SE multiplier")
    abline(v = input$rho, col = "#e74c3c", lwd = 2, lty = 2)
    points(input$rho, 1 / sqrt(1 - input$rho^2), pch = 19, col = "#e74c3c", cex = 1.4)
  })

  output$box <- renderUI({
    s <- sims()
    se1 <- sd(s$b1); se2 <- sd(s$b2)
    vif <- 1 / (1 - input$rho^2)
    cls <- if (input$rho > 0.9) "hot" else ""
    tags$div(class = "eq-box",
      HTML(sprintf(
        "<b>Spread of b1-hat:</b> <span class='%s'>%.2f</span><br><b>Spread of b2-hat:</b> <span class='%s'>%.2f</span><br><br><b>SE inflation</b><br>1/(1-&rho;<sup>2</sup>) = %.1f&times;",
        cls, se1, cls, se2, vif)))
  })
}

shinyApp(ui, server)

What would ideal data look like? Read it straight off the formula \(\mathrm{SE}(\hat{\beta}_j) = \sigma / \sqrt{SST_j(1 - R_j^2)}\): you want a large sample, each \(X\) spread out (large \(SST_j\)), the regressors uncorrelated with one another (every \(R_j^2 \approx 0\), so no inflation), and low noise \(\sigma\). That is essentially a well-designed experiment — plenty of independent variation in each \(X\) and no redundancy between them. Then every variance inflation factor equals 1, and each coefficient is pinned down as tightly as the data allow.

Connections

  • Sampling Distribution of OLS — once you have \(\hat{\beta}_1 = S_{xy}/S_{xx}\), treat it as a random variable: its mean, its variance \(\sigma^2/SST_x\), and where the \(n-2\) comes from.
  • The Algebra Behind OLS — the same estimator in matrix form, \(\hat{\beta} = (X'X)^{-1}X'y\).
  • Residuals & Controls — the residuals \(e_i\) this page produces, and how they behave.