48
$\begingroup$

Let the Black-Scholes formula be defined as the function $f(S, X, T, r, v)$.

I'm curious about functions that are computationally simpler than the Black-Scholes that yields results that approximate $f$ for a given set of inputs $S, X, T, r, v$.

I understand that "computationally simpler" is not well-defined. But I mean simpler in terms of number of terms used in the function. Or even more specifically, the number of distinct computational steps that needs to be completed to arrive at the Black-Scholes output.

Obviously Black-Scholes is computationally simple as it is, but I'm ready to trade some accuracy for an even simpler function that would give results that approximate B&S.

Does any such simpler approximations exist?

$\endgroup$
39
$\begingroup$

This is just to expand a bit on vonjd's answer.

The approximate formula mentioned by vonjd is due to Brenner and Subrahmanyam ("A simple solution to compute the Implied Standard Deviation", Financial Analysts Journal (1988), pp. 80-83). I do not have a free link to the paper so let me just give a quick and dirty derivation here.

For the at-the-money call option, we have $S=Ke^{-r(T-t)}$. Plugging this into the standard Black-Scholes formula $$C(S,t)=N(d_1)S-N(d_2)Ke^{-r(T-t)},$$ we get that $$C(S,t)=\left[N\left(\frac{1}{2}\sigma\sqrt{T-t}\right)-N\left(-\frac{1}{2}\sigma\sqrt{T-t}\right)\right]S.\qquad\qquad(1)$$ Now, Taylor's formula implies for small $x$ that $$N(x)=N(0)+N'(0)x+N''(0)\frac{x^2}{2}+O(x^3).\qquad\qquad\qquad\qquad(2)$$ Combining (1) and (2), we will get with some obvious cancellations that $$C(S,t)=S\left(N'(0)\sigma\sqrt{T-t}+O(\sigma^3\sqrt{(T-t)^3})\right).$$ But $$N'(0)=\frac{1}{\sqrt{2\pi}}=0.39894228...$$ so finally we have, for small $\sigma\sqrt{T-t}$, that $$C(S,t)\approx 0.4S\sigma\sqrt{T-t}.$$ The modified formula $$C(S,t)\approx 0.4Se^{-r(T-t)}\sigma\sqrt{T-t}$$

gives a slightly better approximation.

$\endgroup$
  • 3
    $\begingroup$ Good answer. I wonder if there are any approximations for options that are not at the money? My simple approach would be to assume a delta of 50% for the ATM call option and imply a price for the non ATM option as Option = ATM Price + 0.5*(Strike - Forward). Anyone got anything better? $\endgroup$ – Robert Feb 7 '12 at 12:08
  • $\begingroup$ @Robert I agree with the idea but there seems a sign mistake. Should it be Option = ATM Price + 0.5*(Forward - Strike)? So, out-of-the-money option (Strike>Forward) should be less valuable than the ATM option. $\endgroup$ – Jaehyuk Choi Aug 4 '18 at 13:02
32
$\begingroup$

This one is the best approximation I have ever seen:

If you hate computers and computer languages don't give up it's still hope! What about taking Black-Scholes in your head instead? If the option is about at-the-money-forward and it is a short time to maturity then you can use the following approximation:

call = put = StockPrice * 0.4 * volatility * Sqrt( Time )

Source: http://www.espenhaug.com/black_scholes.html

$\endgroup$
23
$\begingroup$

The Black-Scholes 'normal-vol' formula leads quickly to a similar approximation to the one described by olaker.

Click here for a paper which contains a formal derivation of the call and put prices based on a normal model (ie a brownian motion rather than a geometric brownian motion).

The formula for the call price is:

$$\text{Call} = (F-K)N(d_1) + \frac{\sigma \sqrt{T-t}}{\sqrt{2 \pi}} e^{-\frac{1}{2} d_1^2},$$

where

$$d_1 = \frac{F-K}{\sigma\sqrt{T-t}}.$$

BTW, I work in fixed income, so I always tend to write a version that is suitable for swaptions.

Options which are ATM

You can see that for $F=K$ this becomes $\text{ATMCall} = \frac{\sigma \sqrt{T-t}}{\sqrt{2 \pi}} \approx 0.4 \, \sigma \, \sqrt{T-t}$.

Options which are not ATM

I have recently discovered a generalization of this formula which works very well for strikes which are not at the money too. See my blog for a longer discussion, here are the main points.

The standard decomposition for an option is:

$$\text{Option value} = \text{Intrinsic value} + \text{Time value}.$$

In the Black-Scholes normal formula above, if you investigate the term $(F-K)N(d_1)$ in a spreadsheet, you’ll see that for small levels of volatility and maturity (try, for example, $\sigma=0.0025$, Maturity=1) it is actually quite close to $\max(0,F-K)$ – which is the intrinsic value of the call.

Consequently, the BS normal formula is almost:

$$\text{Call Price} = \text{Intrinsic Value} + \text{ATMPrice} \times e^{-\frac{1}{2} d_1^2},$$

where

$$\text{ATMPrice} = \frac{\sigma \sqrt{T-t}}{\sqrt{2 \pi}}.$$

However, if you compare this approximation to the true BS formula in a spreadsheet, you’ll see that around the strike (especially for larger values of $\sigma$) it gives too much value to the call: basically the term $\text{ATMPrice} \times e^{-\frac{1}{2}d_1^2}$, is too big when $d_1$ is non-zero and small. This is telling us that the difference between $(F-K)N(d_1)$ and $\max(0,F-K)$ gets important near the strike. As I say, have a look in a spreadsheet.

Nonetheless, this simple-but-wrong formula for the Call Price points us in the right direction: it shows that the time value of the option should be written in terms of the price of the ATM option.

Here is my solution.

I call it The Hardy Decomposition:

$$\text{Call Price} = \text{Intrinsic Value} + \text{ATMPrice}\times\text{HardyFactor}$$ where

$$\text{HardyFactor} = e^{-\frac{1}{2} d_1^2} + \frac{d_1}{0.4}N(d_1) – \max(\frac{d_1}{0.4},0).$$

So far this is just a rearrangement of the original Black-Scholes 'normal-vol' formula. The key result is that the $\text{HardyFactor}$ is well approximated by a simple expression:

$$\text{HardyFactor} \approx (1- 0.41 |d_1|) e^{-|d_1|}$$

So you can actually use the following as a pretty-good approximation for call prices:

$$\text{Call Price} = \text{Intrinsic Value} + \text{ATMPrice}\times (1- 0.41 |d_1|) e^{-|d_1|}.$$

A similar result holds for put options.

You can use this Hardy Decomposition to calculate option prices in your head - you only need to remember a few values:

Remember these few values

$\endgroup$
16
$\begingroup$

In addition to what vonjd already posted I would recommend you to look at the E.G. Haug's article - The Options Genius. Wilmott.com. You can find some aproximations of BS not only for vanilla european call and put but even for some exotics. For example:

  • chooser option: call = put = $0.4F_{0} e^{-\mu T}\sigma(\sqrt{T}-\sqrt{t})$
  • asian option: call = put = $0.23F_{0} e^{-\mu T}\sigma(\sqrt{T}+2\sqrt{t})$
  • floating strike lookback call = $0.8F_{0} e^{-\mu T}\sigma\sqrt{T} - 0.25{\sigma^2}T$
  • floating strike lookback put = $0.8F_{0} e^{-\mu T}\sigma\sqrt{T} + 0.25{\sigma^2}T$
  • forward spread: call = put = $0.4F_{1} e^{-\mu T}\sigma\sqrt{T}$
$\endgroup$
  • 1
    $\begingroup$ The wilmott link is not available now. Does anyone know a new link? $\endgroup$ – Jaehyuk Choi Aug 4 '18 at 15:15
6
$\begingroup$

The logistic distribution approximates the normal distribution function used in the Black-Scholes. The drawbacks to the normal cumulative distribution function are that it cannot be computed exactly through elementary functions, it cannot be inverted algebraically (i.e., the inverse bijection cannot be solved algebraically), and it is computationally expensive. An alternative is the logistic distribution which can be computed with elementary functions, can be algebraically inverted, and is computationally inexpsenive. Also, as a result of the inversion property, optimal boundary values of American options can be easily estimated for the logistic distribution.

The general solution to Black-Scholes is given by:

$$V_t = \Phi[d_1]S_t - \Phi[d_2]K e^{-r (T-t)} $$

where $\Phi[X]$ is the CDF of the normal distribution.

It can be shown that a fast (i.e., computationally cheap) numerical approximation is as follows:

$$V_t \approx \left(\frac{1}{2} + \frac{1}{2} \operatorname{Tanh} \left( \frac{\pi \, d_1}{2 \sqrt{3}} \right) \right) S_t - \left(\frac{1}{2} + \frac{1}{2} \operatorname{Tanh} \left( \frac{\pi \, d_2}{2 \sqrt{3}} \right) \right) K e^{-r(T-t)} = \left(1- \frac{1}{1+e^{\pi d_1 / \sqrt{3}}} \right)S_t - \left(1- \frac{1}{1+e^{\pi (d_1-\sigma \sqrt{T-t}) / \sqrt{3}}} \right) K e^{-r(T-t)}$$


The Logic

A logistic distribution $F$ -- which can be expressed as a rescaled hyperbolic tangent -- can closely approximate the normal distribution function $\Phi$. Likewise, its inverse function -- the "logit" function $F^{-1}$ -- can be rescaled to approximate the inverse normal CDF -- the "probit" function $\Phi^{-1}$.

In comparison, the logistic distribution has fatter tails (which may be desirable). Whereas the normal distribution's CDF and inverse CDF ("probit") cannot be expressed using elementary functions, closed form expressions for the logistic distribution's CDF and its inverse are facilely derived and behave like elementary algebraic functions.

The logistic distribution arises from the differential equation $\frac{d}{dx}f(x) = f(x)(1-f(x))$. Intuitively, this function is typically used to model a growth process in which the rate behaves like a bell curve. In physics, it arises as the "limit distribution of a finite-velocity damped random motion described by a telegraph process in which the random times between consecutive velocity changes have independent exponential distributions with linearly increasing parameters."

In comparison, the normal distribution arises from the following differential equation: $ \frac{d \,f(x)}{dx}=f(x)\frac{(\mu-x)}{\sigma^2}$). The normal distribution is commonly used to model diffusion processes. E.g., a Wiener processes is a stochastic process which has independent normally distributed increments with mean $\mu$ and variance $\sigma^2$. In the limit, this is a Brownian Motion.

