Truncated/Censored Regression in the Multivariate Normal Case Assume that a $n$-dimensional vector of real valued variables $\mathbf{Y}=(Y_1,...,Y_n)^T$ is given. We know that they are jointly normal which expectation vector $\boldsymbol{\mu}=(\mu_1,...,\mu_n)^T$  and Variance-Covariance Matrix $\boldsymbol{\Sigma}\neq I_{n}\sigma^2$ (i.e. the observations are dependent). So we have
$$ \mathbf{Y}\sim N_n(\boldsymbol{\mu},\boldsymbol{\Sigma})$$.
Now assume furthermore that there is a share of censored observations (i.e. we know they are smaller than zero but only observe a zero)
$\mathbf{Y}_m=(0,..,0)^T\subset \mathbf{Y}$ and a share of observed observations $\mathbf{Y}_o$. Hence we have $\mathbf{Y}=(\mathbf{Y}_o,\mathbf{Y}_m)^T$
Now if i know that the expectation is a function of some parameters, like $\boldsymbol{\mu}=(X_1^T\beta,...,X_n^T\beta)^T$ and also know the structure of the covariance $\boldsymbol{\Sigma}=\boldsymbol{\Sigma}(\lambda)$  how would I estimate those parameters $\theta=(\beta,\lambda)^T$?
I thought about the following in order to form the likelihood:
$$ f(\mathbf{Y},X,\theta)=  f(\mathbf{Y}_o,X_o,\theta)f(\mathbf{Y}_m,X_m,\theta|\mathbf{Y}_o,\mathbf{Y}_m<0)$$
i.e. I would take a multivariate normal density for the uncensored part ($f(\mathbf{Y}_o,X_o,\theta)$) and then estimate the probability of a multivariate normal cdf with conditional mean ($E[\mathbf{Y}_m|\mathbf{Y}_o,X]$) and variance ($V[\mathbf{Y}_m|\mathbf{Y}_o,X]$) that all $\mathbf{Y}_m$ are smaller than zero. (Conditional Mean and Variance are easily found by applying a partition as in http://www.maths.manchester.ac.uk/~mkt/MT3732%20(MVA)/Notes/MVA_Section3.pdf)
Although I thought that this is a straight generalization of the standard censored regression model with independent observations (see for example https://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf), I have not found any paper or implementation.
Hence my questions: Do you know a solid way to solve this problem? Is there even some implementation, for example in R?
 A: It is not so straightforward when responses are correlated, it turns out. 
The method described in this paper: https://arxiv.org/pdf/1808.10541.pdf is effectively what you are looking for (except they deal with right-censoring, but the modification is straightforward, they suggest). The citation for this paper (per gung's request) is: Molstad, A. J., Hsu, L., and Sun, W. (2018). Gaussian process regression for survival time prediction with genome-wide gene expression. arXiv preprint, arXiv:1808.10541.
The authors show that you can fit such a model using a Monte-Carlo expectation conditional maximization algorithm. The complication is that in performing the $E$-step, the distribution of the censored outcomes is intractable due to the dependence across censored responses, so they use a Gibbs sampler to approximate the $Q$-function (that is, the expected value of the complete data log-likelihood). To modify their method for your setup, you would need to (i) change the hypercube over which they draw samples for the approximate $E$-step, and (ii) change their conditional $M$-step for the variance components according to whatever parameterization you use for your covariance matrix. 
