Thanks a lot for this blog post. My colleagues and I would also like to calculate the gradient of the bivariate standard normal CDF w.r.t. rho, i.e. ∂F/∂ρ. We got: ∂F(x,y,rho)/∂ρ = ∂ ( int_(-infty)^x int_(-infty)^y f( a, b, rho ) db da ) / ∂ρ = int_(-infty)^x int_(-infty)^y ∂ f( a, b, rho ) / ∂ρ ) db da, whereas f( x, y, rho ) is the bivariate standard normal PDF, and we calculated ∂ f( x, y, rho ) / ∂ρ. However, we have not yet managed to solve the double integral so that we can easily compute this integral. Any hints?