1 Introduction

For investors and traders, a key component in their decision making is to assess the risk they run. A common way to do so is to focus on the dispersion of returns. However, this measure has the problem that is computed on past performance and may have little to do with the current level of risk. Within the Black–Scholes framework (Black and Scholes 1973), later extended by Merton (1973), it is possible to identify a relation between the value of an asset and the option written on it. Though such a formula is subjected to strong criticism (Derman and Taleb 2005), and notwithstanding its pitfalls, the formula and its related approximations (Taylor or delta-gamma) are widely used in risk management (Estrella 1995) and “in many respects the story of the establishment of the Black–Scholes–Merton model (BSM) simply marks the emergence of contemporary financial risk management” (Millo and MacKenzie 2009). Many other models have been developed to address some of the BSM weaknesses, such as those proposed by Engle and Mustafa (1992), Heston (1993), Duan (1995) Ritchken and Trevor (1999), Christoffersen and Jacobs (2004) and Duan et al. (2006) as well as models based on neural networks (for a review see Mostafa et al. 2017) but those are beyond the scope of the present paper. This is for the simple reason that we share the same experience of most quants such as Paul Wilmott: when it comes to price and managing the risk of thousands of options in real time “better” models are not a suitable choice. Problems should be streamlined and models should be of practical use. Therefore constant volatility Black–Scholes is a common choice while tail risk is managed through worst-case scenarios. Last but not least, the BSM is, simply put, another way to ensure that “a portfolio consisting of a long position in a call and a short position in a put, valued by the traditional discounted expected value of their payoffs, must statically replicate a forward contract” (Derman and Taleb 2005). Therefore, irrespectively of its shortcomings, it is widely used as an arbitrage free tool by dealers or market-makers in options who need to stay consistent with the value of their “raw supplies” (Derman and Taleb 2005).

As the usage of the BSM formula has become widespread in financial markets (Millo and MacKenzie 2009), options are priced and traded in terms of their risk or, in other words, the so-called implied volatility (i.e. the actual volatility embedded in the option’s price). To be precise, given the value S of an underlying asset, the strike price K, the time to maturity T, the interest rate r and the volatility \(\sigma \), the BSM formula derives the price C of a European call option through the formula:

$$\begin{aligned} C := S\,N(d_1) - X\,N(d_2) , \end{aligned}$$
(1.1)

where \(X= K\,{\mathrm{e}}^{-rT}\) is the present value of the strike price, N(x) is the cumulative distribution function of the standard normal i.e.

$$\begin{aligned} N(x) := \frac{1}{\sqrt{2\pi \,}} \int _{-\infty }^x {\mathrm{e}}^{-t^2/2} \, {\mathrm{d}}t , \end{aligned}$$
(1.2)

and

$$\begin{aligned} d_1 := \frac{\log (S/X)}{\sigma \sqrt{T\,}} + \frac{\sigma }{2}\sqrt{T\,} , \qquad d_2 := \frac{\log (S/X)}{\sigma \sqrt{T\,}} - \frac{\sigma }{2}\sqrt{T\,} \end{aligned}$$

are the first and the second parameter of probability, i.e. respectively, “the factor by which the present value of contingent receipt of the stock, contingent on exercise, exceeds the current value of the stock” (Nielsen 1993) and the risk-adjusted probability of exercise.

Since the BSM formula involves the use of the cumulative distribution function of the standard normal N(x), then we have to face two kinds of problems: (a) how to compute the price of the option for given S, K, T, r and \(\sigma \); (b) given S, K, T, r, how to compute the so-called “implied volatility”, i.e. the volatility corresponding to the current price of the option; in other terms, how to invert the call function, i.e. the function \(C={\mathscr {C}}(\sigma )\) which gives the price of the option as a function of the volatility.

The first problem may be faced in several ways. An approach is to use the standard software based on the power series expansion of the cumulative normal distribution function (cndf). Alternatively N can be approximated with suitable functions, generally of rational form. In this framework, there is a variety of contributions, see for example Bowling et al. (2009), Choudhury (2014) and Eidous and Abu-Shareefa (2020) who focus on the choice of the best coefficients. Hofstetter and Selby (2001) obtained approximations by replacing the cndf with the logistic distribution, Li (2008) developed an heuristic closed-form method based on the rational functions which requires one or two steps of Newton–Raphson algorithm to improve the accuracy. Unfortunately it works only in a limited domain. Jacquier and Lorig (2015) compute the option price by a quadrature of the inverse Fourier transform and the “true implied volatility by fitting the SVI parametrization to it”. Other closed-form approximations that work only in the case \(S=X\) are the Pólya approximation (Pólya 1949; Matić et al. 2017), and the logistic approximation (Pianca 2005). A completely different point of view is due to Fan and Mancini (2009), who derive the price of an European call by a “nonparametric regression”.

Concerning the second problem, the inverse of the call function does not have an analytical representation and therefore the problem can only be approached by numerical methods (Newton–Raphson and other iterative schemes) and, in some instances, even those methods may fail for technical reasons (Orlando and Taglialatela 2017; Dura and Moşneagu 2010; Lorig et al. 2014; Liu et al. 2019). Recently, Liu et al. (2019) and Cao et al. (2019) suggested new numerical methods to reconstruct the volatility through neural networks.

A different approach to the solution of both problems is to approximate directly the call function with a suitable function, and then to use the inverse of such a function to derive the approximated implied volatility. The rationale of such an approach is that “The cost and inconvenience of iterating also motivate the search for explicit formulas. For example, traders often need to plot intra-day implied volatility in real time. In this case, the non-numerical approach such as using an explicit formula is a must” (Li 2005). In this way, practitioners can use spreadsheets to derive the price of the call as a function of the volatility and the implied volatility as a function of the price. This is of particular importance since only volatility matters in option valuation while the direction of the asset price does not. Moreover, when traders agree on the volatility they agree, also, on the option value and personal preferences are not relevant.

