Why does $\mathbb{E}[s^2] \geq \sigma^2$ if we omit relevant variables (linear regression model)? Suppose that the "true" model is $Y = X\beta + Z \gamma + \varepsilon$ with the standard assumptions of the linear model. However, we only perform an OLS-regression on the variables contained in $X$.
I have slides for an econometrics course that claim that, under such an underspecification,
$$
\mathbb{E}[s^2] \geq \sigma^2,
$$
where $\sigma^2$ is the variance of the $\varepsilon_i$ and the sample variance is $s^2 := RSS/{(N-K)}$ (with $N$ as the number of observations and $K$ the number of regressors, including the intercept regressor).
Why is this the case?

To be completely clear we are assuming the structure specified above, along with the following: $\mathbb{E}[\varepsilon] = 0$, $\textrm{Var}(\varepsilon) = \sigma^2 I_N$, and $[X|Z] \in \mathbb{R}^{N \times (K + L)}$ has full rank and is deterministic.

My thoughts:
Suppose we define $\hat{\sigma}^2$ as the sample variance in a correctly specified model (i.e., regressing on $[X|Z]$). We know that $\mathbb{E}\hat{\sigma}^2 = \sigma^2$. So if it were generally true that $s^2 \geq \hat{\sigma}^2$, then the desired result would follow. However, I don't think that $s^2 \geq \hat{\sigma}^2$ is actually true, so I'm not sure where to go next.
 A: The correctly specified ("true") model is $Y = [X|Z]\begin{pmatrix} \beta \\ \gamma \end{pmatrix} + \varepsilon = X\beta + Z \gamma + \varepsilon$. However, our (underspecified) model is $Y = X\beta + \delta$, where in reality $\delta = Z\gamma + \varepsilon$.
Note that our estimator for $\sigma^2$ (given this model underspecification) is then $s^2 := \frac{\hat{\delta}^\intercal \hat{\delta}}{N-K}$, where the residual vector of our model is $\hat{\delta} := Y - X(X^\intercal X)^{-1}X^\intercal Y =: Y - HY = (I - H)Y =: MY$.
It will prove useful to rearrange $\hat{\delta}$:
$$
\hat{\delta} = (I-H)Y = (I-H)(X\beta + Z\gamma + \varepsilon) = MZ\gamma + M \varepsilon = M(Z\gamma + \varepsilon),
$$
since $HX= X$.
Now we examine the expected value of the residual sum of squares of our model:
$$
\mathbb{E}[\hat{\delta}^\intercal \hat{\delta}] = \mathbb{E}[(Z\gamma + \varepsilon)^\intercal M (Z\gamma + \varepsilon)] \quad \text{ (since $M^\intercal = M$ and $M^2 = M$)} \\
= \mathbb{E}[\gamma^\intercal Z^\intercal M Z \gamma] + 2\mathbb{E}[\gamma^\intercal Z^\intercal M \varepsilon] + \mathbb{E}[\varepsilon^\intercal M \varepsilon] \quad \text{(since $\mathbb{E}$ linear)} \\
= \underbrace{\gamma^\intercal Z^\intercal M Z \gamma}_{= \gamma^\intercal Z^\intercal M^\intercal M Z \gamma} + 2 \gamma^\intercal Z^\intercal M \underbrace{\mathbb{E}[\varepsilon]}_{=0} + \mathbb{E}[\varepsilon^\intercal M \varepsilon] \\
= \|M Z \gamma\|_2^2 + \mathbb{E}[\varepsilon^\intercal M \varepsilon]
$$
The last term on the right hand side above is $\sigma^2 (N-K)$, which we'll now show:
$$
\mathbb{E}[\varepsilon^\intercal M \varepsilon] = \mathbb{E}[\text{tr}(\varepsilon^\intercal M \varepsilon)] = \mathbb{E}[\text{tr}(M \varepsilon\varepsilon^\intercal)] \quad \text{(by properties of the trace function)} \\
= \text{tr}( \mathbb{E}[M \varepsilon\varepsilon^\intercal]) \quad \text{ (trace is a summation and $\mathbb{E}$ is linear)} \\
= \text{tr} \big( (I_N - H) \underbrace{\mathbb{E}[\varepsilon \varepsilon^\intercal]}_{= \text{Var}(\varepsilon) = \sigma^2 I_N} \big) = \sigma^2 \big( \underbrace{\text{tr}(I_N)}_{=N} - \underbrace{\text{tr}(H)}_{\substack{ = \text{tr}(X(X^\intercal X)^{-1}X^\intercal)} \\ =  \text{tr}(X^\intercal X (X^\intercal X)^{-1}) \\ = \text{tr}(I_K) = K} \big) = \sigma^2 (N-K).
$$
Hence, we have $\mathbb{E}[\hat{\delta}^\intercal \hat{\delta}] = \sigma^2(N-K) + \| MZ \gamma \|_2^2$, which yields our main result:
$$
\mathbb{E}[s^2] = \mathbb{E}\left[\frac{\hat{\delta}^\intercal \hat{\delta}}{N-K}\right] = \frac{1}{N-K} \mathbb{E}[\hat{\delta}^\intercal \hat{\delta}] = \sigma^2 + \frac{\|M Z\gamma\|_2^2}{N-K} \geq \sigma^2.
$$