Interestingly, the logistic distribution arises in a physical process which is analogous to Brownian motion.

Note that the CDF of the logistic distribution $F$ can be expressed using hyperbolic tangent function:

$F(x;\mu ,s)={\frac {1}{1+e^{{-{\frac {x-\mu }{s}}}}}}={\frac 12}+{\frac 12}\;\operatorname {Tanh}\!\left({\frac {x-\mu }{2s}}\right)$

Given that distribution's variance is ${\tfrac {s^{2}\pi ^{2}}{3}}$, the logistic distribution can be scaled to approximate the normal distribution by multiplying its variance $\frac{3}{\pi ^2}$. The resultant approximation will have the same first and second moments as the normal distribution, but will be fatter tailed (i.e., "platykurtotic").

Also, $\Phi$ is related to the error function (and its complement) by: $\Phi (x)={\frac {1}{2}}+{\frac {1}{2}}\operatorname {erf} \left(x/{\sqrt {2}}\right)={\frac {1}{2}}\operatorname {erfc} \left(-x/{\sqrt {2}}\right)$


Thus, for a standard normal distribution with $\mu =0$ and $\sigma =1$: $$\operatorname{erf}(\frac{x}{\sqrt{2}}) \approx \operatorname{Tanh}\left(\frac{x \, \pi}{2 \sqrt{3}} \right) \equiv \frac{e^{\frac{\pi\,x}{\sqrt{3}}}-1}{e^{\frac{\pi\,x}{\sqrt{3}}}+1} $$

