Background: the Black-Scholes formulae

Recall the celebrated Black-Scholes equation

{\frac {\partial V}{\partial t}}+{\frac {1}{2}}\sigma ^{2}S^{2}{\frac {\partial ^{2}V}{\partial S^{2}}}+(r-q)S{\frac {\partial V}{\partial S}}-rV=0.

Here

$t$ is a time in years; we generally use $t=0$ as now;
$V(S,t)$ is the value of the option;
$S(t)$ is the price of the underlying asset at time $t$ ;
$\sigma$ is the volatility — the standard deviation of the asset's returns;
$r$ is the annualized risk-free interest rate, continuously compounded;
$q$ is the annualized (continuous) dividend yield.

The solution of this equation depends on the payoff of the option — the terminal condition. In particular, if at the time of expiration, $T$ , the payoff is given by $V(S,T)=C(S,T)=:\max\{S-K,0\}$ , in other words, the option is a European call option, then the value of the option at time $t$ is given by the Black-Scholes formula for the European call:

C(S_{t},t)=e^{-r\tau }[F_{t}N(d_{1})-KN(d_{2})]

where $\tau =T-t$ is the time to maturity, $F=S_{t}e^{(r-q)\tau }$ is the forward price, and

d_{1}={\frac {1}{\sigma {\sqrt {\tau }}}}\left[\ln \left({\frac {S_{t}}{K}}\right)+\left(r-q+{\frac {1}{2}}\sigma ^{2}\right)\tau \right]

and

d_{2}=d_{1}-\sigma {\sqrt {\tau }}.

Here we have used $N(x)$ to denote the standard normal cumulative distribution function,

N(x)={\frac {1}{\sqrt {2\pi }}}\int _{-\infty }^{x}e^{-z^{2}/2}\,dz.

Similarly, if the payoff is given by $V(S,T)=P(S,T)=:\max\{K-S,0\}$ , in other words, the option is a European put option, then the value of the option at time $t$ is given by the Black-Scholes formula for the European put:

P(S_{t},t)=e^{-r\tau }[KN(-d_{2})-F_{t}N(-d_{1})].

We will implement the formulae for the European call and European put in q. However, our first task is to implement $N(x)$ .

Task 1: Implementing the standard normal cumulative distribution function

As mentioned in the Handbook of Mathematical Functions, $N(x)$ can be approximated by

\left\{{\begin{array}{ll}1-\phi (x)\left[c_{1}k+c_{2}k^{2}+c_{3}k^{3}+c_{4}k^{4}+c_{5}k^{5}\right],&x\geq 0,\\1-N(-x),&x<0,\end{array}}\right.

where

\phi (x)=\exp(-x^{2}/2)/{\sqrt {2\pi }},

$k=1/(1+0.2316419x),c_{1}=0.319381530,c_{2}=-0.356563782,c_{3}=1.781477937,c_{4}=-1.821255978,c_{5}=1.330274429.$

Can you implement this function in q?

First, we need

pi:acos -1;

One (terse) implementation would be

normal_cdf:{abs(x>0)-(exp[-.5*x*x]%sqrt 2*pi)*t*.31938153+t*-.356563782+t*1.781477937+t*-1.821255978+1.330274429*t:1%1+.2316419*abs x};

Task 2: Implement the Black-Scholes formula for the European call and put

Equipped with our implementation of normal_cdf, can you implement the Black-Scholes formula for the European call?

First, we need

compute_d1:{[S;K;r;q;sigma;tau](log[S%K]+((r-q)+.5*sigma*sigma)*tau)%sigma*sqrt tau};

Then,

call_price:{[S;K;r;q;sigma;tau]
  d1:compute_d1[S;K;r;q;sigma;tau];
  d2:d1-sigma*sqrt tau;
  F:S*exp tau*r-q;
  (exp neg r*tau)*(F*normal_cdf d1)-K*normal_cdf d2};

and

put_price:{[S;K;r;q;sigma;tau]
  d1:compute_d1[S;K;r;q;sigma;tau];
  d2:d1-sigma*sqrt tau;
  F:S*exp tau*r-q;
  (exp neg r*tau)*(K*normal_cdf neg d2)-F*normal_cdf neg d1};

We can test these implementations on a few sets of parameters, for example

call_price[100f;105f;.05;.07;.1;.5]

gives 0.7991363, whereas

put_price[100f;105f;.05;.07;.1;.5]

gives 6.646135.

Background: a simple Monte Carlo model

The Black-Scholes equation may not have analytic solutions for all derivatives that we are interested in. However, it may still be possible to solve it numerically using Monte Carlo methods.

The model for stock price evolution is