Generally, such closed form approximations rely on Taylor approximations or on the power series expansion of the cumulative normal distribution function (cndf) [see for example Manaster and Koehler (1982), Brenner and Subrahmanyam (1998), Bharadia et al. (1995), Chance (1996), Corrado and Miller (1996a, (1996b), Liang and Tahara (2009) and Li (2005)] [for a review of the subject see Orlando and Taglialatela (2017)]. In particular, Li proposes to approximate the call function with its third or second order Taylor polynomial and derives the implied volatility by solving a third or second degree equation. Obviously, the Taylor polynomials have local character, so the approximations are useful only in the case \(S=X\) or S close to X and only for restricted values of \(\sigma \) and C.

In this paper, we propose a different approach; to be precise, for every value of S, X and T, the call function \(C={\mathscr {C}}(\sigma )\) is a sigmoidal function with finite limits at 0 and \(+\infty \); so we propose to approximate globally \({\mathscr {C}}\) with a function which has the same behaviour at 0, at \(+\infty \) and at the inflection point of \({\mathscr {C}}\). Moreover, since we have to use the inverse of the approximating function to approximate the inverse of \({\mathscr {C}}\), we need the approximating function to have an inverse which is easily analytically representable. So we have been led to look for the approximating function by means of the hyperbolic tangent as in Orlando and Taglialatela (2020), which is the prototype of a sigmoidal function with an easy inverse.

The paper is organized as follows. In the first section, we consider the case \(S\ne X\) and we introduce the so-called “standardized call functions”, which is a single-parameter function representing the general family of calls, and we propose the desired approximations of the call functions. In the second section we propose the approximation of the call function in the case \( S = X\). The third section is devoted to deriving the approximations of the implied volatility again in the cases \(S \ne X\) and \( S = X\). In the fourth section, we present some results of numerical simulations and compare such results with the ones given by other authors. Performance of this closed form solution as well as empirical tests on market data are included. Finally, the last section summarizes the work and concludes.

2 Approximating the call function when \(S \ne X\)

2.1 The standardized call function

In order to simplify the presentation we introduce a family of standardized call functions:

$$\begin{aligned} \chi _\alpha (x) := N\biggl (\frac{\alpha }{2} \Bigl (x-\frac{1}{x}\Bigr ) \biggr ) - {\mathrm{e}}^{\alpha ^2/2} \, N\biggl (-\frac{\alpha }{2}\Bigl (x+\frac{1}{x}\Bigr )\biggr ) , \qquad x>0 , \end{aligned}$$

depending on a single parameter \(\alpha >0\).

The following Proposition contains the main properties of the mappings \(\chi _\alpha \).

Proposition 2.1

For all \(\alpha >0\) one has that:

  1. (i)

    \(\lim _{x\rightarrow 0^+} \chi _\alpha (x) = 0\) and \(\lim _{x\rightarrow +\infty } \chi _\alpha (x) = 1\).

  2. (ii)

    \(\chi _\alpha (x)\) is strictly increasing in \(]0,+\infty [\).

  3. (iii)

    \(\chi _\alpha (x)\) is strictly convex in ]0, 1] and strictly concave in \([1,+\infty [\).

Proof of (i)

It is a trivial consequence of the limits at \(\pm \infty \) of N(x).

Proof of (ii)

We have:

$$\begin{aligned} \chi _\alpha '(x)&= \frac{\alpha }{2\,\sqrt{2\pi \,}} \, \exp \biggl [-\frac{\alpha ^2}{8}\Bigl (x-\frac{1}{x}\Bigr )^2\biggr ] \Bigl (1 + \frac{1}{x^2}\Bigr ) \\&\quad + \frac{\alpha \, {\mathrm{e}}^{\alpha ^2/2}}{2\,\sqrt{2\pi \,}} \, \exp \biggl [-\frac{\alpha ^2}{8}\Bigl (x+\frac{1}{x}\Bigr )^2\biggr ] \Bigl (1 - \frac{1}{x^2}\Bigr ) \end{aligned}$$

and, since

$$\begin{aligned} \exp \biggl [-\frac{\alpha ^2}{8}\Bigl (x-\frac{1}{x}\Bigr )^2\biggr ] = {\mathrm{e}}^{\alpha ^2/2} \, \exp \biggl [-\frac{\alpha ^2}{8}\Bigl (x+\frac{1}{x}\Bigr )^2\biggr ] \end{aligned}$$

we get

$$\begin{aligned} \chi _\alpha '(x) = \frac{\alpha }{\sqrt{2\,\pi \,}} \, \exp \biggl [-\frac{\alpha ^2}{8}\Bigl (x-\frac{1}{x}\Bigr )^2\biggr ] , \end{aligned}$$
(2.1)

from which we derive that \(\chi _\alpha (x)\) is strictly increasing in \(]0,+\infty [\). \(\square \)

Proof of (iii)

Differentiating (2.1) we get

$$\begin{aligned} \chi _\alpha ''(x) = \frac{\alpha ^3}{4\,\sqrt{2\,\pi \,}} \, \exp \biggl [-\frac{\alpha ^2}{8}\Bigl (x-\frac{1}{x}\Bigr )^2\biggr ] \frac{1-x^4}{x^3} , \end{aligned}$$
(2.2)

from which we derive that \(\chi _\alpha (x)\) is strictly convex in ]0, 1] and strictly concave in \([1,+\infty [\).

For S, X and T fixed, the relationship between the call function \(C= {\mathscr {C}}(\sigma ) \) and the family of functions \((\chi _\alpha )_{\alpha >0}\) is contained in the following

Proposition 2.2

Let us fix \(S>0\), \(X>0\) and \(T>0\), with \(X \ne S\), and let us put

$$\begin{aligned} \alpha := \sqrt{2\,\bigl |\log (S/X)\bigr |\,} ; \end{aligned}$$

then we have

$$\begin{aligned} {\mathscr {C}}(\sigma ) = {\left\{ \begin{array}{ll} S\,\chi _\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) &{} {\mathrm{if }}\quad X>S , \\ S-X + X\,\chi _\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) &{} {\mathrm{if}}\quad X<S . \end{array}\right. } \end{aligned}$$

Proof

Assume \(X>S\); since \({\alpha ^2}/{2} = \log (X/S) = -\log (S/X)\), and therefore \(X=S{\mathrm{e}}^{\alpha ^2/2}\), we have

$$\begin{aligned} S \, \chi _\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr )&= S \, N\Bigl (-\frac{\alpha ^2}{2\,\sigma \,\sqrt{T\,}} + \frac{\alpha }{2}\,\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) - S\,{\mathrm{e}}^{\alpha ^2/2} \, N\Bigl (-\frac{\alpha ^2}{2\,\sigma \,\sqrt{T\,}} - \frac{\alpha }{2}\,\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) \\&= S \, N\Bigl (\frac{\log (S/X)}{\sigma \,\sqrt{T\,}} + \frac{\sigma \,\sqrt{T\,}}{2}\Bigr ) - X \, N\Bigl (\frac{\log (S/X)}{\sigma \,\sqrt{T\,}} - \frac{\sigma \,\sqrt{T\,}}{2}\Bigr ) = {\mathscr {C}}(\sigma ) . \end{aligned}$$

