The setting of my problem is as follows:
Consider the general linear model $Y = X\boldsymbol{\beta} + \epsilon$, where $\beta = \begin{bmatrix} \beta_0 \\ \boldsymbol{\beta_1} \end{bmatrix}$ and $X = \begin{bmatrix} \mathbf{1} & \mathbf{X_1} \end{bmatrix}$, $\epsilon_i$ are independent with mean zero and variance $\sigma^2$. Consider the case where at least one element of $\boldsymbol{\beta_1}$ is not equal to zero. Prove that $\mathbb{E}[MSR]>\sigma^2$.
What I did so far:
$SSR = \sum(\hat{Y_i}-\bar{Y})^2 = \hat{Y}^T\hat{Y}-n\bar{Y}^2=Y^THY-n\bar{Y}^2$ where $H = X(X^TX)^{-1}X^T$.
Using the identity $\mathbb{E}[Y^THY]=\mathbb{E}[Y]^TH\mathbb{E}[Y]+tr[H*Var(Y)]$ I was able to get:
$\mathbb{E}[SSR] = \beta^TX^TX\beta+p\sigma^2-n\bar{Y}^2$
This is where I am stuck. It would suffice to show $\beta^TX^TX\beta-n\bar{Y}^2>0$, for then $MSR = \frac{SSR}{p-1} > \frac{p\sigma^2}{p-1} > \sigma^2$, but I cannot seem to figure out how to do that.
Your help is much appreciated!
$\mathbf{EDIT}$
My initial attempt was incorrect: I moved $\bar{Y}^2$ past the expectation, which is obviously wrong. Now, I attempted the problem again and I think it may have worked:
$$SSR = \sum(\hat{Y_i}-\bar{Y})^2 = \hat{Y}^T\hat{Y}-n\bar{Y}^2=Y^THY-n\bar{Y}^2$$
$\mathbb{E}[SSR] = \mathbb{E}[Y^THY-n\bar{Y}^2] = \mathbb{E}[Y^THY] -n\mathbb{E}[\bar{Y}^2]$
$=\mathbb{E}[Y]^TH\mathbb{E}[Y]+tr[H*Var(Y)] - n\mathbb{E}[\bar{Y}^2]$
$=\boldsymbol{\beta}^TX^T(X(X^TX)^{-1}X^T)X\boldsymbol{\beta} + p\sigma^2- n\mathbb{E}[\bar{Y}^2]$ (using $tr(H) = p$)
$=\boldsymbol{\beta}^TX^TX\boldsymbol{\beta}+p\sigma^2- n\mathbb{E}[\bar{Y}^2] = \mathbb{E}[Y]^T\mathbb{E}[Y]+p\sigma^2- n\mathbb{E}[\bar{Y}^2]$
$=\sum{\mathbb{E}[Y_i]^2} - \frac{1}{n}\mathbb{E}[(\sum{Y_i})^2]+p\sigma^2$
$=\sum{\mathbb{E}[Y_i]^2} - \frac{1}{n}\big(\mathbb{E}[(\sum{Y_i})^2]-\mathbb{E}[\sum{Y_i}]^2\big) - \frac{1}{n}\mathbb{E}[\sum{Y_i}]^2+p\sigma^2$
$=\sum{\mathbb{E}[Y_i]^2} - \frac{1}{n}Var(\sum{Y_i}) - \frac{1}{n}\big(\sum{\mathbb{E}[Y_i]}\big)^2+p\sigma^2$
Since the $Y_i$ are independent:
$=\sum{\mathbb{E}[Y_i]^2} - \frac{1}{n}\big(\sum{\mathbb{E}[Y_i]}\big)^2+(p-1)\sigma^2$
Now, we can simply use the fact that $\sum{x_i^2} \ge n\bar{x}^2$ with equality iff $\sum{(x_i-\bar{x})^2}=0$, by taking $x_i=\mathbb{E}[Y_i]$ (for a proof see How can I prove that the mean of squared data points is greater than the square of the mean of the data points?):
$$\mathbb{E}[SSR]\ge(p-1)\sigma^2$$ with equality iff $\mathbb{E}[Y_i]=const$.
Assuming uniqueness of the OLS estimates, if at least one element of $\boldsymbol\beta_1$ is non-zero this gives:
$$\mathbb{E}[MSR] = \frac{1}{p-1}\mathbb{E}[SSR] > \sigma^2$$
For a single observation, write $Y-\bar Y=\beta_0+\beta_1 X_1+...+\beta_pX_p+\epsilon-(\beta_0+\beta_1 \bar X_1+...+\beta_p\bar X_p+\bar \epsilon)=\beta_1(X_1-\bar X_1)+...+\beta_p(X_p-\bar X_p)+(\epsilon-\bar \epsilon)$. Then the square of this for a single observation is $(Y-\bar Y)^2=\beta_1^2(X_1-\bar X_1)^2+...+\beta_p^2(X_p-\bar X_p)^2+(\epsilon-\bar \epsilon)^2+2\beta_1\beta_p(X_1-\bar X_1)+2\beta_1(X_1-\bar X_1)(\epsilon-\bar \epsilon)+2\beta_p(X_p-\bar X_p)(\epsilon-\bar \epsilon)+...$
We need two more pieces. First, the expected value of $E[(\epsilon-\bar \epsilon)^2]=E(\epsilon^2-2\epsilon\bar \epsilon+\bar \epsilon^2)=\sigma^2+\frac{\sigma^2}{n}-\frac{2\sigma^2}{n}=\sigma^2\left(\frac{n-1}{n}\right)$. Also, $E(\epsilon-\bar\epsilon)=0$. Taking the sum over all observations, we have that for the sample, $\sum(Y_i-\bar Y)^2=\beta_1^2\sum(X_{i1}-\bar X_1)^2+...+\beta_p^2(X_{ip}-\bar X_p)^2+(\epsilon_i-\bar \epsilon)^2+2\beta_1\beta_p(X_{i1}-\bar X_1)+2\beta_1(X_{i1}-\bar X_1)(\epsilon_i-\bar \epsilon)+2\beta_p(X_{ip}-\bar X_p)(\epsilon_i-\bar \epsilon)+...$
Then the expected value of this is $\begin{split}E(SSR)&=\beta_1^2\sum E(X_{i1}-\bar X_1)^2+...+\beta_p^2\sum E(X_{ip}-\bar X_p)^2+\sum \underbrace{E(\epsilon_i-\bar \epsilon)^2}_{=\sigma^2\left(\frac{n-1}{n}\right)}+2\beta_1\beta_p\sum E(X_{i1}-\bar X_1)+\underbrace{2\beta_1\sum E(X_{i1}-\bar X_1)(\epsilon_i-\bar \epsilon)+2\beta_p\sum E(X_{ip}-\bar X_p)(\epsilon_i-\bar \epsilon)}_{=0}+...\\ &=\beta_1^2\sum E(X_{i1}-\bar X_1)^2+...+\beta_p^2\sum E(X_{ip}-\bar X_p)^2+2\beta_1\beta_p\sum E(X_{i1}-\bar X_1)+...+\sigma^2(n-1)\end{split}$
Now, if exactly one of the $\beta_j$'s was not equal to 0 the terms other than the $\sigma^2(n-1)$ would just be $\beta_i^2\sum E(X_{ij}-\bar X_j)^2>0$. But here you somehow need to make an argument that those terms collectively are positive.