How to calculate the following expectation and variance?

Question

Given $$X_1 \sim N(\mu_1, \sigma_1^2),$$ $$X_2 \sim N(\mu_2, \sigma_2^2),$$ $$X_1 \mathrm{\ and\ } X_2 \mathrm{\ are\ independent}.$$ and $$Y = \max(\max(X_1,0) + X_2, 0),$$ find $$E(Y) \mathrm{\ \ and\ \ } Var(Y).$$ Python expression is preferred.

---

I tried using Maple to calculate the expectation and convert the solution to Python but the solution is wrong and not runnable in Python. Here is the Maple script:

with(CodeGeneration):
p1 := 1/(sqrt(2pi)s1) exp(-(x1-m1)^2 / (2s1^2));
p2 := 1/(sqrt(2pi)s2) exp(-(x2-m2)^2 / (2s2^2));

y1 := piecewise(x1<=0, 0, x1>0, x1);
y2 := y1 + x2;
y3 := piecewise(y2<=0, 0, y2>0, y2);
y3_int := y3 * p1 * p2;

F := int(y3_int, [x1=-infinity..infinity, x2=-infinity..infinity]);
F := optimize(F)
Python(F);

Here is the result:

cg = 0.1e1 / pi / s1 / s2 * sympy.integrate(math.exp(-(x2 - m2) ** 2 / s2 ** 2 / 2) * sympy.integrate(0 if 0 if x1 <= 0 else x1 if 0 < x1 else 0 + x2 <= 0 else 0 if x1 <= 0 else x1 if 0 < x1 else 0 + x2 if 0 < 0 if x1 <= 0 else x1 if 0 < x1 else 0 + x2 else 0 * math.exp(-(x1 - m1) ** 2 / s1 ** 2 / 2), x1 == xrange(-infinity,infinity)), x2 == xrange(-infinity,infinity)) / 2

---

I also tried using Sympy to calculate the expectation with the following script. However, it is very slow. It has run for 3 hours but no result is returned.

import math
import sympy as sp


m1 = sp.Symbol('m1')
s1 = sp.Symbol('s1')
m2 = sp.Symbol('m2')
s2 = sp.Symbol('s2')

x1 = sp.Symbol('x1')
x2 = sp.Symbol('x2')

p1 = 1 / (sp.sqrt(2 * sp.pi) * s1) * sp.exp(-(x1-m1)**2 / (2 * s1**2))
p2 = 1 / (sp.sqrt(2 * sp.pi) * s2) * sp.exp(-(x2-m2)**2 / (2 * s2**2))

y1 = sp.Piecewise((0, x1 <= 0), (x1, True))
y2 = y1 + x2
y3 = sp.Piecewise((0, y2 <= 0), (y2, True))
y3_int = y3 * p1 * p2

expe = sp.integrate(y3_int, (x1, -sp.oo, sp.oo), (x2, -sp.oo, sp.oo))

---

Any help would be appreciated!

Answer 1

Here's as far as I could get. To simplify things let $f_1$ and $f_2$ denote the probability density functions of $X_1$ and $X_2$. Let $X$ denote a random variable with the normal distribution $N(0,1)$. $$E(Y)=\int_0^{\infty}\int_{-x_1}^{\infty}(x_1+x_2)f_1(x_1)f_2(x_2)\:dx_2dx_1 + \int_{-\infty}^0\int_0^{\infty}x_2f_2(x_2)\:dx_2f_1(x_1)dx_1$$ because there are two cases where $Y$ is not $0$.

$1)$ $X_1 > 0$ and $X_2 > -X_1$

$2)$ $X_1\leq 0$ but $X_2 > 0$

The second term simplifies to $$ P(X_1\leq 0)\int_0^{\infty}x_2f_2(x_2)\:dx_2=P(X_1\leq 0)\Big[\int_{-\mu_2/\sigma_1}^{\infty}\frac{\sigma_2}{\sqrt{2\pi}}x_2e^{-\frac{1}{2}x_2^2}dx_2+\mu_2P(X\geq \frac{-\mu_2}{\sigma_2})\Big]$$ $$P(X_1\leq 0)\Big[\frac{\sigma_2}{\sqrt{2\pi}}e^{-\frac{\mu_2^2}{2\sigma_2^2}}+\mu_2P(X\geq \frac{-\mu_2}{\sigma_2})\Big]$$

The first term simplifies to $$\int_0^{\infty}\int_0^{\infty}(x_1+x_2)f_1(x_1)f_2(x_2)\:dx_1dx_2+ \int_0^{\infty}\int_{-x_1}^{0}(x_1+x_2)f_1(x_1)f_2(x_2)\:dx_2dx_1$$ the first term of which simplifies to $$P(X_2\geq 0)\Big[\frac{\sigma_1}{\sqrt{2\pi}}e^{-\frac{\mu_1^2}{2\sigma_1^2}}+\mu_1P(X\geq \frac{-\mu_1}{\sigma_1})\Big]+P(X_1\geq 0)\Big[\frac{\sigma_2}{\sqrt{2\pi}}e^{-\frac{\mu_2^2}{2\sigma_2^2}}+\mu_2P(X\geq \frac{-\mu_2}{\sigma_2}\Big]$$ For given means and variances the probabilities would not be hard to find using a computer and $z$-scores. As for the remaining term, I think it's possible to simplify it, but I haven't tried since the work required to do so is enormous.