Assume \(X<S\); since \(\alpha ^2/2 = \log (S/X) = -\log (X/S)\) and therefore \(S= X {\mathrm{e}}^{\alpha ^2/2}\), we have

$$\begin{aligned} X \, \chi _\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr )&= X \, N\Bigl (-\frac{\alpha ^2}{2\,\sigma \,\sqrt{T\,}} + \frac{\alpha }{2}\,\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) - X\,{\mathrm{e}}^{\alpha ^2/2} \, N\Bigl (-\frac{\alpha ^2}{2\,\sigma \,\sqrt{T\,}} - \frac{\alpha }{2}\,\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) \\&= X \, N\Bigl (-\frac{\log (S/X)}{\sigma \,\sqrt{T\,}} + \frac{\sigma \,\sqrt{T\,}}{2}\Bigr ) - S \, N\Bigl (- \frac{\log (S/X)}{\sigma \,\sqrt{T\,}} - \frac{\sigma \,\sqrt{T\,}}{2}\Bigr ) . \end{aligned}$$

Now, since

$$\begin{aligned} N(x) = 1 - N(-x) , \qquad \text {for any }x\in {\mathbb {R}}, \end{aligned}$$

we get

$$\begin{aligned} X \, \chi _\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr )&= X - X \, N\Bigl (\frac{\log (S/X)}{\sigma \,\sqrt{T\,}} - \frac{\sigma \,\sqrt{T\,}}{2}\Bigr ) - S + S\, N\Bigl (\frac{\log (S/X)}{\sigma \,\sqrt{T\,}} + \frac{\sigma \,\sqrt{T\,}}{2}\Bigr ) \\&= X - S + {\mathscr {C}}(\sigma ) . \end{aligned}$$

\(\square \)

2.2 Construction of the approximating functions

From Proposition 2.2 it follows that, in order to get a good approximation of the call functions \({\mathscr {C}}(\sigma )\), it is sufficient to give, for all \(\alpha >0\), a good approximation \(\widehat{\chi }_\alpha \) of \(\chi _\alpha \).

For the sake of simplicity, let us fix \(\alpha >0\) so that we can omit the index \(\alpha \) and simply denote \(\chi \) and \(\widehat{\chi }\) the mappings \(\chi _\alpha \) and \(\widehat{\chi }_\alpha \).

By Proposition 2.1, we know that \(\chi \) has a sigmoidal shape.

This fact suggests us to look for an approximation based on the hyperbolic tangent

$$\begin{aligned} \tanh (x) = \frac{{\mathrm{e}}^x-{\mathrm{e}}^{-x}}{{\mathrm{e}}^x+{\mathrm{e}}^{-x}} = \frac{{\mathrm{e}}^{2x}-1}{{\mathrm{e}}^{2x}+1} \end{aligned}$$
(2.3)

which has a similar shape and has the advantage of having a very simple inverse function

$$\begin{aligned} {{\,\mathrm{arctanh}\,}}(x) = \frac{1}{2} \, \log \Bigl (\frac{1+x}{1-x}\Bigr ) . \end{aligned}$$
(2.4)

In our case, we look for an approximating function of the form

$$\begin{aligned} \widehat{\chi }(x) := \frac{1}{2} + \frac{1}{2}\,\tanh \bigl (\varphi (x)\bigr ) = \frac{{\mathrm{e}}^{2\,\varphi (x)}}{{\mathrm{e}}^{2\,\varphi (x)}+1} \end{aligned}$$
(2.5)

where \(\varphi :]0,+\infty [ \rightarrow {\mathbb {R}}\) is strictly increasing and satisfies the conditions \(\varphi (0^+) = -\infty \) and \(\varphi (+\infty ) = + \infty \), so that \(\widehat{\chi }\) is strictly increasing and tends to 0 as x tends to 0 and tends to 1 as x tends to \(+\infty \).

For example, we can choose \(\varphi \) of the form

$$\begin{aligned} \varphi (x) := c_1\,x-\frac{c_2}{x}+c_3 \end{aligned}$$
(2.6)

with \(c_1>0\) and \( c_2>0\) so that \(\varphi (x)\) is strictly increasing and has the desired behaviour at 0 and \(+\infty \).

Obviously, we have to choose the constants \(c_1, c_2, c_3\) in such a way that \(\widehat{\chi }\) gives the best approximation of \(\chi \); hence we impose that both \(\widehat{\chi }\) and \(\chi \) have an inflection point at \(x=1\) with the same tangent lines there, i.e. we impose that \(c_1\), \(c_2\) and \(c_3\) have to satisfy the conditions:

$$\begin{aligned} {\hat{\chi }}(1) = \chi (1) , \qquad \qquad {\hat{\chi }}'(1) = \chi '(1) , \qquad \qquad {\hat{\chi }}''(1) = \chi ''(1) = 0 . \end{aligned}$$
(2.7)

As

$$\begin{aligned} \tanh '(x) = 1-\tanh ^2(x) , \end{aligned}$$

we have

$$\begin{aligned} \widehat{\chi }'(x)&= \frac{1}{2} \, \Bigl [1-\tanh ^2\bigl (\varphi (x)\bigr )\Bigr ]\,\varphi '(x) , \\ \widehat{\chi }''(x)&= - \tanh \bigl (\varphi (x)\bigr )\Bigl [1-\tanh ^2\bigl (\varphi (x)\bigr )\Bigr ]\,[\varphi '(x)]^2 + \frac{1}{2} \, \Bigl [1-\tanh ^2\bigl (\varphi (x)\bigr )\Bigr ]\,\varphi ''(x) = \\&= \frac{1}{2} \, \Bigl [1-\tanh ^2\bigl (\varphi (x)\bigr )\Bigr ] \, \biggl [\varphi ''(x) - 2\,\tanh \bigl (\varphi (x)\bigr )\,\bigl [\varphi '(x)\bigr ]^2\biggr ] . \end{aligned}$$

Thus, conditions (2.7) give

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{1}{2} + \frac{1}{2}\,\tanh (c_1-c_2+c_3) = \chi (1) = N(0) - {\mathrm{e}}^{\alpha ^2/2} N(-\alpha ) = \frac{1}{2} - {\mathrm{e}}^{\alpha ^2/2} N(-\alpha ) \\ \frac{1}{2} \, \Bigl [1-\tanh ^2(c_1-c_2+c_3)\Bigr ]\,(c_1+c_2) = \chi '(1) = \frac{\alpha }{\sqrt{2\pi \,}} \\ -2\,c_2 - 2\,\tanh (c_1-c_2+c_3)\,(c_1+c_2)^2 = 0 . \end{array}\right. } \end{aligned}$$

