Programming/Kdb/Labs/Option pricing

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

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;T](log[S%K]+((r-q)+.5*sigma*sigma)*T)%sigma*sqrt[T]};

Then,

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