Expected option return in MATLAB

Question

The expected return of an option is given by its expected payoff under $P$ over its market price under $Q$.

For the Black-Scholes model, expected call option return is given as (see here):

$$ E(R)=\frac{E^P[(S_T-K)^+]}{e^{-rT}E^Q[(S_T-K)^+]}=\frac{e^{\mu \tau}[S_tN(d_1^*)-e^{\mu \tau}KN(d_2^*)]}{C_t(r,T,\sigma,S,K)}-1 $$

$$\text{with }d_1^*=\frac{\ln S_t/K+(\mu+\frac{1}{2}\sigma^2)\tau}{\sigma\sqrt{\tau}},\qquad d_2^*=d_1^*-\sqrt{\tau}\sigma$$

I implemented the $P$-payoff in MATLAB as

    E(R) = exp((mu-d)*T)*blsprice(S, K, mu, T, sigma,d)

and get correct values (comparing with other studies).

However, I also tried to calculate out the expectation integral numerically in MATLAB as follows:

     E(R2) = integral(@(S_T)max(S_T-K,0).*normpdf(log(S_T),mu-d-sigma^2/2,sigma),0,inf)

(with some arbitrary parameters) and I get a different value.

Can someone explain whether there is error in my code for E(R2), or is MATLAB integration just not accurate enough?

E.g. try with

 S=1,T=1,K=1,r=0.01,mu=0.1,sigma=0.05,d=0.02

Answer 1

Under GBM $$ \frac {dS_t}{S_t} = \mu dt + \sigma dW_t $$ we get $$ S_T = S_0 e^{(\mu - \frac{1}{2}\sigma^2)T + \sigma W_T} $$ suggesting that $$ S_T \sim \text{ln}\mathcal {N} ( \tilde {\mu}, \tilde {\sigma}) $$ where \begin{align} \tilde {\mu} &= \ln S_0 + (\mu - \frac{1}{2}\sigma^2)T \\ \tilde {\sigma} &= \sigma \sqrt {T} \end{align}

Now if $X \sim \text{ln}\mathcal {N} (\mu, \sigma)$, the pdf of $X $ reads $$ p (x) = \frac {1}{x \sigma \sqrt {2\pi}} e^{-\frac {(\ln x - \mu)^2}{2\sigma^2}} $$ showing that the lognormal pdf relates to the normal pdf as follows $$ \text {lognormpdf} (x, \mu, \sigma) = \frac {\text {normpdf}(\ln (x), \mu, \sigma)}{x} $$

So finally: $$ I = \mathbb {E}_0 [(S_T-K)^+] = \int_0^\infty \max(S_T-K,0) p (S_T) dS_T $$ should be calculated as

mu_tilde = log(S_0) + (mu - 0.5*sigma^2)*T
sigma_tilde = sigma*sqrt(T)
I = integral(@(S_T) max(S_T-K,0) .* 1./S_T.*normpdf(log(S_T),mu_tilde,sigma_tilde),0,inf)

So you basically have 3 problems, 2 of which are in the expressions of the drift/diffusion coefficients (but these will remain hidden when picking $S_0=T=1$) and the main one being the missing factor $1/S_T $ to get the lognormal pdf.