If we put \(\beta = {\mathrm{e}}^{\alpha ^2/2} N(-\alpha )\), then we can rewrite such equations as

$$\begin{aligned} \tanh (c_1-c_2+c_3)&= -2\beta \, \end{aligned}$$
(2.8)
$$\begin{aligned} (1-4\beta ^2)\,(c_1+c_2)&= \frac{2\alpha }{\sqrt{2\pi \,}} , \end{aligned}$$
(2.9)
$$\begin{aligned} c_2 - 2\,\beta \,(c_1+c_2)^2&= 0 . \end{aligned}$$
(2.10)

From (2.9) we get

$$\begin{aligned} c_1+c_2 = \frac{2\alpha }{\sqrt{2\pi \,} \, (1-4\beta ^2)} ; \end{aligned}$$

thus Eq. (2.10) gives

$$\begin{aligned} c_2 = \frac{4\alpha ^2\beta }{\pi \,(1-4\beta ^2)^2} , \end{aligned}$$
(2.11)

and consequently

$$\begin{aligned} c_1 = \frac{2\alpha }{ \sqrt{2\pi \,} \, (1-4\beta ^2) } - \frac{ 4\alpha ^2\beta }{\pi (1-4\beta ^2)^2} = \frac{ \alpha \, \bigl [ \sqrt{2\pi \,} (1-4\beta ^2) - 4\alpha \beta \bigr ] }{ \pi ( 1-4\beta ^2)^2} . \end{aligned}$$
(2.12)

Finally, from (2.8) it follows

$$\begin{aligned} c_3 = {{\,\mathrm{arctanh}\,}}(-2\,\beta ) - c_1 + c_2 = \frac{1}{2} \log \Bigl (\frac{1+2\beta }{1-2\beta }\Bigr ) + \frac{\alpha \bigl [ 8\alpha \beta -\sqrt{2\pi \,}\,(1-4\beta ^2) \bigr ] }{\pi (1-4\beta ^2)^2} . \end{aligned}$$
(2.13)

Hence, \(c_1\), \(c_2\) and \(c_3\) are uniquely determined and depend only on \(\alpha \). We have only to check that \(c_1\) and \(c_2\) are positive.

The positivity of \(c_2\) follows from (2.11); the proof of the positivity of \(c_1\) is a little more involved.

In fact the sign of \(c_1\) is the same as the sign of

$$\begin{aligned} 2\pi (1-4\beta ^2) - 4\alpha \beta \sqrt{2\pi \,}&= \alpha ^2 + 2\pi - \bigl ( \alpha + 2\sqrt{2\pi \,} \beta \bigr )^2 \nonumber \\&= \bigl (\sqrt{\alpha ^2+ 2\pi } - \alpha - 2\sqrt{2\pi \,} \beta \bigr ) \cdot \bigl (\sqrt{\alpha ^2+ 2\pi } + \alpha + 2\sqrt{2\pi \,} \beta \bigr ) \end{aligned}$$
(2.14)

The second factor is positive; hence, the sign of \(c_1\) is the same as the sign of the first one, i.e.

$$\begin{aligned} \sqrt{\alpha ^2+2\pi \,} - \alpha - 2\sqrt{2\pi \,}\,{\mathrm{e}}^{\alpha ^2/2} N(-\alpha ) = {\mathrm{e}}^{\alpha ^2/2} \cdot \Bigl ( {\mathrm{e}}^{-\alpha ^2/2} \bigl [ \sqrt{\alpha ^2+2\pi \,} - \alpha \bigr ] - 2 \int _{-\infty }^{-\alpha } {\mathrm{e}}^{-t^2/2} \, {\mathrm{d}}t \Bigr ) . \end{aligned}$$
(2.15)

So, in order to prove that \(c_1\) is positive, it is sufficient to prove that

$$\begin{aligned} {\mathrm{e}}^{-\alpha ^2/2} \bigl [ \sqrt{\alpha ^2+2\pi \,} - \alpha \bigr ] - 2 \int _{-\infty }^{-\alpha } {\mathrm{e}}^{-t^2/2} \, {\mathrm{d}}t> 0 \qquad \hbox { for all}\ \alpha >0 . \end{aligned}$$
(2.16)

Such an inequality is an immediate consequence of the Komatsu–Pollak estimate (Komatsu 1955; Pollak 1956). For the reader’s convenience, we provide a proof in “Appendix”.

2.3 The approximating call functions

Hence, for all \(\alpha >0\) we have found a good approximation to \(\chi _\alpha \), namely the function

$$\begin{aligned} \widehat{\chi }_\alpha (x) = \frac{1}{2} \Bigl ( 1 + \tanh \bigl (\varphi _\alpha (x)\bigr ) \, \Bigr ) \quad \text {with} \quad \varphi _\alpha (x) := c_1(\alpha ) \, x - \frac{c_2(\alpha )}{x} + c_3(\alpha ) , \end{aligned}$$
(2.17)

where \(c_1(\alpha )\), \(c_2(\alpha )\) and \(c_3(\alpha )\) are given by (2.12), (2.11) and (2.13).

From this and from Proposition 2.2, we can conclude that for all \(S>0\) and \(X>0\), with \(S\ne X\), a good approximation to \({\mathscr {C}}(\sigma )\) is the function

$$\begin{aligned} \widehat{{\mathscr {C}}}(\sigma ) := {\left\{ \begin{array}{ll} S \, \widehat{\chi }_\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) \qquad &{} \mathrm{if}\quad X>S \\ S-X + X\,\widehat{\chi }_\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) \qquad \qquad &{} \mathrm{if}\quad X<S \end{array}\right. } \end{aligned}$$
(2.18)

with

$$\begin{aligned} \alpha := \sqrt{2\,\bigl |\log (S/X)\bigr |\,} . \end{aligned}$$

3 Approximating the call function when \(S=X\)

In the special case \(S=X\), we have \(d_1(\sigma ) = \frac{\sigma }{2}\,\sqrt{T\,}\) and \(d_2(\sigma ) = - \frac{\sigma }{2}\,\sqrt{T\,}=-d_1\). Hence, \({\mathscr {C}}(\sigma )\) reduces to

