Is there any way to recover a "closed-form-expression" for the cdf of a bivariate normal distribution?
What I have in mind is the following: for a univariate standard normal, we can write $$\text{Pr}(X\leq x)=\dfrac{1}{2}\left(1+\text{erf}\left(\dfrac{x}{\sqrt{2}}\right)\right)$$ in terms of the error function.
When doing the same exercise with a bivariate standard normal with a correlation of, say, $\rho=1/2$, I find that (with $f(x,y)$ being the joint density) $$\int_{-\infty}^{y}f(t,s)\,ds=\dfrac{\exp\left(-\dfrac{t^2}{2}\right)\left(1-\text{erf}\left(\dfrac{\sqrt{6}t}{6}-\dfrac{\sqrt{6}y}{3}\right)\right)}{2\sqrt{2}\pi}$$ such that recovering $$\text{Pr}(X\leq x,Y\leq y)=\int_{-\infty}^x\int_{-\infty}^yf(t,s)\,ds\,dt$$ involves an integration over $$\exp\left(-\dfrac{t^2}{2}\right)\left(1-\text{erf}\left(\dfrac{\sqrt{6}t}{6}-\dfrac{\sqrt{6}y}{3}\right)\right)$$ where I have absolutely no idea how to tackle that. Is there any way to write $\text{Pr}(X\leq x , Y\leq y)$ in terms of error- or other well-known functions?
Thank you very much for your help.