$$\operatorname{erf}(x) \approx \operatorname{Tanh}\left(\frac{x \, \pi}{ \sqrt{6}} \right) \equiv \frac{e^{\pi\,x\frac{2}{\sqrt{3}}}-1}{e^{\pi\,x\frac{2}{\sqrt{3}}}+1} $$

$$\Phi \left( x \right) \approx \frac{1}{2} + \frac{1}{2} \operatorname{Tanh} \left( \frac{\pi \, x}{2 \sqrt{3}} \right) \equiv 1-\frac{1}{1+e^{\pi x \over\sqrt{3}}} $$

And easily, thus: $$x \mapsto \Phi^{-1}\left(p\right) \approx -\frac{2\sqrt{3}\operatorname{ArcTanh}\left( 1-2p \right)}{\pi}$$


Again, the chief advantage to approximating normal with the logistic distribution is that the CDF and Inverse CDF can be easily expressed using elementary functions. I have found this property to be very useful in estimating values of American options or in the general case where the terminal/optimal boundary conditions are not known.

However, it should be noted that the distributions behave differently in "the tails" which could lead to asymptotic residuals if normal is expected.

Comparison of the logit function with a scaled probit (i.e. the inverse CDF of the normal distribution), comparing {\displaystyle \operatorname {logit} (x)} \operatorname {logit} (x) vs. {\displaystyle \Phi ^{-1}(x)/{\sqrt {\frac {\pi }{8}}}} \Phi ^{-1}(x)/{\sqrt {\frac {\pi }{8}}}, which makes the slopes the same at the y-origin.

