Expectation of Maximum of Translated Gaussian

Question

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?

Answer 1

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.