How do I take the partial derivative of bivariate normal cdf and bivariate normal pdf with its arguments (i.e. $x_{1}$ ,$x_{2}$ , and $\rho$ in the following equations)?
\begin{equation} y=\Phi(x_{1},x_{2},\rho) \end{equation}
\begin{equation} z=\phi(x_{1},x_{2},\rho) \end{equation}
where $x_{1}$ is normally distributed with mean 0 and variance 1 and $x_{2}$ is normally distributed with mean 0 and variance 1. $\rho$ is the correlation between $x_{1}$ and $x_{2}$.
\begin{align} y &= \Phi(x_1,x_2,\rho) = \int_{-\infty}^{x_1}\left[\int_{-\infty}^{x_2} \phi(a,b,\rho)\,\mathrm db\right]\,\mathrm da\\ \frac{\partial y}{\partial x_1} &= \frac{\partial}{\partial x_1}\Phi(x_1,x_2,\rho) = \frac{\partial}{\partial x_1}\int_{-\infty}^{x_1} \left[\int_{-\infty}^{x_2} \phi(a,b,\rho)\,\mathrm db\right]\,\mathrm da\\ &= \int_{-\infty}^{x_2} \phi(x_1,b,\rho)\,\mathrm db \end{align} via the rule for differentiating under the integral sign. Similarly, $$\frac{\partial y}{\partial x_2} = \int_{-\infty}^{x_1} \phi(a,x_2,\rho)\,\mathrm da.$$ If you don't recall the rule for differentiating integrals, see for example, the comments following this answer on math.SE.
The derivative with respect to $\rho$ is straightforward to find but messy in its details. We have that $$\phi(x_1,x_2,\rho)=\frac{1}{2\pi\sqrt{1-\rho^2}} \exp\left[-\frac{x^2 -2\rho xy + y^2}{2(1-\rho^2)}\right]$$ whose partial derivative with respect to $\rho$ is left to the OP to find. If $g(x_1,x_2,\rho)$ denotes this partial derivative, then $$\frac{\partial}{\partial \rho}\Phi(x_1,x_2,\rho) = \frac{\partial}{\partial \rho}\int_{-\infty}^{x_1} \left[\int_{-\infty}^{x_2} \phi(a,b,\rho)\,\mathrm db\right]\,\mathrm da = \int_{-\infty}^{x_1} \int_{-\infty}^{x_2} g(a,b,\rho)\,\mathrm db\,\mathrm da$$
While the accepted answer provides great tips for somebody looking to re-derive these quantities, for some work I was doing I actually just needed to know the closed-form partial derivatives of the bivariate normal cdf. Further, the integral for the partial with respect to $\rho$ is non-trivial. Since my google search brought me here (and then down a rabbit hole of statistics papers from decades ago), I thought I'd share my findings for the rest of the internet community:
\begin{align*} \frac{\partial}{\partial x_1} \Phi(x_1, x_2, \rho) &= \phi(x_1)\Phi\bigg(\frac{x_2-\rho x_1}{\sqrt{1-\rho^2}}\bigg), \\ \frac{\partial}{\partial x_2} \Phi(x_1, x_2, \rho) &= \phi(x_2)\Phi\bigg(\frac{x_1-\rho x_2}{\sqrt{1-\rho^2}}\bigg),~\text{and} \\ \frac{\partial}{\partial \rho} \Phi(x_1, x_2, \rho) &= \frac{1}{2\pi\sqrt{1-\rho^2}}\exp\bigg\{\frac{-(x_1^2-2\rho x_1x_2+x_2^2)}{2(1-\rho^2)}\bigg\}. \end{align*}
Here, $\phi(\cdot)$ and $\Phi(\cdot)$ are the standard normal pdf and cdf, respectively. The literature reference for the partial with respect to $\rho$ is:
Z. Drezner and G.O. Wesolowsky (1990). On the computation of the bivariate normal integral. Journal of Statistical Computation and Simulation 35 (1-2): 101-107. Equation (4).
