Second moment of the product of normal cdf & pdf

Question

I'm trying to compute the following integral:

$$ I = \int_{-\infty}^\infty x^2 \Phi(a+bx)\phi(x)dx $$

Where $\Phi(\cdot)$ and $\phi(\cdot)$ are the cdf and pdf respectively of the standard normal.

Wikipedia has a list of related integrals (https://en.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions) but not the one I want to compute. I've tried decomposing $I$ into functions from this list, but without success.

There are also some related MathSE/Mathoverflow posts (e.g. Definite integral of a product of normal pdf and cdf), but they only integrate the cdf/pdf product, and not the second moment.

Does anyone have a way to solve this integral analytically? Thanks!

Answer 1

Another approach could be combining methods outlined here integration-of-the-product-of-pdf-cdf-of-normal-distribution and redefining the integral with a more general definition $\int_{-\infty}^{\infty}\phi(cx)\Phi(a+bx)dx$ and differentating with respect to auxillary parameter $c$, evaluating at $c=1$ we obtain the same integral in the question $\int_{-\infty}^{\infty}\phi(x)\Phi(a+bx)x^2dx$. Then final result becomes, \begin{equation}\int_{-\infty}^{\infty}\phi(x)\Phi(a+bx)x^2dx=\phi\left(\frac{a}{\sqrt{1+b^2}}\right)-\frac{ab^2}{(a+b^2)^{3/2}}\phi\left(\frac{a}{\sqrt{1+b^2}}\right) \end{equation} The derivation is below,

\begin{eqnarray} I=\int_{-\infty}^{\infty}\phi(x)\Phi(a+bx)x^2dx\\ =\frac{\partial\left(\int_{-\infty}^{\infty}\phi(cx)\Phi(a+bx)\right)}{\partial c}|_{c=1}\\ =\frac{\partial\left(\int_{-\infty}^{\infty}\phi(y)\Phi(a+\frac{by}{c})\right)}{\partial c}|_{c=1}\\ =\frac{\partial\left(\frac{\Phi\left(\frac{a}{\sqrt{1+\frac{b^2}{c^2}}}\right)}{c}\right)}{\partial c}|_{c=1}\\ =\phi\left(\frac{a}{\sqrt{1+b^2}}\right)-\frac{ab^2}{(a+b^2)^{3/2}}\phi\left(\frac{a}{\sqrt{1+b^2}}\right) \end{eqnarray}

A quick R implementation shows that the numerical and analytical formula coincides for $a=1,b=2,c=1$

pnorm(a/sqrt(1+b^2/c^2))/c^3-dnorm(a/sqrt(1+b^2/c^2))*(a*b^2/(c^3*(1+b^2/c^2)^1.5))/c^2
[1] 0.5434922
integrate(function(x,a,b,c){pnorm(a+b*x)*dnorm(x)*x^2},-Inf,Inf,a=a,b=b,c=1)$value
[1] 0.5434922

Answer 2

\begin{align} I &= \int_{-\infty}^{\infty} x^2 \Phi (a + b x) \phi(x) \mathrm{d}x \\ &= \int_{-\infty}^{\infty} \left[ -\phi(x) \right]' x \Phi (a + b x) \mathrm{d}x \\ &= \left[ -\phi(x) x \Phi (a + b x) \right]_{-\infty}^\infty + \int_{-\infty}^{\infty} \phi(x) \Phi (a + b x) \mathrm{d}x + b \int_{-\infty}^{\infty} x \phi(x) \phi (a + b x) \mathrm{d}x \\ &= 0 + \Phi \left( \frac{a}{\sqrt{1+b^2}} \right) -\frac{a b^2}{\sqrt{2 \pi}(1+b^2)^{3/2}} \exp\left(-\frac{a^2}{2+2b^2} \right) \\ &= \Phi \left( \frac{a}{\sqrt{1+b^2}} \right) -\frac{a b^2}{(1+b^2)^{3/2}} \phi\left(\frac{a}{\sqrt{1+b^2}} \right) \end{align}