What is wrong with the linear regression variance estimator when p+1>n? According to linear regression formula, assuming the predictors are nonrandom, $${\rm var}[\hat{\beta}]=MSE *(X'X)^{-1}$$ where MSE is an unbiased estimator of the population variance.
Intuitively, when the number of predictors (p + 1) is greater than the number of data, the regression model is overfitting the data, variance of $\hat{\beta}$ should be big.
However, according to the formula, MSE = 0 and we have $${\rm var}[\hat{\beta}]=MSE *(X'X)^{-1} = 0$$
What goes wrong? There must be some implicit assumption with this formula that failed in this scenario.
 A: Nothing is wrong with the formula
No, it is not correct that the estimated variance of $\hat\beta$ is zero. There are a couple of errors in your reasoning.
First, the MSE is NA rather than 0. The MSE is equal to the residual sum of squares divided by the residual degrees of freedom (df).
Numerator and denominator are both zero in this case, so the MSE is undefined.
This makes intuitive sense because there is no way to estimate the true variance $\sigma^2$ here.
Second, the elements of $(X'X)^{-1}$ are infinitely large because $X'X$ is singular. In fact, $p+1-n$ of the eigenvalues of $X'X$ are exactly zero. The existence of zero eigenvalues means that $X'X$ is not invertible, in the sense that it doesn't have a finite inverse.
The fact that $p+1 > n$ also means that $\hat\beta$ is not uniquely determined. There are a range of $\hat\beta$ values that all yield a perfect fit to the data, and the data yields no information whatsoever to as to which one of these values is closest to the true value.
So it makes intuitive sense to think of the variance estimate as being infinite.
In summary, while the variance formula doesn't give a value in this situation, it does guide you as to what is happening as well as any formula can.
So how did you compute the wrong value?
The question still remains though, how did you compute the variance as zero, when no computation is actually possible?
If I had to guess what you've done wrong, it might that you have tried to compute the residual df naively using the textbook formula $n-p-1$, with $p$ equal to the number of predictor variables, excluding the intercept. That gives a negative number in this case, which of course is nonsense because you can't have less than no degrees of freedom. That would indeed give you the incorrect result of MSE$=0/(n-p-1)=0$ that you state in your question.
The correct residual df is actually defined as $n-p'$ where $p'$ is the rank of $X$. Now $p'$ is indeed equal to $p+1$ in the classic textbook context that the number of predictors is $p<n$ and all the predictors are linearly independent of each other and of the intercept. More generally however, $p'$ can never be greater than $n$, because the rank of a matrix cannot be larger than the minimum of the number of rows and the number of columns.
In your case you will have $p'=n$, assuming that there aren't any unnecessary linear dependencies between your predictors.
Of course, it would be easier to identify your error with less guessing if you did explain how you calculated MSE.
I also wonder how you obtained a value for $(X'X)^{-1}$, considering that the inverse doesn't exist and the computation is impossible. Any software package for computing the inverse will return an error.
I'm guessing that perhaps you just assumed, having incorrectly got MSE=0, that the whole variance estimate must be zero, and you didn't even consider what  $(X'X)^{-1}$ might be equal to.
That reasoning overlooks that $0 \times \infty=$NA.
R does it all correctly
Of course, statistical software will do it all correctly for you. Here is a little example with $n=3$ and $p=4$.
> y
[1] -0.0723 -0.4384  0.1519
> X <- cbind(Intercept=1,x1,x2,x3,x4)
> X
     Intercept     x1     x2      x3     x4
[1,]         1 -0.126  3.070  0.2913 -0.290
[2,]         1 -1.150  0.694 -0.1480 -1.614
[3,]         1  1.862 -1.290  0.0913  0.842

The rank of $X$ is 3, not 5:
> qr(X)$rank
[1] 3

The linear regression calculations are all done correctly:
> summary(lm(y~x1+x2+x3+x4))

Call:
lm(formula = y ~ x1 + x2 + x3 + x4)

Residuals:
ALL 3 residuals are 0: no residual degrees of freedom!

Coefficients: (2 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.2096         NA      NA       NA
x1            0.2317         NA      NA       NA
x2            0.0543         NA      NA       NA
x3                NA         NA      NA       NA
x4                NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:     1,      Adjusted R-squared:   NaN 
F-statistic:  NaN on 2 and 0 DF,  p-value: NA

R tells you explicitly that the residual df is 0 and that MSE$^{1/2}$ (aka "Residual standard error") is NaN. The standard errors for the coefficients $\hat\beta$ are of course all NA (not 0).