$$\begin{aligned} {\mathscr {C}}(\sigma ) = S\,N\Bigl (\frac{\sigma }{2}\,\sqrt{T\,}\Bigr ) - S\,N\Bigl (-\frac{\sigma }{2}\,\sqrt{T\,}\Bigr ) = S\Bigl [2\,N\Bigl (\frac{\sigma }{2}\,\sqrt{T\,}\Bigr ) - 1\Bigr ] , \end{aligned}$$
(3.1)

since \(N(-x) = 1-N(x)\) for all \(x\in {\mathbb {R}}\). Recalling the error function defined by

$$\begin{aligned} {{\,\mathrm{erf}\,}}(z) := \frac{2}{\sqrt{\pi \,}} \int _0^z {\mathrm{e}}^{-t^2} \ {\mathrm{d}}t , \end{aligned}$$

which is related to N(x) by the identity

$$\begin{aligned} N(x) = \frac{1}{2} \Bigl [1+{{\,\mathrm{erf}\,}}\Bigl (\frac{x}{\sqrt{2\,}}\Bigr )\Bigr ] \end{aligned}$$
(3.2)

we can also write

$$\begin{aligned} {\mathscr {C}}(\sigma ) = S\,{{\,\mathrm{erf}\,}}\biggl (\sigma \,\sqrt{\frac{T}{8}\,}\biggr ) . \end{aligned}$$
(3.3)

Hence, in order to approximate the call function, we have to approximate the function N(x) or equivalently \({{\,\mathrm{erf}\,}}(x)\). There are a lot of approximated formulae for N(x) as well as for \({{\,\mathrm{erf}\,}}(x)\), see Choudhury (2014), Eidous and Abu-Shareefa (2020) and the references therein. However, not all of these formulas are useful for our purpose, since the inverse is not easy to compute; (e.g. Eidous and Abu-Shareefa 2020) or they are similar to the hyperbolic tangent [see Page (1977) or the almost identical formula in Bowling et al. (2009)].

We propose to approximate \({{\,\mathrm{erf}\,}}(x)\) with the function

$$\begin{aligned} \Theta (x) := \tanh \biggl ( \frac{2}{\sqrt{\pi }} \, x + \frac{8-2\pi }{3\sqrt{\pi ^3}} \, x^3 \, \biggr ) , \end{aligned}$$

which has the same limit at \(+\infty \) and the same Taylor expansion of third order in 0 as the error function, and therefore we propose to approximate \({\mathscr {C}}(\sigma )\) with

$$\begin{aligned} \widehat{{\mathscr {C}}}(\sigma ) = S\,\tanh \biggl (\, \Bigl ( \sigma \,\sqrt{\frac{T}{2\pi }\,} \; \Bigr ) + \frac{4-\pi }{12} \, \Bigl (\sigma \,\sqrt{\frac{T}{2\pi }\,}\; \Bigr )^3\biggr ) . \end{aligned}$$
(3.4)

Notice that \(\Theta (x)\) is quite similar to the functions

$$\begin{aligned} \Theta _F(x)&= \tanh (a\,x + b\,x^3) ,&\quad&\text {with }&a&= 1.129324&\quad&\text {and }&b&= 0.100303 , \end{aligned}$$
(3.5)
$$\begin{aligned} \Theta _P(z)&= \tanh (a'\,x + b'\,x^3) ,&\quad&\text {with }&a'&= 1.1297&\quad&\text {and }&b'&= 0.0997 , \end{aligned}$$
(3.6)

obtained, respectively, in Fairclough (2000) and Page (1977) by an optimization procedure.

Numerical tests show that \(\Theta (x)\), \(\Theta _F(x)\) and \(\Theta _P(x)\) provide approximations of the same order of accuracy (see Sect. 5.2 for details).

4 Approximation of the implied volatility

In order to find the implied volatility we should invert the call function \({\mathscr {C}}\), i.e. we should be able to solve the equation \({\mathscr {C}}(\sigma )=C\). Since \(\widehat{{\mathscr {C}}}\) is a good approximation of \({\mathscr {C}}\), we approximate the implied volatility by solving the equation

$$\begin{aligned} \widehat{{\mathscr {C}}}(\sigma ) = C . \end{aligned}$$
(4.1)

4.1 Case \(S \ne X\)

In this case, by (2.18), equation (4.1) is equivalent to

