Software 5.6.1 (The residual-maker matrix). For the savings data with all four regressors (\(n = 50\), \(k = 5\) including the constant), build \(X\) and \(M\) directly and check Theorem 5.5.1:

> y <- LifeCycleSavings$sr
> X <- cbind(1, as.matrix(LifeCycleSavings[, c("pop15", "pop75",
+                                             "dpi", "ddpi")]))
> M <- diag(50) - X %*% solve(t(X) %*% X) %*% t(X)
> max(abs(M %*% M - M)); max(abs(M - t(M)))
[1] 3.663736e-15
[1] 7.979728e-16
> sum(diag(M))
[1] 45
> table(round(eigen(M, symmetric = TRUE)$values, 10))

 0  1
 5 45
> e <- M %*% y
> fit <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
> max(abs(e - resid(fit)))
[1] 4.831691e-13
> c(sum(e^2), deviance(fit))
[1] 650.713 650.713

Reading the output. \(M^2 - M\) and \(M - M'\) are zero up to rounding errors of order \(10^{-15}\), the limit of double-precision arithmetic. The trace is \(45 = n - k\), and the eigenvalues are \(45\) ones and \(5\) zeros, as Theorem 5.1.5 requires. \(MY\) reproduces the residuals that lm computes, and \(e'e = 650.71\) is the residual sum of squares that R calls the deviance.

One caution: forming \((X'X)^{-1}\) explicitly is fine for teaching but poor numerical practice. lm never does it, using a QR decomposition of \(X\) instead, which is more accurate when the regressors are nearly collinear. The \(50 \times 50\) matrix \(M\) is never formed in real computations either.

Stata: in Mata, M = I(50) - X*invsym(X’X)*X’ and trace(M). EViews: matrix M = @identity(50) - X*@inverse(@transpose(X)*X)*@transpose(X), then =@trace(M).

Software 5.6.2 (Simulating a quadratic form). Theorem 5.4.2(ii) checked by simulation, with a \(3 \times 3\) covariance matrix that is far from diagonal:

> set.seed(5)
> Sigma <- matrix(c(4, 2, 0, 2, 3, 1, 0, 1, 2), 3, 3)
> C <- t(chol(Sigma))                     # Sigma = C C'
> q <- replicate(10000, { x <- C %*% rnorm(3)
+                         drop(t(x) %*% solve(Sigma, x)) })
> c(mean = mean(q), var = var(q))
    mean      var
3.064371 6.386179
> round(quantile(q, c(0.5, 0.9, 0.95, 0.99)), 3)
   50%    90%    95%    99%
 2.398  6.364  8.036 11.766
> round(qchisq(c(0.5, 0.9, 0.95, 0.99), df = 3), 3)
[1]  2.366  6.251  7.815 11.345

Reading the output. The draws \(x = Cz\) with \(z \sim N(0, I)\) have covariance \(CC' = \Sigma \), by Theorem 5.2.1(i). The quadratic form has simulated mean \(3.06\) and variance \(6.39\), against \(3\) and \(6\) for \(\chi ^2_3\), and its quantiles are close to the \(\chi ^2_3\) quantiles. With 10 000 draws, differences of this size are simulation noise. solve(Sigma, x) computes \(\Sigma ^{-1}x\) without forming the inverse, the numerically better way.

Stata: drawnorm x1 x2 x3, cov(Sigma) generates the draws. EViews: nrnd with a Cholesky factor from @cholesky.