Following up on the question (and answers) here, I'm trying to derive $\frac{\partial \Phi(x_1, x_2|\mathbf{\underline{\theta}})}{\partial x_1}$ and $\frac{\partial \Phi(x_1, x_2|\mathbf{\underline{\theta}})}{\partial x_2}$ where $\mathbf{\underline{\theta}} = (\mu_1, \mu_2, \sigma_1^2, \sigma_2^2, \rho)$ and $\Phi$ denotes the bivariate normal CDF under examination.
To clarify, I'm not looking for the final answer/expression, the answer is already available online. I just want to know if I'm deriving it correctly, and if not, where am I going wrong?
My work is as follows:
$$ F(x_1, x_2|\mathbf{\underline{\theta}}) = \int_{-\infty}^{x_1}\int_{-\infty}^{x_2} \phi(a, b|\mathbf{\underline{\theta}}) \;db\;da $$
$$ = \int_{-\infty}^{x_1}\int_{-\infty}^{x_2} \frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}.exp \left\{ -\frac{1}{2(1-\rho^2)}\left[\left(\frac{a-\mu_1}{\sigma_1}\right)^2-2\rho\left(\frac{a-\mu_1}{\sigma_1}\right)\left(\frac{b-\mu_2}{\sigma_2}\right)+\left(\frac{b-\mu_2}{\sigma_2}\right)^2\right] \right\}\;db\;da $$
$$ = \int_{-\infty}^{x_1}\frac{1}{\sqrt{2\pi}\sigma_1}.exp \left\{ -\frac{1}{2(1-\rho^2)}\left(\frac{a-\mu_1}{\sigma_1}\right)^2\right\}\int_{-\infty}^{x_2}\frac{1}{\sqrt{2\pi}\sigma_2\sqrt{1-\rho^2}}.exp \left\{ -\frac{1}{2(1-\rho^2)}\left[-2\rho\left(\frac{a-\mu_1}{\sigma_1}\right)\left(\frac{b-\mu_2}{\sigma_2}\right)+\left(\frac{b-\mu_2}{\sigma_2}\right)^2\right] \right\}\;db\;da $$
Let:
$$ \left(\frac{a-\mu_1}{\sigma_1}\right)=X \;;\;\left(\frac{b-\mu_2}{\sigma_2}\right)=Y $$
Therefore the inner integral is now:
$$ = \int_{-\infty}^{x_2}\frac{1}{\sqrt{2\pi}\sigma_2\sqrt{1-\rho^2}}.exp \left\{ -\frac{1}{2(1-\rho^2)}\left(-2\rho XY+Y^2\right) \right\}\;db\;da $$
Focusing on the inner integral, specifically the term inside the exponential; completing the square yields:
$$ exp \left\{ -\frac{1}{2(1-\rho^2)}\left(Y-\rho X\right)^2 \right\}.exp \left\{ \frac{1}{2(1-\rho^2)}\left(\rho X\right)^2 \right\} $$
Taking out the second term in the product (which is constant in Y) yields:
$$ \int_{-\infty}^{x_1}\frac{1}{\sqrt{2\pi}\sigma_1}.exp \left\{ -\frac{1}{2(1-\rho^2)}X^2\right\}.exp \left\{ \frac{1}{2(1-\rho^2)}\left(\rho X\right)^2 \right\}\int_{-\infty}^{x_2}\frac{1}{\sqrt{2\pi}\sigma_2\sqrt{1-\rho^2}}.exp \left\{ -\frac{1}{2(1-\rho^2)}\left(Y-\rho X\right)^2 \right\}\;db\;da $$
$$ = \int_{-\infty}^{x_1}\phi(a)\int_{-\infty}^{x_2}\frac{1}{\sqrt{2\pi}\sigma_2\sqrt{1-\rho^2}}.exp \left\{ -\frac{1}{2(1-\rho^2)}\left(Y-\rho X\right)^2 \right\}\;db\;da $$
Undoing our compact notation in the inner integral yields:
$$ = \int_{-\infty}^{x_1}\phi(a)\int_{-\infty}^{x_2}\frac{1}{\sqrt{2\pi}\sigma_2\sqrt{1-\rho^2}}.exp \left\{ -\frac{1}{2}\left(\frac{b-(\mu_2+\frac{\sigma_2}{\sigma_1}a\rho-\frac{\sigma_2}{\sigma_1}\mu_1\rho)}{\sigma_2\sqrt{1-\rho^2}}\right)^2 \right\}\;db\;da $$
Rewriting the inner integral as a CDF yields:
$$ = \int_{-\infty}^{x_1}\phi_{X_1}(a)\Phi_{X_2^*}(x_2|a)\;da $$
Note that I'm using * to denote the fact that the CDF has a different mean and variance than the marginal cdf of $X_2$.
Now here's where things get a little confusing for me. The answer I referenced at the beginning of this question implies that I need to use Leibniz's rule for differentiating under the integral sign to find $\frac{\partial \Phi(x_1, x_2|\mathbf{\underline{\theta}})}{\partial x_1}$. While I am familiar with Leibniz's rule, I can't help shake this feeling that I'm applying it here incorrectly. Here's the rest of my work:
$$ \frac{\partial \Phi(x_1, x_2|\mathbf{\underline{\theta}})}{\partial x_1}=\frac{\partial}{\partial x_1} \int_{-\infty}^{x_1}\phi_{X_1}(a)\Phi_{X_2^*}(x_2| a)\;da $$
$$ = \int_{-\infty}^{x_1}\frac{\partial}{\partial x_1}\phi_{X_1}(a)\Phi_{X_2^*}(x_2| a)\;da + \phi_{X_1}(x_1)\Phi_{X_2^*}(x_2| x_1)\frac{\partial x_1}{\partial x_1} - \phi_{X_1}(x_1)\Phi_{X_2^*}(x_2| (-\infty))\frac{\partial (-\infty)}{\partial x_1} $$
$$ = \int_{-\infty}^{x_1}\frac{\partial}{\partial x_1}\phi_{X_1}(a)\Phi_{X_2^*}(x_2| a)\;da + \phi_{X_1}(x_1)\Phi_{X_2^*}(x_2| x_1).1 - \phi_{X_1}(x_1)\Phi_{X_2^*}(x_2| (-\infty)).0 $$
$$ = \int_{-\infty}^{x_1}\frac{\partial}{\partial x_1}\phi_{X_1}(a)\Phi_{X_2^*}(x_2| a)\;da + \phi_{X_1}(x_1)\Phi_{X_2^*}(x_2| x_1) $$
Now I know the answer is:
$$ \phi_{X_1}(x_1)\Phi_{X_2^*}(x_2| x_1) $$
However, the only way I can think of to zero the first term is:
$$ \int_{-\infty}^{x_1}\frac{\partial}{\partial x_1}\phi_{X_1}(a)\Phi_{X_2^*}(x_2| a)\;da $$
$$ = \int_{-\infty}^{x_1}0\;da = 0 $$
However, the above feels incredibly wrong to me. My source of confusion arises from the integration dummy, a. I'm not sure if I should treat a as $x_1$ when taking the derivative with respect to $x_1$ or just treat it as a constant. Any help is appreciated.
