# Covariance matrix of multivariate normal when negative values are made zero

Let $$x$$ be $$n$$ dimensionally multivariate normally distributed with mean $$\mu$$ and covariance matrix $$\Sigma$$. Now let $$y$$ be random variables defined by

$$\begin{equation} y_i= \begin{cases} 0, & x_i\leq 0 \\ x_i, & x_i>0 \end{cases} \end{equation}$$

What is the mean and covariance matrix of $$y$$? Can we get an analytic solution? If not, what about in some special cases. For example, where $$\mu=\bf{0}$$ and $$\Sigma$$ has a diagonal of one and all other entries the same value $$\rho$$?

I think just dealing with the bivariate case is all that is necessary as you're only interested in covariances. I also think that there is likely no analytic solution when the means of the $$X$$'s are not zero. Here is an approach using Mathematica:

For the bivariate case $$Y_1$$ and $$Y_2$$ both have means that can be calculated by integrating over 0 to $$\infty$$ using the marginal densities of $$X_1$$ and $$X_2$$, respectively.

Find the means of $$Y_1$$ and $$Y_2$$:

mean1 = Integrate[y1 PDF[NormalDistribution[0, σ1], y1], {y1, 0, ∞}, Assumptions -> σ1 > 0]
(* σ1/Sqrt[2 π] *)

mean2 = Integrate[y2 PDF[NormalDistribution[0, σ2], y2], {y2, 0, ∞}, Assumptions -> σ2 > 0]
(* σ2/Sqrt[2 π] *)


Now the covariance depending on if $$\rho$$ is positive or negative:

pdf = PDF[BinormalDistribution[{0, 0}, {σ1, σ2}, ρ], {y1, y2}];

covρPositive = FullSimplify[Integrate[y1 y2 pdf, {y1, 0, ∞}, {y2, 0, ∞},
Assumptions -> {σ1 > 0, σ2 > 0, 0 < ρ < 1}] - mean1 mean2,
Assumptions -> {σ1 > 0, σ2 > 0, 0 < ρ < 1}]
(* (σ1 σ2 (-1 + π ρ + Sqrt[1 - ρ^2] - ρ ArcCos[ρ]))/(2 π) *)

covρNegative = Integrate[y1 y2 pdf, {y1, 0, ∞}, {y2, 0, ∞},
Assumptions -> {σ1 > 0, σ2 > 0, -1 < ρ <= 0}] - mean1 mean2 // FullSimplify
(* (σ1 σ2 (-1 + π ρ + Sqrt[1 - ρ^2] - ρ ArcCos[ρ]))/(2 π) *)


As a check one can set values for $$\sigma_1$$, $$\sigma_2$$, and $$\rho$$ and take random samples:

n = 10000000;
parms = {σ1 -> 1, σ2 -> 7, ρ -> -1/2};
x1x2 = RandomVariate[BinormalDistribution[{0, 0}, {σ1, σ2}, ρ] /. parms, n];
x1x2[[All, 1]] = Max[0, #] & /@ x1x2[[All, 1]];
x1x2[[All, 2]] = Max[0, #] & /@ x1x2[[All, 2]];

(* Sample estimate *)
Covariance[x1x2][[1, 2]]
(* -0.7329134893798569 *)

(* True covariance formula *)
(σ1 σ2 (-1 + π ρ + Sqrt[1 - ρ^2] - ρ ArcCos[ρ]))/(2 π) /. parms // N
(* -0.7325923679884647 *)


The results are consistent for this example.