Challenges in approximating the Black and Scholes call formula with hyperbolic tangents

In this paper we introduce the concept of standardized call function and we obtain a new approximating formula for the Black and Scholes call function through the hyperbolic tangent. This formula is useful for pricing and risk management as well as for extracting the implied volatility from quoted options. The latter is of particular importance since it indicates the risk of the underlying and it is the main component of the option's price. Further we estimate numerically the approximating error of the suggested solution and, by comparing our results in computing the implied volatility with the most common methods available in literature we discuss the challenges of this approach.


Preface
This is a reduced version of the paper without tables and figures. 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 recur to the dispersion of returns. However this measure has the problem that is computed on past performances and may have little to do with the current level of risk. Within the Black-Scholes (BS) framework [Black and Scholes, 1973], later extended by Merton [Merton, 1973], it is possible to identify a relation between the value of an asset and the option written on it. This is expressed through a link between the option's price and some factors such as time, rates, dividends and, above all, volatility [Hull, 2006]. Since then, as the usage of the BS formula has became widespread in financial markets, options are traded and priced in terms of their risk or, in other terms, of the so-called implied volatility (i.e. the actual volatility embedded in the option's price).
For the reasons above mentioned a key issue has become the analytical tractability of the said formula for pricing, hedging and inversion (needed to find the BS implied volatility).
On the other hand, given the structure of the BS formula that cannot be analytically inverted, the implied volatility can be found only through numerical approximation methods. However, in some instances, even those methods may fail for technical reasons [Orlando and Taglialatela, 2017].
In this work we show a closed form formula for approximating the BS call formula by means of the hyperbolic tangent. This allows us to determine both the value of the call for any change of the key variables and to derive the implied volatility at once for all possible combinations of underlying, strike, time, etc. To achieve the aforementioned result, we introduce the so called "standardized call function", which is a single-parameter function representing the general family of calls. The paper is organized as follows. The first section provides some background on both the topic and the literature and defines the notations used in the rest of the paper. The second and third sections propose some approximations of the call function in both the cases S = X and S = X. The fourth section is devoted to derive the approximations of the implied volatility again in the cases S = X and S = X. In the fifth section we present some results of numerical simulations and compare such results with the ones given by other authors. Finally the last section summarizes this work and draws the line for future research.
Finally Figure 1.1 below shows the overall daily distribution of the VIX [CBOE, 2017] since inception up to 2015. Throughout these years median volatility has been ∼ 18%, mode ∼ 13% and average is ∼ 20% and high values are rare.
For computing both the call price as well as the implied volatility several numerical algorithms are available [Dura and Moşneagu, 2010], [Orlando and Taglialatela, 2017]. However there is a cost and some drawbacks for those techniques which motivate, where possible, the search for closed form approximations (see for example Manaster and Koehler (1982) [Manaster and Koehler, 1982]). Notwithstanding in using those approximations one must be careful. For example Estrella (1995) [Estrella, 1995] showed that, for a plausible range of parameter values, the Taylor series for the BS formula diverges and "even when the series converges, finite approximations of very large order are generally necessary to achieve acceptable levels of accuracy". The alternative of using the exact formula with predetermined changes of the underlying, while it provides more accurate results, presents some drawbacks in terms of computations and loss of flexibility [Estrella, 1995]. Therefore "simple solutions are desirous because they have two very attractive properties. They are easy to implement, and provide very fast computational algorithm" [Hofstetter and Selbya, 2001].
Other closed-form approximations that work only in the At-The-Money case are the Pólya approximation [Pólya, 1949], [Matic et al.,0], and the logistic approximation [Pianca, 2005]. For these, it has been shown that the first is remarkably accurate for a very large range of parameters, while the second one is accurate only for extreme maturities [Pianca, 2005].
Finally, Hofstetter & Selby (2001) [Hofstetter and Selbya, 2001] obtained approximations by replacing the cndf by the logistic distribution and Li (2008) [Li, 2008] developed a closed-form method based on the rational functions.
1.2. The Black and Scholes Formula. The BS formula for deriving the price C of a European call option is described by where • S is the value of the underlying, • X = K e −rT is the present value of the strike price, • K is the strike price, • r is the interest rate, • T is the time to maturity in terms of a year, • N(x) is the cumulative distribution function of the standard normal i.e.
T is the first parameter of probability i.e. "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], T is the second parameter of probability which represents the risk-adjusted probability of exercise, • σ is the volatility. It is worth noting that, given the parameters S, X, T , (or, equivalently S, K, r, T ) the price C of the call is a function C = C (σ) of the volatility σ; the implied volatility is obtained by inverting such a function. In the following sections we shall obtain a suitable approximation C of C for all S, X, T . The implied volatility will be approximated by inverting C .
2. Approximating the call function when S = X 2.1. The standardized call function. In order to simplify the presentation we introduce a family of standardized call functions: depending on a single parameter α > 0.
The following Proposition contains the main properties of the mappings χ α .
For S, X and T fixed, the relationship between the call function C = C (σ) and the family of functions (χ α ) α>0 is contained in the following Proposition 2. Let us fix S > 0, X > 0 and T > 0, with X = S, and let us put α := 2 log(S/X) ; then we have Proof. Assume X > S; since α 2 /2 = log(X/S) = − log(S/X), and therefore X = Se α 2 /2 , we have Assume X < S; since α 2 /2 = log(S/X) = − log(X/S) and therefore S = Xe α 2 /2 , we have 2.2. Construction of the approximating functions. From Proposition 2 it follows that, in order to get a good approximation of the call functions C (σ), it is sufficient to give, for all α > 0, a good approximation χ α of χ α .
For the sake of simplicity, let us fix α > 0 so that we can omit the index α and simply denote χ and χ the mappings χ α and χ α .
By Proposition 1, we know that χ has a sigmoidal shape. This fact suggests us to look for an approximation based on the hyperbolic tangent − 1 e 2x + 1 which has a similar shape and has the advantage of having a very simple inverse function To be precise we look for an approximating function of the form where ϕ : ]0, +∞[ → R is strictly increasing and satisfies the conditions ϕ(0 + ) = −∞ and ϕ(+∞) = +∞, so that χ is strictly increasing and tends to 0 as x tends to 0 and tends to 1 as x tends to +∞.
For example we can choose ϕ of the form with c 1 > 0 and c 2 > 0 so that ϕ(x) is strictly increasing and has the desired behaviour at −∞ and +∞. Obviously, we have to choose the constants c 1 , c 2 , c 3 in such a way that χ gives the best approximation of χ; hence we impose that both χ and χ 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: Thus conditions (2.7)-(2.9) give The first equation gives c 1 − c 2 + c 3 = arctanh 2 χ(1) − 1 and we can rewrite the second and third equations as (1) ; thus (2.11) gives (2.12) and consequently (2.14) Hence c 1 , c 2 and c 3 are uniquely determined and depend only on χ(1) and χ ′ (1); from this it follows that c 1 , c 2 and c 3 depend only on α, since We have to check that c 1 and c 2 are positive. We begin with c 2 . According to (2.15) we have (2.17) 1 − 2 χ(1) = 2 e α 2 /2 N(−α) > 0 ; from this and from (2.12) it follows c 2 > 0. The proof of the positivity of c 1 is a little more involved. First of all we remark that Since the first factor is positive, in order to prove that c 1 > 0, we need to show that the second factor is positive too, that is: for any α > 0 .
2.3. The approximating Call functions. Hence, for all α > 0 we have found a good approximation to χ α , namely the function where c 1 (α), c 2 (α) and c 3 (α) are given by (2.13), (2.12) and (2.14) with χ replaced by χ α . From this and from Proposition 2, we can conclude that for all S > 0 and X > 0, with S = X, a good approximation to C (σ) is the function with α := 2 log(S/X) .

Approximating the call function when S = X
In the special case S = X, we have where erf is the error function defined by Also in this case, we consider an approximation of erf(z) which is based on the hyperbolic tangent (2.3). For example, following Ingber [Ingber, 1982], a good approximation of erf(z) may be Θ 0 (z) := tanh 2 √ π z which has the same limit at infinity and the same derivative in 0. A better approximation of erf(z) is the function which has the same Taylor expansion of order 3 in 0 of the error function, or the function Θ 2 (z) := tanh(a z + b z 3 ) , with a = 1.129324 and b = 0.100303 , obtained in [Fairclough, 2000] by an optimization procedure. Numerical tests show that Θ 1 (z) and Θ 2 (z) provide approximations of the same order of accuracy, (see §5.2 for details).
Hence we have that C (σ) can be approximated by with a = 1.129324 and b = 0.100303.

Approximation of the implied volatility
In order to find the implied volatility we should invert the call function C , i.e. we should be able to solve the equation C (σ) = C. Since C is a good approximation of C , we approximate the implied volatility by solving the equation According to (2.5) such an equation is equivalent to where [S − X] + = max(S − X, 0) is the pay-off. On the other hand, for any λ ∈ R the equation has a unique positive solution, given by Thus the implied volatility can be approximated by where Λ = 1 2 log C − [S − X] + S − C and c 1 (α), c 2 (α) and c 3 (α) are given by (2.13), (2.12) and (2.14).
Electronic copy available at: https://ssrn.com/abstract=3266556 4.2. Case S = X. In this case a first approximation of C is given by C 0 defined by (3.1) and therefore the approximating equation C 0 (σ) = C has the solution Better approximations of C are given by C 1 and C 2 defined by (3.2) and (3.3). Then the equations C 1 (σ) = C and C 2 (σ) = C are equivalent to the equations and (4.6) a σ T 8 respectively. Now it is well known that for all p > 0 the equation has a unique real solution given by Hence the unique solution to equation (4.5) is given by with p := 4 4 − π and q := 3 4 − π log S + C S − C .

5.1.
Call when S = X. By a numerical inspection we can see that the approximating function χ α (x) replicates well the standardized call function χ α (x), except from the part on the far right hand side of the inflection point where the curvature is higher (and the volatility unrealistic).
Therefore, to see better where the differences are and their magnitude we studied the comparisons between χ α (x) and the approximating function χ α (x) in (2.19) for one of the most common maturity (i.e. T = 0.25). 5.1.1. Monte Carlo analysis. Given T = 0.25, we considered the difference χ α (x) − χ α (x) for 10,000 moneyness and 500 instances of σ uniformly distributed between [0 , 1.25]. We first take the whole interval [0 , 1.25] and then we divide it into five parts (each containing 100 σ) to illustrate where the differences are higher or lower.
We noticed that the biggest errors are when σ ∈ [0.75 , 1.25] i.e. for those cases in which the volatility is extremely rare (see Fig. 1.1).

5.2.
Call when S = X. Here we start by plotting the graph of the error function and the approximating functions Θ 0 , Θ 1 and Θ 2 as defined in Section 3 then we compare numerically some significant values.
To show the approximation's error we plot the graph of the differences between the erf(z) and the functions Θ 0 , Θ 1 and Θ 2 as defined in Section 3, using a different scale for the y-axis.
Next we compute some statistical values of the differences between the error function and the functions Θ 0 , Θ 1 and Θ 2 . 5.2.1. Monte Carlo analysis. Now we consider the lattice L composed of the couples (σ, T ) with σ = 10 −4 j, with j = 1, . . . , 10 4 and T = k/12 with k = 1, . . . , 24. From it we randomly extract 10,000 samples in each interval to derive z. The statistics on the errors confirm the good quality of the approximation. 5.3. Implied volatility. As already mentioned there are several methods available in literature for deriving the implied volatility through an approximated formula. In Orlando and Taglialatela (2017) [Orlando and Taglialatela, 2017] we compared the results derived with Brenner & Subrahmanyam [Brenner and Subrahmanyam, 1998], Corrado & Miller [Corrado and Miller, 1996b], [Corrado and Miller, 1996a] and Li [Li, 2005] formulae. As we found that the latter is the most accurate, we consider Li formula as our benchmark.
5.3.1. Case when S = X. We compare the results obtained by Li formula denoted by σ L with those obtained with formula (4.3) denoted by σ. The prices of all calls have been generated with the BS 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, σ is available even when σ L is not and approximates better the implied volatility for all maturities, moneyness and level of σ except the part where it is very high (even though this occurrence is quite unlikely as mentioned in Section 1.1). Moreover the error for the σ formula is more consistent.

5.3.2.
Case when S = X. Let us now consider the ATM case. We compare the implied volatility σ L calculated with Li formula [Li, 2008], σ 1 and σ 2 as defined in (4.7) and (4.8). The prices of all calls are generated with the BS model for a given volatility ranging from 15% to 125% and maturity comprised between 0.1 and 1.5 years. We also perform some statistics on the error between the true volatility and the one estimated.

Conclusions
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. Even though calculations can be easily performed by using numerical methods we share the Li opinion [Li, 2005] that "to simplify some applications such as spreadsheet, it may be useful to have an approximation formula if that formula is reasonably simple, accurate and valid for a wide range of cases. The cost and inconvenience of iterating also motivate the search for explicit formulas".
Here a standardized call function has been introduced to represent the whole family of calls and to simplify the calculations. In addition we have shown how the approximation of the aforementioned standardized call can be performed through hyperbolic tangents instead of the usual Taylor truncation. This allows a greater accuracy for extreme values of σ, 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 superior to the ones proposed in literature so far and are valid regardless of option's moneyness. Therefore because of the higher accuracy and flexibility, this approach could replace current methods with little additional effort. Further research will address the cases when there is a marked difference between the approximation C (σ) and the BS call C (σ) in order to provide a better approximation.