$$\begin{aligned} \widehat{\chi }_\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) = C^* , \quad \text {where} \quad C^* := {\left\{ \begin{array}{ll} C/S &{} \mathrm{if}\quad X>S , \\ (C-S+X)/X &{} \mathrm{if}\quad X<S . \end{array}\right. } \end{aligned}$$

According to (2.5) such an equation is equivalent to

$$\begin{aligned} \varphi _\alpha \Bigl (\frac{\sigma \,\sqrt{T\,}}{\alpha }\Bigr ) = \frac{1}{2} \, \log \Bigl (\frac{C^*}{1-C^*}\Bigr ) = \frac{1}{2} \, \log \Bigl (\frac{C - [S-X]^+}{S-C}\Bigr ) , \end{aligned}$$
(4.2)

where \([S-X]^+ = \max (S-X,0)\) is the pay-off.

On the other hand, for any \(\lambda \in {\mathbb {R}}\) the equation

$$\begin{aligned} c_1\,x - \frac{c_2}{x} + c_3 = \lambda \, \qquad \text {i.e.} \qquad c_1\,x^2 - (\lambda -c_3)\,x-c_2 = 0 \end{aligned}$$

has a unique positive solution, given by

$$\begin{aligned} x = \frac{1}{2\,c_1} \Bigl [\lambda - c_3 + \sqrt{(\lambda -c_3)^2 + 4\,c_1\,c_2 \, }\Bigr ] . \end{aligned}$$

Thus, the implied volatility can be approximated by

$$\begin{aligned} \widehat{\sigma } = \frac{\alpha }{2\,c_1(\alpha ) \, \sqrt{T\,}} \left[ \Lambda - c_3(\alpha ) + \sqrt{\bigl (\Lambda - c_3(\alpha )\bigr )^2 + 4\,c_1(\alpha ) \,c_2(\alpha ) \, }\right] , \end{aligned}$$
(4.3)

where \(\Lambda = \frac{1}{2} \, \log \Bigl (\frac{C - [S-X]^+}{S-C}\Bigr )\) and \(c_1(\alpha )\), \(c_2(\alpha )\) and \(c_3(\alpha )\) are given by (2.12), (2.11) and (2.13).

4.2 Case \(S=X\)

In this case, the equation \(\widehat{{\mathscr {C}}}(\sigma ) = C\) is equivalent to the equation

$$\begin{aligned} \biggl (\sigma \,\sqrt{\frac{T}{2\pi }\,}\biggr ) + \frac{4-\pi }{12}\,\biggl (\sigma \,\sqrt{\frac{T}{2\pi }\,}\biggr )^3 = \frac{1}{2} \, \log \Bigl (\frac{S+C}{S-C}\Bigr ) . \end{aligned}$$
(4.4)

Now it is well known that for all \(p>0\) the equation

$$\begin{aligned} x^3 + 3\,p\,x = 2\,q \end{aligned}$$

has a unique real solution given by

$$\begin{aligned} x = \root 3 \of { \sqrt{ p^3 + q^2 \,} + q } \; - \; \root 3 \of { \sqrt{ p^3 + q^2 \,} - q \,} . \end{aligned}$$

Hence, the approximating value of the implied volatility, i.e. the unique solution to Eq. (4.4) is given by

$$\begin{aligned} \widehat{\sigma } := \sqrt{\frac{2\pi }{T}\,} \biggl [\root 3 \of { \sqrt{ p^3 + q^2 \,} + q } \; - \; \root 3 \of { \sqrt{ p^3 + q^2 \,} - q \,}\biggr ] , \end{aligned}$$
(4.5)

with

$$\begin{aligned} p := \frac{4}{4-\pi } \quad \text {and} \quad q := \frac{3}{4-\pi } \, \log \Bigl (\frac{S+C}{S-C}\Bigr ) . \end{aligned}$$

5 Numerical results

Volatility in the markets has been calculated by many. Mo and Wu (2007), on a sample from January 3, 1996, to August 14, 2002, of 346 weekly observations for each index, reported that the implied volatility on the S&P 500 Index (SPX), the FTSE 100 Index (FTS), and the Nikkei-225 Stock Average (NKY) ranges from 15 to 35% with S/X comprised between 80 and 120% and maturity 1 m, 3 m, 6 m and 12 m. Glasserman and Wu (2011), on a sample consisting of 2049 data points from August 9, 2001 to June 16, 2009, found that the implied volatility ranges from 5 to 43% with currency options on EURUSD, GBPUSD and USDJPY for at 10-delta put (P10d), 25-delta put (P25d), At-The-Money (ATM), 25-delta call (C25d), and 10-delta call (C10d).

Figure 1 shows the overall daily distribution of the VIX (CBOE 2020) since inception up to 2015. Throughout these years median volatility has been \(\sim 18\%\), mode \(\sim 13\%\) and average is \(\sim 20\%\) and high values are rare.

Fig. 1
figure 1

Overall daily distribution of the VIX. Distribution is heavily skewed to the left. Source Mehta (2015). Volatilities above \(\sim 30\%\) belong to the 10th decile while the remaining \(\sim 90\%\) of the VIX distribution is comprised between 0 and \(\sim 30\%\)

Full size image

5.1 Call when \(S \ne X\)

5.1.1 On the estimation error

First of all notice that for all \(\alpha > 0\) the approximation error

$$\begin{aligned} E_\alpha (x) = \chi _\alpha (x) - \widehat{\chi }_\alpha (x) \end{aligned}$$

corresponds to the moneyness \(S/X={\mathrm{e}}^{\alpha ^2/2}>1\) as well as to the moneyness \(S/X={\mathrm{e}}^{-\alpha ^2/2}<1\).

We start by showing (Fig. 2) the graph of \(E_\alpha (x)\) for \( \alpha = 0.3124\) corresponding to the moneyness \(S/X = 1.05\), as well as \(S/X=0.95\). There it can be seen that the error of the approximation \(\bigl |E_\alpha (x)\bigr |\) is greater than 0.01 only for \(x>2\), e.g. for \(T=0.25\), when \(\sigma > 1.2495\) which is an unrealistic level of the volatility.

Fig. 2
figure 2

Graph of the difference between the standardized call function \(\chi _\alpha (x)\) and its approximating function \(\widehat{\chi }_\alpha (x)\) for \(\alpha = 0.3124\) corresponding to \(S/X=1.05\), or \(S/X=0.95\)

Full size image

More generally, Tables 1 and 2 report the error \(\bigl |E_\alpha (x)\bigr |\) for some significant values of x and \(\alpha \) in absolute and percentage terms. Notice that the chosen values of \(\alpha \) correspond approximately to S/X equal to 1.03, 1.13, 1.32, 1.65, 2.18 as well as S/X equal to 0.97, 0.88, 0.75, 0.61, 0.46 respectively. These tables show that the approximation is sufficiently good, especially for the small values of x, which correspond to realistic values of sigma.

Table 1 Difference between \(\chi _\alpha (x)\) and \(\widehat{\chi }_\alpha (x)\)
Full size table
Table 2 Error between \(\chi _\alpha (x)\) and \(\widehat{\chi }_\alpha (x)\)
Full size table

We have also tried to compute the maximum value of the approximation error as a function of \(\alpha \). Indeed for every \(\alpha >0\) the function \(\bigl |E_\alpha (x)\bigr |\) is smooth in \(]0,+\infty [\), and tends to 0 as x tends to 0 and to \(+\infty \). Thus, the maximum value of \(\bigl |E_\alpha (x)\bigr |\) in the interval \(]0,+\infty [\) exists and can be located by finding the zeros of its derivative,

$$\begin{aligned} E_\alpha ^{'}(x) = \frac{\alpha }{\sqrt{2\,\pi \,}} \, \exp \biggl [-\frac{\alpha ^2}{8}\Bigl (x-\frac{1}{x}\Bigr )^2\biggr ] - \frac{2\,{\mathrm{e}}^{2\,\varphi _\alpha (x)}\,\varphi _\alpha '(x)}{({\mathrm{e}}^{2\,\varphi _\alpha (x)}+1)^2} , \end{aligned}$$

where \(\varphi _\alpha (x):= c_1(\alpha )\,x-\frac{c_2(\alpha )}{x}+c_3(\alpha )\) and \(c_1(\alpha )\), \(c_2(\alpha )\) and \(c_3(\alpha )\) are given by (2.12), (2.11) and (2.13). Unfortunately, finding the zeros of \(E_\alpha '\) is an almost impossible task. The same applies to finding a function \(\Phi =\Phi (\alpha )\) such that for every \(\alpha >0\) one has \(\bigl |E_\alpha (x)\bigr | \le \Phi (\alpha )\) for all x.

For the said reason, for a given \(\alpha >0\), in order to compute

$$\begin{aligned} G_{\alpha }=\max \Bigl \{ \ \bigl |E_\alpha (x)\bigr | \ \Bigm | \ x>0 \ \Bigr \} , \end{aligned}$$

we resort to numerical methods, namely to the Nelder–Mead simplex algorithm (Nelder and mead 1965; Gao and Han 2012), which, even though it is slower than gradient-based algorithms, works well when the call function to approximate is steep or flat [for issues related to those cases see Orlando and Taglialatela (2017)]. Table 3 reports the maximum value of the approximation error for the same values of \(\alpha \) considered in Tables 1 and 2.

Table 3 Maximum error between \(\chi _\alpha (x)\) and \(\widehat{\chi }_\alpha (x)\). Nelder–Mead method
Full size table

5.1.2 Monte Carlo analysis on the estimation error

Notice that for all \(\alpha > 0\), \(T > 0\), \(\sigma > 0\) the approximation error

$$\begin{aligned} \Bigl | E_\alpha \bigl ( \sigma \sqrt{T\,}/\alpha \bigr ) \Bigr | = \Bigl | \chi _\alpha \bigl ( \sigma \sqrt{T\,} /\alpha \bigr ) - \widehat{\chi }_\alpha \bigl ( \sigma \sqrt{T\,} /\alpha \bigr ) \Bigr | \end{aligned}$$

represents the approximation error between the call \( {{\mathscr {C}}} (\sigma )\) and its approximation \(\widehat{{\mathscr {C}}}(\sigma )\) for the case when \(X=1\) and \(S= {\mathrm{e}}^{\alpha ^2/2} >1 = X\) as well for the case \(S=1\) and \(X= {\mathrm{e}}^{\alpha ^2/2} >1 = S\).

The following Table 4 shows that such an approximation error is very good for the most traded combination of T, S/X and \(\sigma \).

We have also performed a Monte Carlo analysis on the error

$$\begin{aligned} \Bigl | E_\alpha \bigl (\sigma \sqrt{T}/\alpha \bigr )\Bigr | = {\mathscr {C}}(\sigma ) -\widehat{{\mathscr {C}}}(\sigma ) \end{aligned}$$

for 10, 000 values of \(\alpha \) and 500 instances of \(\sigma \), uniformly distributed in the interval ]0, 1.25]. We first take the whole interval ]0, 1.25] and then we divide it into five parts (each containing 100 \(\sigma \)) to illustrate where the differences are higher or lower.