$\endgroup$
  • $\begingroup$ -1. You seem to simply approximate based on the profile of the two functions. There is neither numerical nor theoretical justification. In fact, as you can see from your plot, the approximation is not good. In fact, for large positive $x$, $1-\Phi(x)\sim e^{-\frac{x^2}2}$ rather than $e^{-ax}$ as your $1-\tanh(x)$ would have it. You have the asymptotics wrong. There are fast and theoretically concretely justified approximations. Look up the Mills ratio. $\endgroup$ – Hans Jul 17 '18 at 23:17
  • $\begingroup$ Here's a theoretical justification for modelling stock price dynamics as following telegraph processes: link.springer.com/article/10.1023%2FA%3A1009437108439 $\endgroup$ – David Addison Aug 6 at 2:54
  • $\begingroup$ Interesting paper. It is behind a paywall though. However, how does the paper address my objection which is that your answer or the linked model is not based on the Black-Scholes model which is the concern of OP's question, rather than some other model? You can raise an infinite number of other models. It is irrelevant to a question regarding one particular model, which in this case is Black-Scholes. $\endgroup$ – Hans Aug 6 at 16:16
2
$\begingroup$

As other people have said, you need to approximate the cumulative. The problem is that wherever you look you will find that to approximate it you will need to use exponentials or trigonometric functions which are also very expensive. What you can do is to build yourself a cubic spline with pre-cached values for the cumulative and calculate the value at other points x by (cubic) interpolation. That will make it much faster. Possibly you will also call this method with a series of ordered values so that you can avoid doing the binary search to locate the interval. You will have a cached index for the last interval located and look around that one.

$\endgroup$
  • 1
    $\begingroup$ There is no need for (cubic) spline. There are good numerical approximations for the error function. Look up the Mills ratio. $\endgroup$ – Hans Jul 17 '18 at 23:11
2
$\begingroup$

To answer the question by Robert about existing approximations for the Black-Scholes formula in the non-ATM case: yes there exists a simple and fast converging series representation for the Black-Scholes call, namely \begin{equation}\label{series} Call(S,K,r,\sigma,\tau) \, = \, \frac{Ke^{-r\tau}}{2} \sum\limits_{\substack{n = 0 \\ m = 1}}^{\infty} \, \frac{(-1)^n}{n!\Gamma(1+\frac{m-n}{2})} \, \left( Z^2 - [\log] \right)^{n} Z^{m-n} \end{equation} where we have denoted $[\log]:= \log\frac{S}{K} \, + \, r\tau$ and $Z:= \frac{\sigma\sqrt{\tau}}{\sqrt{2}} $. A complete proof can be found in https://arxiv.org/pdf/1710.01141.pdf . Of course it recovers the Brenner-Subrahmanyam approximation, because in the ATM case then $[\log]=0$ and the series becomes a simple power series: \begin{equation} \label{Brenner_series} Call^{(ATM)}(S,K,r,\sigma,\tau) \, = \, \frac{S}{2} \sum\limits_{\substack{n = 0 \\ m = 1}}^{\infty} \, \frac{(-1)^n}{n!\Gamma(1+\frac{m-n}{2})} \, Z^{n+m} \, = \, \frac{S}{2}\frac{1}{\Gamma(\frac{3}{2})} \, Z \, + O(Z^3) \end{equation} Recalling $\Gamma(\frac{3}{2}) = \frac{\sqrt{\pi}}{2}$, we are thus left with \begin{equation} Call^{(ATM)}(S,K,r,\sigma,\tau) \, = \, \frac{S}{\sqrt{2\pi}} \, \sigma\sqrt{\tau} \, + O\left((\sigma\sqrt{\tau})^3\right) \end{equation} which is the Brenner-Subrahmanyam approximation.

