Covariance matrix of multivariate normal when negative values are made zero

Question

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$?

Answer 1

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.