In Tables 5, 6 and 7, we show such an analysis, respectively, for the cases where time is equal to 1 month, one quarter and one semester (respectively).

Table 4 Statistics on the difference between \(\chi _\alpha (\sigma \sqrt{T}/\alpha )\) and \(\widehat{\chi }_\alpha (\sigma \sqrt{T}/\alpha )\)
Full size table
Table 5 Statistics on the difference between \({\mathscr {C}}(\sigma ) \) and \(\widehat{{\mathscr {C}}}(\sigma )\), case \(T=1/12\)
Full size table
Table 6 Statistics on the difference between \({\mathscr {C}}(\sigma ) \) and \(\widehat{{\mathscr {C}}}(\sigma )\), case \(T=0.25\)
Full size table
Table 7 Statistics on the difference between \({\mathscr {C}}(\sigma ) \) and \(\widehat{{\mathscr {C}}}(\sigma )\), case \(T=0.5\)
Full size table

Notice that the biggest errors are when \(\sigma \in [0.75,1.25]\) i.e. for those cases in which the volatility is extremely rare (see Fig. 1).

5.2 Call when \(S = X\)

To show the approximation’s error, in Fig. 3 we plot the graph of the difference between the functions \({{\,\mathrm{erf}\,}}(x)\) and \(\Theta (x)\) as defined in Sect. 3.

Fig. 3
figure 3

The graph of \({{\,\mathrm{erf}\,}}(x)-\Theta (x)\)

Full size image

As previously mentioned, there is a large number of approximations for the standard normal distribution or, equivalently, for the error function. In Table 8, we compare our approximations with Bowling et al. (2009) [which is equivalent to Page (1977)], and Eidous and Abu-Shareefa (2020). As shown, the approximations are very close and \(\Theta \) has the smallest RMSE.

Table 8 Comparison between percent change of the error due by approximating \({{\,\mathrm{erf}\,}}\) with the function \(\Theta \) and with the Fairclough, Page and Eidous et al. formulas for some values of x
Full size table

5.2.1 Monte Carlo analysis on the estimation error

Now we consider the lattice \({\mathfrak {L}}\) composed of the couples \((\sigma , T)\) with \(\sigma \in [0, 1.25]\) and \(T =k/12 \) with \(k=1,\dotsc , 24\). From it, we randomly extract 10,000 samples. The statistics on the errors reported in Table 9 confirm the good quality of the approximation.

Table 9 Statistics on the difference \({{\,\mathrm{erf}\,}}(\sigma \sqrt{T/8}) - \Theta (\sigma \sqrt{T/8})\), 10,000 simulations
Full size table
Table 10 Comparison of S. Li and (4.3) estimated implied volatilities for Out-The-Money calls
Full size table

5.3 Implied volatility

