Several threads (e.g. [Integration of the product of pdf & cdf of normal distribution](https://mathoverflow.net/questions/101469/integration-of-the-product-of-pdf-cdf-of-normal-distribution) ) have shown that 


$E[\Phi(x)]=\Phi(\mu/\sqrt{\sigma^2+1})$ when $x\sim N(\mu,\sigma^2)$.


I'd like to compute $Var(\Phi(x))$ for such an $x$, ideally without numerically integrating. Does anyone have any ideas?



Variance of the normal CDF
I think I figured it out, actually. Just as for the mean we can reinterpret $E[\Phi(x)]$ as $E[\mathbf 1[y<x]]$ for $y\sim N(0,1)$, we can reinterpret $E[\Phi(x)^2]$ as $E[\mathbf 1[y<x,y^\prime<x]]$ for independent standard normals $y,y^\prime$. Therefore,


$Var(\Phi(x)) = E[\Phi(x)^2]-E[\Phi(x)]^2=\Phi(\frac{\mu}{\sqrt{\sigma^2+1}},\frac{\mu}{\sqrt{\sigma^2+1}},\frac{\sigma^2}{\sigma^2+1})-\Phi(\frac{\mu}{\sqrt{\sigma^2+1}})^2$


Where $\Phi(h,k,r)$ is the bivariate normal CDF with correlation $r$.


