Expectation of Maximum of Translated Gaussian Given $X\propto N(\mu,\psi^2)$ and $K$ constant, find $E(\max\{K-X,0\})$.
My attempt:
$$
\begin{align}
E(\max\{K-X,0\})\\
&=\int_{-\infty}^{\infty}\max\{K-x,0\}f_X(x)dx\\
&=\int_{K}^{\infty}(K-x)f_X(x)dx\\
&=\int_{K}^{\infty}[(K-\mu)-(X-\mu)]f_X(x)dx\\
&=\int_{K}^{\infty}(K-\mu)f_X(x)dx-\int_{K}^{\infty}(x-\mu)f_X(x)dx\\
&=(K-\mu)(1-\Phi(K))-\int_{K}^{\infty}(x-\mu)f_X(x)dx\\
&=(K-\mu)(1-\Phi(K))-\int_{K}^{\infty}(x-\mu)\frac{1}{(2\pi \psi^2)^{1/2}}e^{-\frac{(x-\mu)^2}{2\psi^2}}dx\\
&=(K-\mu)(1-\Phi(K))-\int_{\frac{(K-\mu)^2}{2\psi^2}}^{\infty}\psi^2\frac{1}{(2\pi \psi^2)^{1/2}}e^{-z}dz\\
&=(K-\mu)(1-\Phi(K))+\frac{\psi}{\sqrt{2\pi}}e^{-z}|_{\frac{(K-\mu)^2}{2\psi^2}}^{\infty}\\
&=(K-\mu)(1-\Phi(K))-\frac{\psi}{\sqrt{2\pi}}e^{-\frac{(K-\mu)^2}{2\psi^2}}
\end{align}
$$
where I have used the change of variables $z=\frac{(x-\mu)^2}{2\psi^2}$.
However, when I check my answer numerically, I see that the integral is not correct. Here is working Python code to check:
import numpy as np
from scipy.stats import norm

# simulation parameters for checking answers numerically
trials = 100000
K = 0.25
mu = np.pi/2
psi = 3.0

simulated = np.mean([np.max([K-x, 0])
                     for x in np.random.normal(mu, psi, trials)])
numerical = (K-mu)*(1-norm.cdf(K, mu, psi)) - np.exp(-((K-mu)**2)/(2*psi**2))*psi/np.sqrt(2*np.pi) 
print(simulated)
print(numerical)


This gives:
simulated = 0.6448298784565597
numerical = -1.9713797586737627
What is the issue with my integration?
 A: The bounds on the original integral are not correct. The requirement $x\le K$ means the bounds must run from $-\infty$ to $K$.
Here are the correct steps:
$$
\begin{align}
E(\max\{K-X,0\})\\
&=\int_{-\infty}^{\infty}\max\{K-x,0\}f_X(x)dx\\
&=\int_{-\infty}^{K}(K-x)f_X(x)dx\\
&=\int_{-\infty}^{K}[(K-\mu)-(X-\mu)]f_X(x)dx\\
&=\int_{-\infty}^{K}(K-\mu)f_X(x)dx-\int_{-\infty}^{K}(x-\mu)f_X(x)dx\\
&=(K-\mu)\Phi(K)-\int_{-\infty}^{K}(x-\mu)f_X(x)dx\\
&=(K-\mu)\Phi(K)-\int_{-\infty}^{K}(x-\mu)\frac{1}{(2\pi \psi^2)^{1/2}}e^{-\frac{(x-\mu)^2}{2\psi^2}}dx\\
&=(K-\mu)\Phi(K)-\int_{\infty}^{\frac{(K-\mu)^2}{2\psi^2}}\psi^2\frac{1}{(2\pi \psi^2)^{1/2}}e^{-z}dz\\
&=(K-\mu)\Phi(K)+\frac{\psi}{\sqrt{2\pi}}e^{-\frac{(K-\mu)^2}{2\psi^2}}
\end{align}
$$
Here is working Python code to check:
import numpy as np
from scipy.stats import norm

# simulation parameters for checking answers numerically
trials = 100000
K = 0.25
mu = np.pi/2
psi = 3.0

simulated = np.mean([np.max([K-x, 0])
                     for x in np.random.normal(mu, psi, trials)])
numerical = (K-mu)*(norm.cdf(K, mu, psi)) + np.exp(-((K-mu)**2)/(2*psi**2))*psi/np.sqrt(2*np.pi) 
print(simulated)
print(numerical)

This gives:
simulated = 0.6526939144762599
numerical = 0.650583431878866
where the difference in value arises because of the finite number of trials used for the simulation.