Now, the generic non-ATM case, the series is also very simple and efficient to implement. It can conveniently be represented under the form of an (n,m)-matrix: for instance below are the numerical values for the partial terms of the (n,m)-partial terms of the series where S=3800, K=4000, r=1%, $\sigma$=20% and $\tau=1Y$. enter image description here The sub-totals in the last line of the matrix show the quick convergence of the series to the Black-Scholes price of 235.514 (precision of $10^{-3}$ reached after only some iterations). Of course the convergence will be as quick as the option is close to the money, or the maturity (and therefore the $Z$-variable) is small. But in any case it will remain efficient and absolutely convergent.

I hope it helps!

$\endgroup$
1
$\begingroup$

For Java, combining @vjond answer (as the initial starting estimate for implied Volatility), with a basic Cumulative Density Function for Normal distribution (CDN), and the Black Scholes Model applied to Newton Raphson, below is a basic Implied Volatility calculation for Java :

For the NR code (C), plz refer http://finance.bi.no/~bernt/gcc_prog/recipes/recipes.pdf

For the Black Scholes price calculator (Java and other languages) and Cumulative Density function, please refer http://www.espenhaug.com/black_scholes.html


public double BlackScholes(boolean callFlag, double S, double K, double T, double r, double v)
{
    //S is underlyng spot price
    //K is the strike price
    //T is time left to expiry
    // r is the risk free rate
    // v is volatility
    //callFlag is true for Call Option, false for put

    double d1, d2;

    d1=(log(S/K)+(r+v*v/2)*T)/(v* sqrt(T));
    d2=d1-v* sqrt(T);

    if (callFlag){
        return S*CND(d1)-K* exp(-r*T)*CND(d2);
    }else{
        return K* exp(-r*T)*CND(-d2)-S*CND(-d1);
    }
}

// The cumulative normal distribution function
public double CND(double X)
{
    double L, K, w ;
    double a1 = 0.31938153, a2 = -0.356563782, a3 = 1.781477937, a4 = -1.821255978, a5 = 1.330274429;

    L = Math.abs(X);
    K = 1.0 / (1.0 + 0.2316419 * L);
    w = 1.0 - 1.0 / sqrt(2.0 * Math.PI) * exp(-L *L / 2) * (a1 * K + a2 * K *K + a3
            * Math.pow(K,3) + a4 * Math.pow(K,4) + a5 * Math.pow(K,5));

    if (X < 0.0)
    {
        w= 1.0 - w;
    }
    return w;
}


public double implVolNR(double S, double K, double r, double time, double option_price, boolean callFlag) {
    if (option_price < 0.99 * (S - K * exp(-time * r))) { // check for arbitrage violations. Option price is too low if this happens
        return 0.0;
    }
    final int MAX_ITERATIONS = 100;
    final double ACCURACY = 1.0e-2;
    double t_sqrt = Math.sqrt(time);
    double sigma = (option_price / S) / (0.398 * t_sqrt); // find initial value
    for (int i = 0; i < MAX_ITERATIONS; i++) {
        double price = BlackScholes(callFlag, S, K,time, r,sigma);



        //option_price call black scholes(S,K,r,sigma,time);
        double diff = option_price - price;
        if (Math.abs(diff) < ACCURACY) return sigma;
        double d1 = (log(S / K) + r * time) / (sigma * t_sqrt) + 0.5 * sigma * t_sqrt;
        double vega = S * t_sqrt * CND(d1);
        sigma = sigma + diff / vega;
    }
    return Double.NaN; // return error values
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.