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

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?

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.

• 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? Feb 7 '12 at 12:08
• @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. Aug 4 '18 at 13:02

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 )

• Why do we require short time to maturity? The approximation fails when time to maturity is long? Also, why would call and put have the same value? Dec 4 '19 at 13:07
• @Idonknow The approximation is based on the first order Taylor Series approximation of a function of d1. When d1 is "large", the second order and higher terms of the Taylor Series become non-negligible, which makes the first order approximation "bad". Since d1 for an at-the-forward option is 0.5*sigma*sqrt(time), large time makes d1 large, which makes the approximation bad. Large volatility will have the same effect. May 21 '20 at 15:36

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:

• I don't get why, ATM, there is no the underlying spot value in your formula ? Dec 10 '19 at 13:31
• Because I am using a normal vol, not a lognormal vol. In vonjd's answer his 'volatility' term is a lognormal vol, so (S x sigma) is the corresponding normal vol. May 17 '21 at 15:54

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}$
• The wilmott link is not available now. Does anyone know a new link? Aug 4 '18 at 15:15

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.

• -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.
– Hans
Jul 17 '18 at 23:17
• Here's a theoretical justification for modelling stock price dynamics as following telegraph processes: link.springer.com/article/10.1023%2FA%3A1009437108439 Aug 6 '19 at 2:54
• 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.
– Hans
Aug 6 '19 at 16:16

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 $$\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}$$ 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: $$\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)$$ Recalling $\Gamma(\frac{3}{2}) = \frac{\sqrt{\pi}}{2}$, we are thus left with $$Call^{(ATM)}(S,K,r,\sigma,\tau) \, = \, \frac{S}{\sqrt{2\pi}} \, \sigma\sqrt{\tau} \, + O\left((\sigma\sqrt{\tau})^3\right)$$ 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!

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.

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

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
}


The premium on a put option may be estimated very easily using the following equation that I have discovered, which is fairly accurate compared to BS when strike price is on the x-axis and premium is on the y-axis:

$$2\sigma e^{-\frac{\left|k-S\right|}{4\sigma}}\ +\left|\frac{k-S}{2}\right|+\frac{k-S}{2}$$

For a call option, reflect the equation in the y axis and apply transformations to get it to the right strike price.

• Can you perhaps provide a reference or derivation where that approximation comes from? Oct 6 '21 at 15:49
• I find it odd that there’s no dependence on time to maturity. Jan 4 at 0:53
• @DaneelOlivaw if you use the standard deviation corresponding to a particular time period (e.g. when there are three days left to expiration you use the three-day standard deviation) the time to maturity will be baked into the formula. This is because standard deviation will tend to 0 as time to expiration does, and so the formula will behave in the same way as if the time were added. Jan 5 at 6:08
• @Kevin sure. The term on the left is just the laplacian distribution. The two terms on the right are the option intrinsic value. The coefficient $2\sigma$ is added to make sure that the curve is smooth. I used the laplacian just because it happened to have a similar shape to the graph of the time value on the option Jan 5 at 6:34