I wish to compute the following integral.
$$ I_1 = \underbrace{\int_{-\infty}^{\infty} \dots\int_{-\infty}^\infty}_{n} \frac{\partial^n}{\partial x_1 \dots \partial x_n}\Bigg(\max \Big( \sum_i^n a_i \max(x_i,0) \Big),0 \Bigg) f(x_1,\dots,x_n) dx_1 \dots dx_n \\ = \prod_i^n a_i \underbrace{\int_{-\infty}^{\infty} \dots\int_{-\infty}^\infty}_{n} \delta^{n - 2} \Big(\sum_i^n a_i \max(x_i,0) \Big) f(x_1,\dots,x_n)u(x_1) \dots u(x_n) dx_1 \dots dx_n \\ = \prod_i^n a_i \underbrace{\int_0^{\infty} \dots\int_0^\infty}_{n} \delta^{n - 2}\left(\sum_i^na_i x_i\right) f(x_1,\dots,x_n) dx_1 \dots dx_n \\ $$
The function $f(x_1,\dots,x_n)$ is the joint gaussian distribution such that $(x_1,\dots,x_n)^\top \sim \mathcal{N}(\mathbf{0},\Sigma)$. Note that it belongs to Schwarz class. Also that $\delta^{n-2}$ is the (n-2)$^{th}$ derivative of the dirac. $u(x_i)$ is the heavyside step function.
I managed to show the following:
Lemma #1: For any even number (n) of jointly gaussian RVs $\sim \mathbb{N}(\mathbf{0},\Sigma)$ where $\mathbb{E}[x_i x_j] = \Sigma(i,j)$ the following is true:
$$ I_2 = \prod_i^n a_i \underbrace{\int_{-\infty}^{\infty} \dots\int_{-\infty}^\infty}_{n} \delta^{n - 2}\left(\sum_i^na_i x_i\right) f(x_1,\dots,x_n) dx_1 \dots dx_n = \prod_{i=0}^{\frac{n}{2}-1}\left(\frac{1}{2} - i\right) \frac{2^{\frac{n}{2}}}{\sqrt{2 \pi}} \prod_i^n a_i \Big(\sum_i^n a_i^2 \Sigma(i,i) + 2 \sum_i^n \sum_j^{i-1}a_i a_j \Sigma(i,j) \Big)^{\frac{1}{2} - \frac{n}{2}} $$
To give you an insight on what this weird looking solution is about, It is in fact the $\frac{n}{2}$ derivatives of the covariances of the following quantity:
$$I_2 = \frac{\partial^{\frac{n}{2}} }{\prod_i^{\frac{n}{2}} \partial \Sigma(i,i+1)} \frac{1}{\sqrt{2 \pi}} \Big(\sum_i^n a_i^2 \Sigma(i,i) + 2\sum_i^n \sum_j^{i-1}a_i a_j \Sigma(i,j) \Big)^{\frac{1}{2}}$$
The question is how to extend $I_2$ to the case where the lower bounds of the integral start from zero instead of $-\infty$ as in $I_1$?.
For further clarification, cosinder the case n = 2. How to find $I_1$ given that $I_2$ is known:
$$ I_1 = a_1 a_2 \int_0^{\infty}\int_0^\infty \delta \left(a_1 x_1 + a_2 x_2\right) f(x_1,x_2) dx_1 dx_2 \\ I_2 = a_1 a_2 \int_{-\infty}^{\infty}\int_{-\infty}^\infty \delta \left(a_1 x_1 + a_2 x_2\right) f(x_1,x_2) dx_1 dx_2 = \frac{a_1 a_2}{\sqrt{2 \pi}} \frac{1}{\sqrt{a_1^2 \Sigma(1,1)^2 + a_2^2 \Sigma(2,2)^2 + 2 a_1 a_2 \Sigma(1,2)}}\\ $$