Multivariate Regression - Proof regarding Constant

Question

A colleague of mine thinks that the constant in a multivariate regression is equal to the mean of the independent variable, usually denoted by $\bar{y}$. I disagree with my colleague, yet I somehow fail to show the opposite. I tried to prove my argument via the Frisch–Waugh–Lovell (FWL) theorem. However, I somehow get stuck and thought that the wisdom of the crowd might be able to help.

My data generating proces (DGP) is given by Equation 1: $$ \mathbf{y} = \beta_{0} \cdot \boldsymbol{\iota} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon} $$

with $\mathbf{y}$ being of dimension $[n \times 1]$, $\boldsymbol{\iota}$ being of dimension $[n \times 1]$, $\mathbf{X}$ being of dimension $[n \times k]$, $\boldsymbol{\beta}$ being of dimension $[k \times 1]$, and $\boldsymbol{\epsilon}$ being of dimension $[n \times 1]$.

By multiplying the equation above by the residual maker matrix $\mathbf{M}_{X} = \mathbf{I}_{N} - \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ from the left, I arrive at Equation 2: $$ \mathbf{M}_{X}\mathbf{y} = \beta_{0} \cdot \mathbf{M}_{X}\boldsymbol{\iota} + \mathbf{M}_{X}\boldsymbol{\epsilon} $$

The OLS estimator for $\beta_{0}$, i.e. $\hat{\beta}_{0}$, in Equation 2 is given by: $$ \hat{\beta}_{0} = (\boldsymbol{\iota}'\mathbf{M}_{X}'\mathbf{M}_{X}\boldsymbol{\iota})^{-1}\boldsymbol{\iota}'\mathbf{M}_{X}'\mathbf{M}_{X}\mathbf{y} $$

with $\mathbf{M}_{X}$ being symmetric and idempotent, we get Equation 3: $$ \hat{\beta}_{0} = (\boldsymbol{\iota}'\mathbf{M}_{X}\boldsymbol{\iota})^{-1}\boldsymbol{\iota}'\mathbf{M}_{X}\mathbf{y} $$

And this is where I am stuck. I tried to run a simulation in R, using the code below, however dividing the sum of the one residuals by the sum of the other residuals, yielded me some weird result, which I have not fully understood yet.

# define the number of individuals
n = 10000

# draw the epsilon shocks
eps <- rnorm(n, sd = 5)

# draw random Xs
covariate <- runif(n, min = 10, max = 20)

# define iota 
iota <- rep(1, n)

# define/compute y
y <- 5 + 3.25 * covariate + eps

# run the normal regression
summary(lm(y ~ covariate))

# compute the residuals using the FWL theorem
fwl_theorem_lm1 <- lm(iota ~ covariate)
fwl_theorem_lm2 <- lm(y ~ covariate)

# grab the residuals
res_lm1 <- fwl_theorem_lm1$residuals
res_lm2 <- fwl_theorem_lm2$residuals

sum(res_lm2) / sum(res_lm1)

Edit: I just realized that the regressions in the code should not include a constant. I adjusted the code accordingly, and now it the estimated value of the constant is returned from the division.

# define the number of individuals
n = 10000

# draw the epsilon shocks
eps <- rnorm(n, sd = 5)

# draw random Xs
covariate <- runif(n, min = 10, max = 20)

# define iota 
iota <- rep(1, n)

# define/compute y
y <- 5 + 3.25 * covariate + eps

# run the normal regression
summary(lm(y ~ covariate))

# compute the residuals using the FWL theorem
fwl_theorem_lm1 <- lm(iota ~ 0 + covariate)
fwl_theorem_lm2 <- lm(y ~ 0 + covariate)

# grab the residuals
res_lm1 <- fwl_theorem_lm1$residuals
res_lm2 <- fwl_theorem_lm2$residuals

sum(res_lm2) / sum(res_lm1)

Answer 1

You can disprove this using a counterexample in software.

set.seed(2022)
N <- 10
x <- runif(N)
y <- cbind(rnorm(N), rnorm(N))
L <- lm(y ~ x)
c(mean(y[,1]), mean(y[,2]))
summary(L)

I get that the intercept terms are $0.1459$ and $0.6199$, while the marginal means of the bivariate $y$ are $-0.281152500774767$ and $0.291819492757056$, which are not equal to the intercept terms from the regression.

If you want to do it by hand, feel free to calculate the OLS solution.