In Orlando and Taglialatela (2017), we compared the results derived with Brenner and Subrahmanyam (1998), Corrado and Miller (1996a, (1996b) and Li (2005) formulae. As we found that the latter is the most accurate, we consider S. Li formula as our benchmark. Notice that, for reader’s convenience, ranges are similar to the ones published by AMC (see Li 2005). This because we wanted to make easier for the reader a comparison with one of the most precise approximations available in the literature.

5.3.1 Case when \(S \ne X\)

In Table 10, we compare the results obtained by S. Li formula denoted by \(\widehat{\sigma }_{L}\) with those obtained with formula (4.3) denoted by \(\widehat{\sigma }\). The prices of all calls have been generated with the BSM model and, then, the implied volatility has been derived by using the inversion formulae. Each column provides the results of the said formulae for maturities T from 0.1 to 1.5 versus the true volatility.

As shown, \(\widehat{\sigma }\) is available even when \(\widehat{\sigma }_{L}\) is not and approximates better the implied volatility for all maturities, moneyness and level of \( \sigma \) except the part where it is very high (but this occurrence is very unlikely as shown in Fig. 1). Moreover the error for the \(\widehat{\sigma }\) formula is more consistent.

When S is close to X, the proposed approximating function is structurally less precise than the S. Li formula. This because his formula is derived from the Taylor polynomial with center \(S/X =1\). Nevertheless, in Table 11 we report some results for the case \(S/X = 1.05\) and low levels of \(\sigma \), where our formula produces a result close to the reality for the short term traded options (i.e. the most traded), while the benchmark fails to provide any value.

Table 11 Comparison of S. Li and (4.3) estimated implied volatilities for Out-The-Money calls
Full size table

5.3.2 Case when \(S = X\)

Let us now consider the ATM case. In Table 12, we compare the implied volatility \(\widehat{\sigma }\) as defined in (4.5) with \(\widehat{\sigma }_{L}\) calculated with S. Li formula (Li 2005). Notice that precision of our formula is of \(10^{-6}\).

Table 12 Comparison of Li and (4.5) estimated implied volatilities for At-The-Money calls
Full size table

Furthermore, in Table 13, we show a Monte Carlo simulation where the prices of all calls are generated with the BSM model for a given volatility ranging from 15 to 125% and maturity comprised between 0.1 and 1.5 years. Displayed on the left of Table 13 there are some statistics on the error between the true volatility and the estimated one.

Table 13 Statistics with respect to T on the estimation error for At-The-Money calls
Full size table

5.3.3 Empirical results on SPDR S&P 500 ETF TRUST

In Table 14 we show how Eq. (4.3) compares versus (Li 2005) and the widely used Brenner and Subrahmanyam (1998) (\(\widehat{\sigma _B}\)) formula. As shown, S. Li formula often does not provide any value, whereas Brenner and Subrahmanyam and our formula give an estimate in all cases. Unfortunately, Brenner and Subrahmanyam formula is precise only for ATM options which may exist only at the inception.

Table 14 Comparison between the SPX true volatility and the implied volatilities computed by means of Brenner et al., Li (2005) and (4.3) formulas
Full size table

5.4 Computational performance of the algorithm: hybrid method versus closed form solution

With regard to computational performance, we recall that the scope of this work is limited to analytical approximations. Among the reasons of this choice is the practical need to obtain a fast result. In order to asses that, we take as a benchmark a hybrid numeric algorithm able to quickly compute the implied volatility in a large number of cases such as the one proposed in Orlando and Taglialatela (2017). While we do not want to discuss all potential alternatives in this instance, we limit to mention that there are other solutions such as the built-in blsimpv available in Matlab (2020) or the heuristic hybrid method proposed by Li (2008) as implemented by Whirdy (2020).

The idea in both papers Li (2008) and Orlando and Taglialatela (2017) is to derive numerically the implied volatility from a suitable starting point with the following characteristics: as close as possible to the true value and able to “guide” a gradient based algorithm to not diverge. Apart from the said similarity, the two approaches are quite different in their implementation. In fact, the first one calculates Black–Scholes implied volatility by using M. Li’s rational function approximations for the initial estimate, followed by 3rd-order Householder’s root finder (i.e. using Vega, Vomma and Ultima) (Li 2008; Whirdy 2020). The second adopts the Newton–Raphson algorithm by taking as a starting point the inflection of the call function (Orlando and Taglialatela 2017). The problem with the first approach is that M. Li formula (Li 2008) works only in a very strict domain while the second approach (Orlando and Taglialatela 2017) works almost every time (unless the objective function is flat or vertical). Table 15, for a given true volatility of 10%, displays the implied volatility obtained with the two above mentioned approaches. As shown, once more, the Li (2005) formula in some instances does not provide a solution and the Whirdy algorithm (Li 2008; Whirdy 2020) shoots off.

Table 15 Comparison of implied volatilities between numerical algorithms approximation and analytical formulas (true volatility = 10%)
Full size table

As the chosen parameters are not uncommon to find in the current market environment, we discard (Li 2008; Whirdy 2020) in favour of Orlando and Taglialatela (2017).

Having demonstrated that the best available hybrid solution is Orlando and Taglialatela (2017), in this work, we compare the said algorithm with the proposed analytic formula for ATM, ITM and OTM options. In Table 16 we show the performance of the Newton–Raphson method as described in Orlando and Taglialatela (2017) versus \(\widehat{\sigma }\) on a test run over 10,000 options. The closed form solution, executed in VBA, is quite fast in finding the 10,000 volatilities for the corresponding ATM, ITM and OTM options and, as expected, it is faster than its corresponding numerical algorithm.

Table 16 Performance of the algorithm in seconds
Full size table

6 Stress tests

Stress testing is a key activity for managing the risk. Risk is represented by \(\sigma \) which is the only variable that may jump suddenly. By sampling the approximating functions in different segments of \(\sigma \), stress testing boils down to taking one of those intervals. Tables in Sects. 5.1.2, 5.2.1 and 5.3.1 list the possible outcomes for both: volatility and strike. Alternative formulas (Corrado and Miller 1996a, b; Li 2005, 2008) for some values cannot provide any solution.

7 Conclusions

In regard to the BSM formula “in many respects the story of the establishment of the Black–Scholes–Merton model simply marks the emergence of contemporary financial risk management” (Millo and MacKenzie 2009). Despite the many pitfalls of the BSM formula, there are several reasons that keep it alive. For example “the BSMs biggest strength is the possibility of estimating market volatility of an underlying asset generally as a function of price and time without, for example, direct reference to expected yield, risk aversion measures or utility functions. The second greatest strength aspect is its self-replicating strategy or hedging: explicit trading strategy in underlying assets and risk-less bonds whose terminal payoff, which equals the payoff of a derivative security at maturity. In other words, theoretically an investor can continuously buy and sell derivatives in the strategy and never incur loss. It is also simple and mathematically tractable as compared to some of its more recent variations” (Yalincak 2012). Therefore, even though in many respects BSM is outdated and has a number of shortcomings, it is still widely used especially for real time computations and for extracting the implied volatility.

In this work, we have recalled the importance of calculating the value of the call for pricing as well as for inferring the implied volatility. To be precise a standardized call function has been introduced to represent the whole family of calls and to simplify the calculations. Then we have shown how the approximation of the aforementioned standardized call can be performed through the hyperbolic tangent instead of the usual Taylor truncation. This allows a greater accuracy for extreme values of \(\sigma \), which makes this approach particularly suitable for stress testing and hedging purposes.

Finally, we have derived some explicit formulae for approximating the implied volatility that seem to be better than the ones proposed in the literature so far and are valid regardless of an option’s moneyness. Therefore because of the higher accuracy and flexibility, this approach could replace current methods with little additional effort.