# Expected option return in MATLAB

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

• You can use Gaussian quadrature methods. – user16651 Jun 24 '16 at 13:44
• $$d_2^*=d_1^*-\sqrt{\tau}\sigma$$ – user16651 Jun 24 '16 at 13:59
• @BehrouzMaleki Thanks for correcting the typo, this was only in the question text but not in the MATLAB code of function blsprice. – emcor Jun 24 '16 at 14:39
• The function "normpdf(log(S_T),mu-d-sigma^2/2,sigma),0,inf)" should likely be " normpdf(log(S_T),(mu-d-sigma^2/2) T,sigma),0,inf)". – Gordon Jun 24 '16 at 20:58
• @Gordon Yes that $T$ could be added. I had $T=K=S=1$ set in my code, but the results are still unequal. I believe it is due to numerical errors in the integration. – emcor Jun 24 '16 at 21:29

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.

• I checked that indeed using the lognormal-distribution solves the error. – emcor Jun 30 '16 at 12:32
• EOR_MATLAB = integral(@(S_T)max(S_T-K,0).*normpdf(log(S_T),mu-d-sigma^2/2,sigma)./S_T,0,inf)/blsprice(S, K, r, T,sigma,d)-1 – emcor Jun 30 '16 at 12:32
• Glad I could help – Quantuple Jun 30 '16 at 13:01