MLE when variance of residuals is null (y is a linear combination of x) Suppose I have the following model to be estimated via MLE assuming normal errors $y_{t}=x_{t} \beta +e_{t}$ with $e=N(0,\sigma^{2})$, where $y, x$ are matrixes and $\beta$ is a vector, so $\sigma^{2}$ is the the covariance matrix of residuals. We have to estimate $\beta$ and $\sigma ^{2}$. So the loglikelihood has the traditional shape where we find $log(|\sigma^{2}|)$ and the term $(y_{t}-x_{t} \beta)^{T}\sigma^-2(y_{t}-x_{t} \beta)$, so that $\sigma^2$ must be positive definite otherwise its determinant is 0 and it is singular.  
Suppose now that $y$ is exactly replicated by $\beta x$ ,i.e. each row of $y$ is a linear combination of the corresponding row in $x$ (for example assume $y=3*Ix$ where $I$ is the unit matrix). Then numerically the MLE solver will likely fail because as the estimated $\beta$ converges to the true $\beta$ then the estimated $\sigma^{2}$ goes toward the null matrix, which is not pd and therefore not invertible. However, I am doing nothing wrong, as I am finding the true coefficients of the data generating process because, as said, y is a linear combination of x. For example in OLS regression, where we just want to minimize the sum of squared residuals, then this does not create any problem (simply the sum of squared residuals will go to 0!), while it does for the MLE.
So my questions are:
1) is it analytically wrong to assume ex-ante that $\sigma^{2}$ is unknown and, in the end, put estimated $\sigma^{2}=0$? Is it right to do so even if, for estimated $\sigma^{2}=0$, the PDF of the likelihood does not exist?
2) How do you solve this problem in practice from the numerical point of view? Especially from the solver point of view?
Thanks
 A: From the analytical perspective no problem. Indeed, if $y_{t} $ is a linear combination of $X_{t} $ with population coefficients $\beta^{*}$ for each t, then, assuming to minimize the minus log likelihood and taking the gradient of the - log likelihood with respect to $\beta$ and setting it equal to zero will yield (after the due simplifications) the first order condition for each t (I will skip the t suffixes):
$$ y = X \beta$$
which, since $y$ is a linear combination of $X$ with parameters $\beta^{*}$ for each t, is equivalent to $y=X\beta^{*}=X \beta$, where $\beta$ is the unknown of the set of equations. As you see, the set of equations is is easily solved for $\beta=\beta^{*}$, which is the maximum likelihood estimate: indeed, it is easily seen that if you set $\beta$ at the "true" value $\beta^{*}$ such that $y-X\beta^{*} =0$ (where $0$ denotes a vector here), then the minus log likelihood function is minimized (i.e. you will have a global minimum with respect to the parameter $\beta$). 
Then, taking the gradient with respect to $\Sigma$ and setting it equal to 0, will yield $\Sigma= \sum_{t} (y- X\beta^{*} )(y-X\beta^{*} )^{T}$, where we know from the point above that $y-X\beta^{*} =0$ for each t, which implies the solution $\Sigma^{*}=0$ (where 0 denotes a matrix here).
For the numerical issues with the available solvers, see for example pag.11 of this source relating to MARSS models. In particular see pag. 11 where they talk about degenerate covariance matrixes and invite to use the "degenerate" option (allow.degen=TRUE), which is also allowed by Python Scipy library when dealing with Multivariate Normal Density (option "allow_singualr=True").
