What are some useful approximations to the Black-Scholes formula?

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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:

Answer 4

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}$

Answer 5

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.

Answer 6

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$. 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!

Answer 7

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.

Answer 8

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
}

Answer 9

Arguably the most useful approximation to the Black-Scholes formula would be the approximation made by any universal approximator such as a neural network, granted that it is trained on data following the BS assumptions. Buehler et al (2018) wrote an interesting paper on this concept (see https://arxiv.org/abs/1802.03042).

Simply put the idea is as follows: simulate a market environment using Black-Scholes assumptions, specify a loss function and let the universal approximator approximate the function: in this case that would be the BS formula.

The beauty (and usefulness) arises from the following: Given that you are able to simulate market dynamics that are more sophisticated than that of Black-Scholes, one can use this same universal approximator to 'find' the Black-Scholes formula in more realistic market dynamics; which is arguably a lot more useful in practice since obviously the BS assumptions don't hold in real life.