A model-free backward and forward nonlinear PDEs for implied volatility

TL;DR

This paper introduces new nonlinear PDEs for modeling implied volatility of convex payoff contingent claims, providing a model-free approach with numerical solutions.

Contribution

It derives backward and forward nonlinear PDEs for implied volatility applicable to a broad class of contingent claims, including convex payoffs, without relying on specific models.

Findings

Derived PDEs for implied volatility applicable to convex payoffs

Proposed iterative finite-difference numerical method for solving the PDEs

Discussed initial and boundary conditions for practical implementation

Abstract

We derive a backward and forward nonlinear PDEs that govern the implied volatility of a contingent claim whenever the latter is well-defined. This would include at least any contingent claim written on a positive stock price whose payoff at a possibly random time is convex. We also discuss suitable initial and boundary conditions for those PDEs. Finally, we demonstrate how to solve them numerically by using an iterative finite-difference approach.

Tables

Table 1. Table 1: Parameters of the test.

$S$	$r$	$q$	$T_{max}$	$K_{\max }$	Call/Put	$N_u$	$N_v$	$R_c$
100.00	0.05	0.01	1.00	400.00	Call	200	100	20

Table 2. Table 2: Convergence of the Picard iterations for various steps in time T 𝑇 T . The time step Δ ​ T Δ 𝑇 \Delta T = 0.01.

Step in time	Iterations per step	$\epsilon$	$t$ elapsed, secs
1	20	1.5e-4	16
2	3	2.3e-6	2.7
3	3	1.9e-6	2.8
5	3	1.4e-6	2.8
10	3	8.3e-7	2.8
20	3	3.8e-7	2.8

Equations

$$\frac{\partial BS(t,S,\sigma )}{\partial t}+(r-q)S\frac{\partial BS(t,S,\sigma )}{\partial S}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}BS(t,S,\sigma )}{\partial S^2}=rBS(t,S,\sigma ),$$

$$BS(S,\Sigma (S,t,Z),t)=Z.$$

$$\Bigg{\{}\frac{\partial{BS(t,S,\sigma)}}{\partial t}+(r-q)S\frac{\partial{BS(t,S,\sigma)}}{\partial S}+\frac{1}{2}\sigma^{2}S^{2}\frac{\partial^{2}{BS(t,S,\sigma)}}{\partial S^{2}}-r{BS(t,S,\sigma)}\Bigg{\}}\biggr{|}_{\sigma=\Sigma(S,t,Z)}=0.$$

$$B{S}_t(S,\Sigma (S,t,Z),t)+B{S}_{\sigma }(S,\Sigma (S,t,Z),t){\Sigma }_t(S,t,Z)=0.$$

$$B{S}_S(S,\Sigma (S,t,Z),t)+B{S}_{\sigma }(S,\Sigma (S,t,Z),t){\Sigma }_S(S,t,Z)=0.$$

$$B{S}_{SS}(S,\Sigma ,t)+2B{S}_{S\sigma }(S,\Sigma ,t){\Sigma }_S+B{S}_{\sigma \sigma }(S,\Sigma ,t){\Sigma }_S^2+B{S}_{\sigma }(S,\Sigma ,t){\Sigma }_{SS}=0.$$

$$-B{S}_{\sigma }(S,\Sigma ,t){\Sigma }_t$$

$$+(r-q)SB{S}_S(S,\Sigma ,t)=rBS(S,\Sigma ,t).$$

$$B{S}_{\sigma }(S,\Sigma ,t){\Sigma }_Z=1.$$

$$B{S}_{S\sigma }(S,\Sigma ,t){\Sigma }_Z+B{S}_{\sigma \sigma }(S,\Sigma ,t){\Sigma }_S{\Sigma }_Z+B{S}_{\sigma }(S,\Sigma ,t).{\Sigma }_{SZ}=0.$$

$$B{S}_{\sigma \sigma }(S,\Sigma ,t){\Sigma }_Z^2+B{S}_{\sigma }(S,\Sigma ,t){\Sigma }_{ZZ}=0.$$

$$B{S}_{\sigma \sigma }(S,\Sigma ,t)=-B{S}_{\sigma }(S,\Sigma ,t)\frac{{\Sigma }_{ZZ}}{{\Sigma }_Z^2}.$$

$$B{S}_{S\sigma }(S,\Sigma ,t)=B{S}_{\sigma }(S,\Sigma ,t)\frac{{\Sigma }_S}{{\Sigma }_Z^2}{\Sigma }_{ZZ}-B{S}_{\sigma }(S,\Sigma ,t)\frac{{\Sigma }_{SZ}}{{\Sigma }_Z}.$$

$$0$$

$$-(r-q)SB{S}_S(S,\Sigma ,t)+rZ.$$

$$0$$

$${\Sigma }_t$$

$$B{S}_S(S,\Sigma ,t)=-\frac{{\Sigma }_S}{{\Sigma }_Z}.$$

$${\Sigma }_t+\frac{1}{2}{\Sigma }^2{S}^2\left[B{S}_S^2(S,\Sigma ,t){\Sigma }_{ZZ}+2B{S}_S(S,\Sigma ,t){\Sigma }_{SZ}+{\Sigma }_{SS}\right]+(r-q)S{\Sigma }_S+rZ{\Sigma }_Z=0.$$

$$B{S}_S(S,\Sigma ,t)=N(d_1(S,\Sigma ,t)),$$

$$d_1(S,\Sigma ,t)\equiv \frac{\ln \frac{S}{K}+\frac{{\Sigma }^2}{2}T}{\Sigma \sqrt{T}},$$

$${\Sigma }_t+\frac{1}{2}{\Sigma }^2{S}^2\left[{\mathcal{N}}^2(d_1(S,\Sigma ,t)){\Sigma }_{ZZ}+2\mathcal{N}(d_1(S,\Sigma ,t)){\Sigma }_{SZ}^2+{\Sigma }_{SS}\right]+(r-q)S{\Sigma }_S+rZ{\Sigma }_Z=0,$$

$$BS(K,\Sigma (K,T,Z),T)=Z,$$

$$\frac{\partial BS(K,T,\sigma )}{\partial T}$$

$$0$$

$$+(r-q)KB{S}_{\sigma }(K,\Sigma ,T){\Sigma }_K-qBS(K,\Sigma ,T).$$

$$1$$

$$0$$

Full text

A model-free backward and forward nonlinear PDEs for implied volatility

Peter Carr1

  

Andrey Itkin1 and Sasha Stoikov2

1 Tandon School of Engineering, New York University, 12 Metro Tech Center, 26th floor, Brooklyn NY 11201, USA

2 Cornell University, School of Operations Research and Information Engineering, 2W Loop Road, New York, NY 10044, USA

**We derive a backward and forward nonlinear PDEs that govern the implied volatility of a contingent claim whenever the latter is well-defined. This would include at least any contingent claim written on a positive stock price whose payoff at a possibly random time is convex. We also discuss suitable initial and boundary conditions for those PDEs. Finally, we demonstrate how to solve them numerically by using an iterative finite-difference approach. **

Introduction

As this is well-known, given an option contract written, e.g., on some stock, the implied volatility is derived from an option price and shows what the market implies about the underlying stock volatility in the future. For instance, the implied volatility is one of six inputs used in a simple option pricing (Black-Scholes) model, but is the only one that is not directly observable in the market. The standard way to determine it by knowing the market price of the contract and the other five parameters, is solving for the implied volatility by equating the model and market prices of the option contract.

There exist various reasons why traders prefer considering option positions in term of the implied volatility, rather than the option price itself, see e.g., [Natenberg, 1994]. However, in this paper our goal is not to discuss the importance of this concept. Instead we focus on the way how the implied volatility is computed, and namely - on performance of this process.

As was mentioned, the traditional method of computing the Black-Scholes implied volatility is relatively simple. Perhaps, the main problem with such an approach occurs if one needs to simultaneously compute the implied volatility for a wide set of the model (contract) parameters. Then, the iterative root solver can converge slowly because i) there is no a universal good initial guess to kick off the iterations, and ii) for some regions of the model parameters the solution either doesn’t exist at all, or is hard to be found as the sensitivity of the model price to changes in the implied volatility is very low.

Therefore, a surprising idea would be to replace an algebraic equation Eq.(2) used to determine the implied volatility with another object, for instance, a partial differential equation (PDE). Here the surprise means that the PDE is a more complicated object than the algebraic equation, and solving the PDE (especially nonlinear) requires some special methods usually much more complicated than a simple root solver. However, in contrast to the algebraic equation that should be solved independently at every given point e.g., in the space of option strikes and maturities, the PDE approach allows finding the solution simultaneously in many such the points by using a series of marching sweeps in time. Also, solving a linear PDEs doesn’t require iterations, and hence, there is no need for the good initial guess.

In this paper we demonstrate how such a PDE could be derived by using just a general definition of the implied volatility. Moreover, we further extend this idea by deriving several types of these PDEs. For instance, we obtain quasilinear backward and forward parabolic PDEs, and also nonlinear hyperbolic PDEs of the first order. As these PDEs are new and yet have not been investigated in the literature, we need to discuss suitable initial and boundary conditions that provide the solution of these PDE to exist and be unique. We also develop an iterative numerical method to solve the PDEs by using a finite-difference approach. The method is of second order of approximation in both space and time, is unconditionally stable and preserves positivity of the solution. Using this method we compute the PDE implied volaitlity and find that our intuition behind the main idea of the paper is correct. In other words, performance of the finite-difference solver exceeds that of the traditional approach by factor of forty. However, this result is subject to some details which are highlighted in the last Section.

Despite in this paper the PDE approach is considered as an alternative to the traditional method of getting the implied volatility by using the root solver, the latter in turn was seriously improved within the last few years. First, in [Jackel, 2015] modified Newton iterations were proposed to solve Eq.(2). When doing so the entire input domain is decomposed into four areas. Rational approximations are used to provide initial guesses and reduce the number of iterations. Also this method provides an accuracy close to the machine precision. However, it breaks possible vectorization (but we don’t consider this topic in this paper anyway).

In [Glau et al., 2018] the approach of [Jackel, 2015] is improved by using a bivariate Chebyshev interpolation of the implied volaitlity. This algorithm allows vectorization and consists of two steps. At the offline phase polynomial weights of a low-rank Chebyshev interpolation of the implied volatility surface are computed and stored for four different input domains, using the algorithm of [Jackel, 2015]. This step is only performed once. During the online phase the input data is split into four domains, and the Chebyshev interpolation is applied to each domain by choosing the precomputed nodes from the offline step.

From this prospective, this problem can also be solved by building an artificial neural network (ANN) and training it on a given data set which will include all parameters of the contract (the spot price, the time to maturity, the strike, option price, etc.) and output the Black-Scholes implied volatility. The method of [Jackel, 2015] can be used when training this ANN as the golden "pricer". As such, the approach proposed in this paper becomes, apparently, less important since the ANN training could be done offline, while performance of the "pricer" is less important, rather than the accuracy. This is completely opposite to the proposed PDE approach. However, using the ANN still requires solving various problems, such as no-arbitrage, the existence of derivatives, etc., see [Itkin, 2019]. Also, the advantage of using the ANN approach vs traditional methods, is, perhaps, a property inherent to many financial models, and not only to the model we consider in this paper.

Aside of practical importance, the existence of the model-free quasilinear parabolic PDE (or the other hyperbolic PDE) governing the implied volatility is an interesting fact, which yet has not been discovered in the literature. Therefore, the novelty of this paper consists in the derived equations together with the proposed numerical method of their solution.

The rest of the paper is organized as follows. In Section 1 we derive a backward quasilinear PDE for the implied volatility. In Section 2 a similar approch is used to derive an analogous forward PDE. Section 3 discusses how one can set the appropriate initial and boundary conditions for these PDEs. In Section 4 we also derive some other PDEs for the implied volatility that have a different type (nonlinear hyperbolic rather than quasilinear parabolic equations). Then we construct a numerical method to solve these equations which is described in detail in Section 5 and two Appendices. The results of our numerical experiments and comparison of two methods: i) the traditional method of finding the implied volatility by using a root solver with the Black-Scholes equation, and ii) the proposed method of solving the derived PDEs numerically, are provided in Section 6. In Section 7 we discussed the results obtained and conclude.

1 Backward PDE

Consider a class of contingent claims which have a single payoff and which are written on the spot price path of a single underlying risky asset. The Black-Scholes implied volatility, [Hull, 1997], can be defined for any such contingent claim, as long as the claim’s value in the Black-Scholes model is monotonic in the volatility of the underlying asset. Standard examples include European options and American options that are optimally held alive. If the claim’s value is decreasing in volatility, then the value of a short position in that claim is increasing in volatility. Hence, there is no loss in generality if implied volatility is only defined for short or long positions in claims whose value increases with volatility. We refer to any such claim as an option. This section derives a new nonlinear partial differential equation (PDE) that relates the Black-Scholes implied volatility to the option price, to the underlying stock price, and to calendar time.

The Black-Scholes implied volatility is a useful measure, as it is a market practice instead of quoting the option premium in the relevant currency, the options are quoted in terms of the Black-Scholes implied volatility. Over the years, option traders have developed an intuition in this quantity. However, it can be further generalized by using a similar concept, but replacing the Black-Scholes framework with another one. For instance, in [Corcuera et al., 2009] this is done under a Lévy framework, and, therefore, based on distributions that match more closely historical returns. In this paper we don’t consider these generalizations, and are concentrated only on the Black-Scholes implied volatility.

Consider a fixed time horizon $t\in[0,T]$ where $T>0$ is the time to maturity. We assume that the underlying stock price $S$ is strictly positive over this horizon. Let $BS(S,K,\sigma,r,q,t)$ denote the function relating the Black-Scholes model value $BS(S,K,\sigma,r,q,t)\in\mathbb{R}$ of a contingent claim to the underlying stock price level $S>0$, the strike $K>0$, the (constant) volatility $\sigma>0$, calendar time $t\in[0,T]$, the constant interest rate $r$ and continuous dividend $q$. Let subscripts denote partial derivatives. We assume that the function $BS(S,K,\sigma,r,q,t)$ is governed by the Black-Scholes backward PDE. To make the notation easier, further in the expression $BS(S,K,\sigma,r,q,t)$ we will drop all variables not relevant to our particular context. Therefore, instead, in this Section we denote the function value as $BS(S,\sigma,t)$. With this notation, the corresponding Black-Scholes PDE reads, [Hull, 1997]

[TABLE]

defined on some open region ${\cal R}$ of space and time. For example, for a European Call or Put, the region ${\cal R}$ would be the Cartesian product ${\cal R}_{E}\equiv[0,\infty)\times[0,T)$. As a second example, for a down-and-out Call with a lower barrier $L\in(0,S_{0})$, the region ${\cal R}$ would be the Cartesian product ${\cal R}_{L}\equiv(L,\infty)\times[0,T)$.

We furthermore assume that the Vega of the claim is always strictly positive, i.e. $BS_{\sigma}(S,\sigma,t)>0$ on ${\cal R}$. This clearly holds for a European Call or Put on ${\cal R}_{E}$, but it can also be shown that it holds for a down-and-out Call on ${\cal R}_{L}$. We also assume that one can directly observe the market price $Z(T,K)$ of such a claim, and that this market price lies in some arbitrage-free open interval ${\cal A}\subset\mathbb{R}$. For a European Call struck at $K$, ${\cal A}_{C}\equiv((Q_{T}S-D_{T}K)^{+},SQ_{T})$, where $D_{T}=e^{-rT}$ is the discount factor, and $Q_{T}=e^{-qT}$ is the price appreciation factor. For a European Put struck at $K$, ${\cal A}_{P}\equiv((D_{T}K-Q_{T}S)^{+},D_{T}K)$. For a down-and-out Call struck at $K$, ${\cal A}_{D}\equiv((D_{T}K-Q_{T}S)^{+},Q_{T}(S-L))$, [Hull, 1997].

For any such claim, we can implicitly define a function $\Sigma(S,t,K,T,r,q,Z)$ that relates the implied volatility $\Sigma$ of the contingent claim to $(S,t)\in{\cal R}$ and to the market price of the contingent claim $Z\in{\cal A}$. Again, for easy of notation we drop $K,T,r,q$ from the list of independent variables, and then the definition of $\Sigma(S,t,Z)$ reads

[TABLE]

Let ${\cal D}$ denote the domain for this function, which is defined as the Cartesian product of ${\cal R}$ and ${\cal A}$. Since the PDE Eq.(1) holds for any level of $\sigma>0$, it holds in particular if we set $\sigma=\Sigma(S,t,Z)$:

[TABLE]

We now show that Eq.(2) and Eq.(3) can be used to generate a nonlinear PDE governing $\Sigma(S,t,Z)$ on ${\cal D}$.

First, differentiate Eq.(2) with respect to $t$:

[TABLE]

Second, differentiate Eq.(2) with respect to $S$ instead:

[TABLE]

In what follows, we will also drop the 3 arguments of $\Sigma$ for notational ease, as it shouldn’t give rise to any confusion. Now differentiate Eq.(5) with respect to $S$:

[TABLE]

Substituting Eq.(4) and Eq.(6) into Eq.(3) implies:

[TABLE]

We now show that we can express Vanna $BS_{S\sigma}(S,\Sigma,t)$ and Volga $BS_{\sigma\sigma}(S,\Sigma,t)$ in terms of Vega $BS_{\sigma}(S,\Sigma,t)$. First, differentiate Eq.(2) with respect to $Z$:

[TABLE]

Next, differentiate Eq.(8) with respect to $S$:

[TABLE]

Now differentiate Eq.(8) with respect to $Z$ instead:

[TABLE]

Solving Eq.(10) for Volga $BS_{\sigma\sigma}(S,\Sigma,t)$ yields

[TABLE]

Substituting Eq.(11) into Eq.(9) and solving for Vanna $BS_{S\sigma}(S,\Sigma,t)$ implies:

[TABLE]

Substituting Eq.(11) and Eq.(12) into Eq.(7) yields

[TABLE]

Also from Eq.(8) we can express the option Vega via $\Sigma_{Z}$, from Eq.(5) - the option Delta, and substitute these expressions into Eq.(13) to obtain

[TABLE]

Therefore, if the option Vega is non-zero, the backward PDE for the implied volatility $\Sigma$ finally takes the form

[TABLE]

This is a quasilinear111A partial differential equation is said to be quasilinear if it is linear with respect to all the highest order derivatives of the unknown function, [Pinchover and Rubinstein, 2005]. PDE of the parabolic type.

When the Black-Scholes model Delta of a contingent claim is known in closed form, there is an alternative quasilinear PDE that governs the function $\Sigma(S,Z,t)$ on ${\cal D}$. Solving Eq.(8) for Vega $BS_{\sigma}(S,\Sigma,t)$, substituting the result into Eq.(5), and then solving for Delta $BS_{S}(S,\Sigma,t)$ implies:

[TABLE]

Hence Eq.(13) can also be written as:

[TABLE]

To illustrate, it is well known that for a European Call, Delta $BS_{S}(S,\Sigma,t)$ is known in closed form:

[TABLE]

where:

[TABLE]

and $N(x)$ is the normal CDF. Hence, for a European Call, the quasilinear PDE Eq.(17) can also be represented as another quasilinear PDE:

[TABLE]

whose domain is the Cartesian product of ${\cal R}_{E}$ and ${\cal A}_{C}$, which is $\mathcal{D}=(0,\infty)\times(0,T)\times((Q_{T}S-D_{T}K)^{+},Q_{T}S).$ Similar quasilinear PDE’s can be derived for European Puts or down-and-out Calls.

For an American put option with a positive early exercise premium, in the continuation region the relevant PDE is Eq.(15), since neither the value nor the Delta is known in closed form.

2 Forward equation

In practice it is often necessary to simultaneously compute the Black-Scholes implied volatility of many options written on the same underlying, but having different strikes and maturity. This is to construct so-called the implied volatility surface, given market quotes $Z(K,T)$ and the market values of $S,r,q$. However, it is well-known that solving this problem by using the backward PDE is time-consuming since for every pair $K,T$ a separate backward PDE has to be solved. In contrast, the forward equation can be used to efficiently do this in one sweep, see e.g., [Dupire, 1994, Gatheral, 2006].

Having this in mind, in this Section we derive a forward PDE for the implied volatility $\Sigma(S,T,K,r,q,Z)$. For easy of notation we will write it as $\Sigma(T,K,Z)$, thus dropping parameters that don’t change. The derivation can be done in a way similar to that in Section 1. Therefore, here we omit a detailed algebra, and provide just short explanations.

Again, we can start with the definition of the implied volatility $\Sigma(K,T,Z)$ in Eq.(2), which in new variables read

[TABLE]

where $\Sigma(K,T,Z)\in[0,\infty)$, and the domain of definition for $(K,T,Z)$ is ${\cal R}_{F,E}\equiv[0,\infty)\times(0,\infty)\times\mathcal{A}_{C}$.

To proceed, we need a forward PDE for the European option price $BS(K,T,\sigma)$ which is also known as Dupire’s equation, [Dupire, 1994]

[TABLE]

Since the PDE Eq.(22) holds for any level of $\sigma>0$, it holds, in particular, if we set $\sigma=\Sigma(K,T,Z)$. In what follows, we will drop the 3 arguments of $\Sigma$ for notational ease.

Differentiating Eq.(21) with respect to $T$, $K$ and twice with respect to $K$ yields

[TABLE]

Substituting Eq.(23) into Eq.(22) implies:

[TABLE]

Again, derivatives $BS_{K\sigma}(K,\Sigma,T)$ and $BS_{\sigma\sigma}(K,\Sigma,T)$ can be expressed in terms of $BS_{\sigma}(K,\Sigma,T)$. Differentiating Eq.(21), first with respect to $Z$, and then with respect either to $K$, or to $Z$, yields

[TABLE]

The last equation in this system can be solved for $BS_{\sigma\sigma}(K,\Sigma,T)$, and then the second one - for $BS_{K\sigma}(K,\Sigma,T)$. Substituting these solutions into Eq.(2) yields

[TABLE]

Using the first line in Eq.(25), this could be re-written as

[TABLE]

If the option Vega is non-zero, we finally obtain the forward PDE for the implied volatility $\Sigma(K,T,Z)$

[TABLE]

This PDE has to be solved subject to some initial and boundary conditions which are discussed in Section 3.

Since we consider the European options, the Black-Scholes model strike delta of a contingent claim is known in closed form. Then, similar to the previous Section, the PDE in Eq.(28) could be slightly simplified. Indeed, as follows from the second line of Eq.(23) and first line of Eq.(25)

[TABLE]

Hence Eq.(28) can also be re-written as:

[TABLE]

For instance, for a European Call, the strike delta $BS_{K}(K,\Sigma,T)$ can be found in closed form. Indeed, since the Black-Scholes European Call option price is a homogeneous function of order 1 in $(S,K)$ it is easy to find that

[TABLE]

where $d_{2}=d_{1}-\sigma\sqrt{T}$.

3 Initial and boundary conditions

The non-linear PDE Eq.(28) describes evolution of the function $\Sigma(K,T,Z)\in[0,\infty)$ at the domain $\mathcal{R}_{F,E}$, and has to be solved subject to some initial and boundary conditions. However, this immediately brings some problems.

Indeed, the initial condition for Eq.(28) has to be set at $T=0$. From the Black-Scholes formula we know that the option Vega $BS_{\sigma}(K,\Sigma,T)$ vanishes at $T=0$. Therefore, the expression in curly braces in Eq.(26) could have any value. Also, the equation $BS_{\sigma}(K,\Sigma,T)=0$ doesn’t have a solution for $\Sigma$. This is because at $T=0$ function $BS(K,\Sigma,0)$ gets its intrinsic value (e.g. for the Call option - $BS(K,\Sigma,0)=(S-K)^{+}$) which doesn’t depend on $\Sigma$. Accordingly, for the ITM options there is no $\Sigma$ which would make the intrinsic value zero, and for the OTM options any $\Sigma$ could be chosen.

To resolve this, we make a change of the dependent variable $\Sigma(K,T,Z)\mapsto{\mathfrak{S}}(K,T,Z)=\Sigma(K,T,Z)\sqrt{T}$. This pursues two goals. First, at $T=0$ a natural initial condition implies ${\mathfrak{S}}(K,0,Z)=0$. Second, if we know ${\mathfrak{S}}(K,T,Z)$ given some values of $(K,T,Z)$, the value of $\Sigma(K,T,Z)$ can be easily restored.

The forward PDE for the new dependent variable ${\mathfrak{S}}(K,T,Z)$ can be derived from Eq.(28), which yields

[TABLE]

Since Eq.(32) is a 2D parabolic PDE (quasilinear), the boundary conditions must be set $\forall T>0$.

Suppose that we consider an European Call defined at the domain $\Omega:[K,Z]\in[0,\infty]\times((Q_{T}S-D_{T}K)^{+},Q_{T}S)$ which is depicted in Fig. 1. A special attention should be drawn to two of these boundaries: $Z=0$ and $Z=Q_{T}S-D_{t}K$ - to decide whether the PDE needs any boundary condition at them.

As mentioned in [Itkin and Carr, 2011], the correct boundary conditions are determined by the speed of the diffusion term as we approach the boundary in a direction normal to the boundary. To illustrate, consider a PDE

[TABLE]

where $C=C(t,x)$ is some function of the time $t$ and the independent variable $x\in[0,\infty)$, $a(x),b(x),c(x)\in C^{2}$ are some known functions of $x$. Then, as shown by [Oleinik and Radkevich, 1973], no boundary condition is required at $x=0$ if $\lim_{x\rightarrow 0}[b(x)-a_{x}(x)]\geq 0$. In other words, the convection term at $x=0$ is flowing upwards and dominates as compared with the diffusion term. A well-known example of such consideration is the Feller condition as applied to the Heston model, see, e.g., [Lucic, 2008].

To make it clear, no boundary condition means that instead of the boundary condition at $x\rightarrow 0$ the PDE itself should be used at this boundary with coefficients $a(0),b(0),c(0)$.

In Section 5 we show how to solve Eq.(32) by using temporal discretization that gives rise to the following equation

[TABLE]

Here ${\mathfrak{S}}={\mathfrak{S}}(T+\Delta T,K,Z)$, $\tilde{{\mathfrak{S}}}=\tilde{\mathfrak{S}}(T,K,Z)$, and $\Delta T$ is the step of the temporal grid in $T$. The Eq.(34) is a linear PDE since all coefficients are already known (they are either given, e.g., $r,q,T$, or known from the previous time step). Thus, Eq.(34) is a 2D version of Eq.(33).

Accordingly, in the $K$ direction with allowance for Eq.(23),Eq.(31), the condition $\lim_{K\to 0}b(K)-a^{\prime}(K)\geq 0$ reads

[TABLE]

because $\lim_{K\to 0}\tilde{{\mathfrak{S}}}(K,T,Z)=0,\ T>0$. Therefore, as per [Oleinik and Radkevich, 1973], no boundary condition should be set at $K=0$ as the PDE itself serves as the correct boundary condition.

However, our domain $\Omega$ includes the only point with $K=0$, while the left boundary is a line $Z=Q_{T}S-D_{T}K$. It can be easily checked that along this line ${\mathfrak{S}}=0$ which is the required boundary conditions.

In the $Z$ direction, again using Eq.(25), we obtain

[TABLE]

From Eq.(25), ${\Sigma}_{Z}=1/BS_{\sigma}(K,{\tilde{{\Sigma}}},T)>0$. Therefore,

[TABLE]

Since at $Z=0$ the domain $\Omega$ is a straight line $KD_{T}\geq SQ_{T}$, it can be checked that the value of $d_{1}(K,{\tilde{{\Sigma}}},T)$ along this line is always negative, while ${\tilde{{\Sigma}}}\neq 0$ except the point $K=SQ_{T}/D_{T}$. Therefore, what remains to be checked is the behavior of the option Vega as $Z\to 0$. The implied volatility ${\tilde{{\Sigma}}}$ which solves the equation $Z=0$ can be found numerically, and is not zero at $KD_{T}>SQ_{T}$. Therefore, Vega is positive (despite very small). Thus, $\lim_{Z\to 0}\left[b(Z)-a^{\prime}(Z)\right]<0$. Therefore, as per [Oleinik and Radkevich, 1973], we need a boundary condition at the boundary $Z=0$. This condition can be found by solving the equation

[TABLE]

with respect to $\Sigma$. A typical example of the behavior of this solution is given in Fig. 2 which is computed using the following input data: $r=0.02,q=0.01$ and two cases with $(S=1.0,T=0.1)$ and $(S=100,T=1.0)$.

Boundary conditions at $K\to\infty$:

To set this boundary condition we need to know an asymptotic behavior of the implied volatility when $K\to\infty$. This was a subject of an intensive research for the last decade, and several useful results are available in the literature, see, e.g., [Lee, 2004, Benaim and Friz, 2009, Gulisashvili, 2010] and references therein. In particular, the following asymptotic formula was obtained in [Gulisashvili, 2010]

[TABLE]

as $K\to\infty$. Here, $C(K)$ is the Call option price corresponding to the given $K,T$, and $F=SQ_{t}/D_{T}$ is the forward price. One can observe, that in our case $C(K)=Z$. Therefore, as the boundary condition at $K\to\infty$ we can use

[TABLE]

Since $Z\leq SQ_{T}$, the logarithm $\log(SQ_{T}/Z)$ is always non-negative.

However, this asymptotics cannot be used along the whole vertical line of the domain $\Omega$ at $K\to\infty$ in Fig. 1. For instance, at $Z=SQ_{T}$ the reminder of the approximation in Eq.(37) is large. In this case one can solve numerically the equation

[TABLE]

Boundary conditions at $Z=SQ_{T}$:

The condition $Z=Q_{T}S$ implies that $N(d_{1}(K,T,{\mathfrak{S}}))=1$ and $N(d_{2}(K,T,{\mathfrak{S}}))=0$. This means that ${\mathfrak{S}}\to\infty$. To resolve this when constructing a numerical solution, one can move this boundary to $Q_{T}(S-\Delta S)$, where $\Delta S\ll S$. Then, using a Taylor series expansion of $BS(K,{\mathfrak{S}},T)$ around ${\mathfrak{S}}\to\infty$, we can find the solution of the equation

[TABLE]

Taking into account first 8 terms in this expansion, we arrive at the following equation

[TABLE]

where $k=SQ_{T}/(KD_{T}),\ \delta=\Delta SQ_{T}/(KD_{T})$, so $\delta\ll k$. This equation has three roots, and usually one is negative. Among the other two positive roots we need the smallest one, which approximately corresponds to the solution of Eq.(41) with $\delta=0$. Then this root aligns with the asymptotics of the solution at large $K$ considered in the previous paragraph. Setting $\delta=0$ in Eq.(41) gives rise to a cubic algebraic equation which can be solved analytically. Alternatively, Eq.(40) can be solved numerically.

4 Another type of PDE

The Eq.(32) can be further transformed into another type of the PDE. Indeed, the expression in square brackets in Eq.(32) can be represented as

[TABLE]

From Eq.(29), and Eq.(31) we have

[TABLE]

Substituting Eq.(43) into Eq.(42), and taking into account that ${\mathfrak{S}}={\mathfrak{S}}(K,T,Z)$, we obtain

[TABLE]

With this expressions, the PDE in Eq.(32) takes the form

[TABLE]

In contrast to Eq.(32), this is a nonlinear hyperbolic PDE. It can be used as an alternative to Eq.(32) which is a quasilinear parabolic PDE. Accordingly, the PDE in Eq.(45) requires just the boundary conditions at lines $K\in[SQ_{t}/D_{T},\infty)$ and $Z=(SQ_{T}-D_{T}K)^{+}$, see Fig. 1.

The Eq.(45) can be further simplified by using Eq.(43) that yields

[TABLE]

where $\Theta(T,K,\sigma)$ is the option Theta. This, of course, also follows from Eq.(23), Eq.(25), so Eq.(46) can be obtained directly from these equations. The Eq.(32) is a one-dimensional non-linear PDE with $K$ being a dummy variable. In other words, for any fixed $K$ this is a one-dimensional non-linear PDE with $T,Z$ be the independent variables.

5 Numerical method

In this Section we describe a numerical method used to solve either the PDE in Eq.(32) or in Eq.(45). The method is constructed similar to how this is done in [Itkin, 2018].

First, we re-write Eq.(32) and Eq.(45) in the operator form by introducing an operator ${\mathcal{L}}(T,K,Z):[0,\infty)\mapsto[0,\infty)$. With this notation Eq.(32) can be represented as

[TABLE]

where

[TABLE]

Similarly, Eq.(45) can be re-written in the same form with

[TABLE]

Suppose we now discretize the time to maturity $T\in[0,\infty)$ by creating a uniform grid in time with step $\Delta T$, so our discrete time is defined at the points $[0,\Delta T,2\Delta T,\ldots)$. By using Taylor series expansion it can be shown that the discrete solution of the Eq.(47) up to $O((\Delta T))$ could be represented in the form

[TABLE]

where $e^{\Delta T\mathcal{L}(T,K,Z)}$ is the operator exponent, [Hochschild, 1981].

Again, using the Taylor series expansion, one can verify that the following scheme

[TABLE]

also approximates Eq.(47) up to in $O((\Delta T))$. Combining them together we obtain the scheme

[TABLE]

which approximates Eq.(47) up to $O((\Delta T)^{2})$.

Now we setup an iterative algorithm to solve this discrete equation. This algorithm is a version of the fix-point Picard iteration algorithm, and relies on the Banach fixed-point theorem, [Granas and Dugundji, 2003]:

At the first iteration we start by setting $\mathcal{L}(T+\Delta T,K,Z)=\mathcal{L}(T,K,Z)$ as the initial guess. Thus, Eq.(52) transforms to Eq.(50), which could be represented in the form

[TABLE]

where in ${\mathfrak{S}}^{(k)}$ the superscript (k) marks the value of ${\mathfrak{S}}$ found at the $k$-th iteration of the numerical procedure. As the exponent $\mathcal{L}(T,K,Z)$ in the right hands size of Eq.(53) is computed at the known time $T$, it is a linear operator with coefficients being the known functions of $(K,T,Z)$. In particular, when $T=0$ these functions are determined by the initial condition. If $T>0$ they are determined by the solution at the previous time step.

The Eq.(53) could be solved by either computing the discrete matrix exponential defined on some grid in $(K,Z)$ (which could be done with complexity $O(N^{3}M+NM^{3})$, $N$ is the number of the grid points in $K$ direction, and $M$ is the number of the grid points in $Z$ direction), or by using any sort of Padé approximation of the exponential operator. For instance, the Padé approximation $(1,1)$ leads to the well-known Crank-Nicholson scheme which could be solved with the linear complexity in $N$ and $M$, and provides $O((\Delta T)^{2})$ approximation in time, see [Itkin, 2017b]. 2. 2.

To proceed we represent Eq.(52) in the form

[TABLE]

Therefore, the next approximation ${\mathfrak{S}}^{(k)}(T+\Delta T,K,Z)$ to ${\mathfrak{S}}(T+\Delta T,K,Z)$ can be found again by either computing the matrix exponential, or by using some Padé approximation.

Having this machinery in hands we can proceed in the same manner until the entire procedure converges, i.e. when the condition

[TABLE]

is reached after $k$ iterations with $\varepsilon$ being the method tolerance, and $\|\cdot\|$ being some norm, e.g., $L^{2}$.

Note, that Eq.(54) with the accuracy $O(({\Delta T})^{2})$ can be re-written in the form

[TABLE]

where tilde in the new notation $\tilde{\mathcal{L}}^{(k-1)}$ means that we take a half of sum of the operators $L$ at times $T$ and $T+\Delta T$, and (k-1) means that at time level $T+\Delta T$ we take the value of the operator at the $(k-1)$-th iteration.

5.1 Coordinate transformation

Since our PDEs are defined in the domain $\Omega$ which has a trapezoidal shape, we cannot directly apply a finite-difference (FD) method for solving Eq.(54), while finite elements or Radial Basis Functions methods can be used here with no problem. However, a simple trick makes application of the FDM possible. The trick utilizes the fact that the right vertical boundary of $\Omega$ lies at $K\to\infty$. Thus, based on the asymptotic formula Eq.(38), in this limit ${\mathfrak{S}}\to\infty$ as well, regardless of the value of $Z\in[0,SQ_{T}/D_{T}]$. Therefore, without any loss of accuracy, the domain $\Omega$ could be replaced with the domain $\tilde{\Omega}$ depicted in Fig. 3.

By construction, the domain $\tilde{\Omega}$ is a parallelogram, and $K\in[0,K_{\max}]$, where $K_{\max}$ is some fixed strike. We introduce it to truncate an infinite strip in $K$ to some fixed computational area.

Now, the parallelogram could be transformed to rectangle by using, e.g., the following affine transformation (for introduction to this topic see, e.g., [Katsumi; and Sasaki, 1994])

[TABLE]

This transformation converts the domain $\tilde{\Omega}$ into a unit square with the following map of the boundaries

[TABLE]

Accordingly, we need to re-write the operator ${\mathcal{L}}(T,K,Z)$ in new variables $(u,v)$. For instance, for the parabolic operator Eq.(48) this yields

[TABLE]

Since $\tilde{\Omega}$ - the domain of definition of the operator ${\mathcal{L}}(T,v,u)$ in Eq.(58) - is a unit square, now the PDE can be solved by using a FD method.

5.2 Splitting in spatial dimensions

Since computation of matrix exponent is time-consuming, especially in case of a two-dimensional operator $L$, we proceed with using a splitting technique, see [Lanser and Verwer, 1999, Duffy, 2006, Itkin, 2017b] and references therein. For instance, in [Itkin, 2017a] an ADI (alternative direction implicit) scheme proposed in [In’t Hout and Welfert, 2007] for the solution of a backward PDE was extended twofold. First, a similar ADI scheme has been proposed for the forward equation which equalizes the results obtained by solving backward and forward PDEs. Second, it replaces the first explicit step of the ADI, applied to the mixed derivative term of the PDE, with the fully implicit step. This is especially important for the forward PDE as it increases the stability of the solution and guarantees its positivity.

However, for the PDE in Eq.(55) using this ADI approach is not a trivial problem, mainly because of various problems with a stable approximation of the mixed derivative term at the second sweep of the ADI. Therefore, instead of the ADI method, here we use the Strang splitting, [Lanser and Verwer, 1999, Duffy, 2006, Itkin, 2017b]. The main idea of this approach is as follows.

With allowance for Eq.(58), let us represent the Eq.(55) in the new coordinates in the form

[TABLE]

Thus, the operator ${\mathcal{L}}_{u}$ includes all differentials with respect to $u$ and also a half of the killing term in Eq.(58). This is a linear convection-diffusion operator with state and time dependent coefficients. The operator ${\mathcal{L}}_{v}$ includes all differentials with respect to $v$ and a half of the killing term in Eq.(58), and also is a linear convection-diffusion operator with state and time dependent coefficients. Finally, the operator ${\mathcal{L}}_{uv}$ includes only a mixed derivatives term in Eq.(58), and also is a linear diffusion operator with state and time dependent coefficients.

Equation Eq.(60) can be solved numerically by using the Strang splitting scheme

[TABLE]

This scheme can also be represented as a sequence of fractional steps, which read

[TABLE]

We underline that the Strang splitting provides the second order approximation in ${\Delta T}$ only if the PDE at every fractional step in Eq.(62) is solved also with the second order approximation in $\Delta T$. To achieve that at every fractional step we use the Padé approximation $(1,1)$ to obtain a Crank-Nicholson scheme in time. For instance, for the first line in Eq.(62) it gives the following scheme

[TABLE]

5.3 Spatial discretization

To apply this algorithm we construct an appropriate discrete nonuniform grid ${\mathbf{G}}(x)$ in the $(u,v)$ space, see [Itkin, 2017b]. In each spatial direction $i\in[u,v]$ a nonuniform grid contains $N_{i}+1$ nodes ($p_{0},p_{1},...,p_{N_{i}}$) with the steps $h_{1}=p_{1}-p_{0},...,h_{N_{i}}=p_{N_{i}}-p_{N_{i}-1}$. Accordingly, in the operator Eq.(58) we replace continuous derivatives with their finite-difference approximations, and denote these discrete operators defined on ${\mathbf{G}}(x)$ as $\triangledown$. Thus, we replace ${\partial}_{x}$ with $\triangledown_{x}$, etc.

5.3.1 Approximation of derivatives

Once a non-uniform grid is constructed, the first and second derivatives can be approximated on this grid. The corresponding expressions are obtained by using the Taylor series expansions. Here for the reference we provide just the final results, see, e.g., [Itkin, 2017b, In’t Hout and Foulon, 2010]:

Let $f:\mathbb{R}\to\mathbb{R}$ be any given function, let ${x_{i}},i\in\mathbb{Z}$ be any given increasing sequence of mesh points, and $h_{i}=x_{i}-x_{i-1},\forall i$. To approximate the first derivatives of $f(x)$ employ the following formulas.

[TABLE]

Here the super index 1X means the first derivatives with the $X$ approximation, $X\in[B,F,C]$ could be backward (B), forward (F) and central (C), and the coefficients $\alpha,\beta,\gamma$ read

[TABLE]

As follows from Eq.(5.3.1), these expressions provide approximation of the first derivatives with the second order $O(h^{2})$.

To approximate the second derivatives of $f(x)$ one can employ the following formulas.

[TABLE]

These expressions provide approximation of the second derivatives also with the second order $O(h^{2})$. For instance, for $\triangledown^{2C}_{x}$ we have in the explicit form

[TABLE]

where

[TABLE]

Let $f:\mathbb{R}^{2}\to\mathbb{R}^{2}$ be any given function of two variables $(x,y)$. To approximate the mixed derivative $\partial_{x,y}f(x_{i},y_{j})$ at any point $(x_{i},y_{j})$ one can consequently apply either one-sided approximations in each direction, or apply the central difference approximation first in $x$ and then in $y$ (or vice versa). The latter results in the FD approximation of the mixed derivative using a 9-points stencil while still providing the second order of approximation. However, as shown in [Itkin, 2017a], this discretization could be unstable. Therefore, instead a fully implicit scheme was proposed which, as applied to our problem, is discussed in the next Section.

For the future we need some additional notation. We denote a matrix of $\triangledown^{1F}_{x}$ as $A^{2F}_{1,x}$; the matrix of $\triangledown^{1B}_{x}$ - as $A^{2B}_{1,x}$; the matrix of $\triangledown^{1C}_{x}$ - as $A^{2C}_{1,x}$ 222Hence, e.g. for $A^{2F}_{1,x}$, the superscript means: backward second order approximation, and subscript means: first order derivative.. For the first order approximations we use the following definitions. Define a one-sided forward discretization of $\triangledown$, which we denote as $A^{1F}_{1,x}:\ A^{1F}_{1,x}C(x)=[C(x+h,t)-C(x,t)]/h$. Also define a one-sided backward discretization of $\triangledown$, denoted as $A^{1B}_{1,x}:\ A^{1B}_{1,x}C(x)=[C(x,t)-C(x-h,t)]/h$. These discretizations provide first order approximation in $h$, e.g., $\partial_{x}C(x)=A^{1F}_{1,x}C(x)+O(h)$. Also $I_{x}$ denotes a unit matrix.

5.3.2 Approximation of the operator ${\mathcal{L}}(T,v,u)$

Here we use the same approach as in [Itkin, 2018]. Denote $A^{2C}_{2,x}$ a matrix of the discrete operator $\triangledown^{2C}_{x}$ on the grid ${\mathbf{G}}(x)$. Observe, that $A^{2C}_{2,x}$ is the Metzler matrix, see [Berman and Plemmons, 1994, Itkin, 2017b]. Indeed, it has all negative elements on the main diagonal, and all nonnegative elements outside of the main diagonal. Also $A^{2C}_{2,x}$ is a tridiagonal matrix.

By the properties of the Metzler matrix $M$ its exponent is a nonnegative matrix. Therefore, operation ${e^{M}}{\mathfrak{S}}$ preserves positivity of vector ${\mathfrak{S}}$. Also all eigenvalues of the Metzler matrix have a negative real part. Therefore, the spectral norm of the matrix follows

[TABLE]

Thus, if in Eq.(54) the operator $\mathcal{M}=\tilde{\mathcal{L}}^{(k-1)}$ on the grid ${\mathbf{G}}(x)$ is represented by the Metzler matrix $A_{\mathcal{M}}$, the map $e^{A_{\mathcal{M}}}$ is contractual, and, hence, the entire scheme Eq.(54) is stable.

By the definition of the operator $\tilde{\mathcal{L}}^{(k-1)}$ in Eq.(60), it consists of several terms.

Second derivatives.

The first term $\frac{1}{2T}{\mathfrak{S}}^{2}K_{1}^{2}{\mathfrak{S}}^{2}_{u}{\partial}^{2}_{v}$ can be approximated on the grid as

[TABLE]

As $T>0$, obviously $A_{2,v}$ is the Metzler matrix. Same is true for the matrix $A_{2,u}$:

[TABLE]

The mixed derivative.

The mixed derivative term in Eq.(60) reads: $-\frac{1}{T}{\mathfrak{S}}^{2}K_{1}^{2}{\mathfrak{S}}_{u}{\mathfrak{S}}_{v}{\partial}_{u}{\partial}_{v}$. Observe, that from Eq.(56) we have

[TABLE]

from Eq.(25)

[TABLE]

and from Eq.(29), Eq.(31)

[TABLE]

With the allowance for these expressions, it follows that

[TABLE]

where $R(x)\equiv\frac{\displaystyle 1-N(x)}{\displaystyle N^{\prime}(x)}$ is the Mills ratio for the standard normal distribution. It can be efficiently computed by the particularly simple continued fraction representation at $x>1$

[TABLE]

or by using Taylor series expansion at $0\leq x\leq 1$, see [Gasull and Utzet, 2013].

Thus, ${\mathfrak{S}}_{u}{\mathfrak{S}}_{v}\geq 0$. Therefore, to have $A_{2,uv}$ to be the Metzler matrix, we need the matrix of ${\mathfrak{S}}_{u}{\mathfrak{S}}_{v}$ to be the negative of the Metzler matrix 333Rigorously speaking, the matrix $A_{2,uv}$ should be the negative of the EM-matrix, where the latter stays for the Eventually M-Matrix. For the introduction into the EM matrices and their properties see [Itkin, 2017b]., again see [Itkin, 2017b]. However, the standard methods fail to provide discretization of the second order that guarantees this feature. Therefore, we need another approach, similar to what was proposed in [Itkin, 2017a]. We discuss this in Appendix A.

As compared with [Itkin, 2017a], the method in this paper has two innovations. First, in [Itkin, 2017a] a first order approximation in time step was constructed, while here we provide the second order approximation. Second, in [Itkin, 2017a] a mixed derivative term has the form $\rho_{u,v}f(v)g(u)\partial_{u}\partial_{v}$, where $f(u)$ is some function of $u$ only, and $g(v)$ is some function of $v$ only. Here we deal with a more general term $\rho_{u,v}f(u,v)\partial_{u}\partial_{v}$, and show that the approach of [Itkin, 2017a] can be applied in this case as well.

Convection terms.

The next step is to construct a suitable approximation for the convection terms. The term $qu\partial_{u}$ admits the forward approximation $\triangledown^{1F}_{u}$ as $q\geq 0$ and $u\in[0,1]$. Approximating the term $-[(r-q)v+B]{\partial}_{v}$ is a more delicate issue and depends on the sign of the function $F(u,v)=(r-q)v+B$. Suppose $r>q$. A typical behavior of function $F(u,v)$ at $T=0.01$ is presented in Fig. 4 at the values of model parameters given in Table 1.

As can be seen, $F(u,v)\geq 0$ for all values of $(u,v)\in[0,1]\times[0,1]$. Therefore, when approximating the term $-F(u,v){\partial}_{v}$. a backward approximation should be used to obtain the Metzler matrix, i.e. $-F(u,v)A^{1B}_{2,v}$. It can be checked that in case $r<q$ the forward approximation should be used as then $F(u,v)\leq 0$.

Finally, the resulting matrices $A_{\tilde{\mathcal{L}}_{u}},A_{\tilde{\mathcal{L}}_{v}},$ are, by construction, also the Metzler matrices, as each of them is a sum of the Metzler matrices, [Berman and Plemmons, 1994].

It is easy to prove that the matrix $\exp(A_{\tilde{\mathcal{L}}^{(k-1)}})$ is the second order approximation of the operator $e^{\mathcal{M}}$ in the spatial steps $h_{u},h_{v}$ on the grid, i.e. $O(h_{u}^{2}+h_{v}^{2}+h_{u}h_{v})$. That follows from the fact that the derivatives are approximated by using Eq.(5.3.1), where every line is the second order approximation to the corresponding partial derivative.

5.4 Convergence of Picard iterations

The remaining point to prove is the convergence of the fixed-point Picard iterations. Again, in doing so we follow [Itkin, 2018]. According to Eq.(54), change of variables in Eq.(56), and Eq.(55) we run the iterations

[TABLE]

By construction, at every iteration $k$ the operator $\tilde{\mathcal{L}}^{(k)}$ is linear, continuous and bounded at the computational domain. Also, the operator ${\mathcal{E}}=\exp[{\Delta T}\tilde{\mathcal{L}}^{(k)}]$ is Lipschitz

[TABLE]

since we constructed discretization of ${\mathcal{E}}$ as $e^{M}$ with $M$ being the Metzler matrix, and it was shown earlier that in the spectral norm $\|e^{M}\|\leq 1$.

Now, to prove convergence of the Picard iterations denote by $\hat{{\mathfrak{S}}}(T+\Delta T,u,v)$ the exact solution of Eq.(47) at time $T+\Delta T$. Then, by the mean-value theorem for operators, we have

[TABLE]

where $\xi^{(k)}$ is a convex combination of ${\mathfrak{S}}^{(k-1)}(T+\Delta T,u,v)$ and ${\mathfrak{S}}^{(k)}(T+\Delta T,u,v)$, and $\mathbb{D}$ denotes the Frèchet derivative of the operator ${\mathcal{E}}$ at the space of all bounded non-linear operators, see, e.g., [Hutson et al., 2005, Li et al., 2010]. Thus, the iterations converge if $\|\mathbb{D}({\mathcal{E}})\|<1$. We prove that actually this is the case for the operators in Eq.(60) in Appendix B. In particular, we show there that the iterative scheme is "almost" unconditionally stable, where "almost" means for the values of the model parameters reasonable from the practical point of view.

All the results obtained in this Section can be further summarized in the following proposition

Proposition 5.1**.**

The scheme Eq.(54) with approximation Eq.(67) is a) (almost) unconditionally stable, b) preserves non-negativity of the solution, c) provides second order approximation of Eq.(54) in space on the grid $\mathbb{G}$, and d) converges.

6 Numerical results and comparison

In our numerical test we use the input parameters shown in Table 1.

We build a non-uniform grid ${\mathbf{G}}(x)$ in the $(u,v)$ space, which contains $N_{u}$ nodes in the $u$ direction and $N_{v}$ nodes in the $v$ direction. Since the gradients of the implied volatility ${\mathfrak{S}}(T,u,v)$ are high at the boundaries of the computational domain, i.e. close to $u=0,1$ and $v=0,1$ it makes sense to make the grid denser close to these areas. Since from practical purposes the far boundary is of a less interest, we compress the grid close to the lines $u=0,\ v=0$. This is obtained by using the compression ratio $R_{c}=50$, see the description in [Itkin, 2017b]. Thus constructed grid is represented in Fig. 6. With the initial parameters in Table 1, one step in the $v$ direction approximately corresponds to a change in strike by 3$, and one step in the$u$ direction approximately corresponds to a 50 cents change in both the options price and the strike.

We solve the forward equation starting at $T=0$ and going forward in time to $T_{\max}$ with the time step $\Delta T=0.01$. The boundary conditions are computed as this is described in Section 3 using numerical solutions of Eq.(36), Eq.(39) and Eq.(40). Thus found boundary values corresponding to $T=0.01$ are presented in Fig. 6 (along the lines $u=0$, $u=1$, $v=0$ and $v=1$).

The standard way to get the implied volatility by solving Eq.(2) numerically is used for comparison, and further is called "traditional". We run this traditional approach at every point of the grid ${\mathbf{G}}(x)$ by using Matlab function blsimpv, or a root solver at the boundaries where blsimpv provides unsatisfactory results. All tests were run on a PC with Intel i7-4790 CPU (8 Cores, 3.6GHz) as a sequential loop, i.e. no parallel tools have been used. The total elapsed time was 113 secs for 20000 points on the grid, or, on average, 6 msc per a single point.. The results are presented in Fig. 3.

The grid ${\mathbf{G}}(x)$ doesn’t change in time. On the other hand, functions $d_{1}({\mathfrak{S}}),d_{2}({\mathfrak{S}})$ in the Black-Scholes formula (see Eq.(19)) also almost don’t depend on $T$ (since $T$ is now embedded into ${\mathfrak{S}}$). Therefore, the only dependence on $T$ in the numerical implementation of Eq.(2) is the dependence of $Q_{T},D_{T}$ which is weak. Thus, Fig. 6 changes with time very slow. Also, for the purpose of a better illustration of the behavior of ${\mathfrak{S}}$ at intermediate values of $(u,v)$, we zoom-in Fig. 6 by excluding the boundaries and setting $T=0.5$. Thus obtained plot is represented in Fig. 7.

6.1 First step in time

The first step in time is the most challenging problem with our approach. Indeed, according to Section 3, the initial condition for the PDE in Eq.(55) is ${\mathfrak{S}}(0,u,v)=0$ for all $u$ and $v$. This, in turn, means that all gradients ${\mathfrak{S}}_{u},{\mathfrak{S}}_{v}$ tend to infinity at $T=0$. Also,

[TABLE]

where ${\mathcal{L}}^{\mathrm{NL}}$ denotes the non-linear part of the operator ${\mathcal{L}}$. Therefore, the iterative scheme in Eq.(55) becomes incorrect. At the moment we don’t have a good resolution of this problem. Therefore, instead, for $T={\Delta T}$ we proceed by computing ${\mathfrak{S}}(\Delta T,u,v)$ in the traditional way, and denote thus obtained implied volatility as ${\mathfrak{S}}_{tr}(\Delta T,u,v)$.

Again, the main problem here is that at $T={\Delta T}$ we don’t have a good initial guess for ${\mathfrak{S}}$ as it cannot be taken from the previous time step. But, in principle, if we would have a good initial guess for ${\mathfrak{S}}$, we could run our iterative algorithm even at $T={\Delta T}$. To validate this, we further use ${\mathfrak{S}}_{tr}(\Delta T,u,v)$ as the initial guess for the iterative scheme. But, based on Eq.(68), for this time step we also need to modify the scheme in Eq.(55) to be

[TABLE]

Then we solve this equation iteratively, but it is easy to check that this scheme provides only the first order approximation in $\Delta T$.

The results are as follows. We do 20 iterations of the method to obtain the relative error

[TABLE]

to be 1.5 bps. The elapsed time for doing so is 16 secs, or 0.8 secs per iteration. The difference between ${\mathfrak{S}}_{tr}(\Delta T,u,v)$ and our numerical solution ${\mathfrak{S}}(\Delta T,u,v)$ is presented in Fig. 8. The graph is truncated up to $u=v=0.8$ because very large values of strikes don’t have practical interest.

6.2 Next steps in time

Having good (not all zeros) values of ${\mathfrak{S}}(\Delta T,u,v)$ we can proceed based on the general scheme described in Section 5. We increase time to $T=2\Delta T$, iterate to get the numerical solution, then compute ${\mathfrak{S}}_{tr}(2\Delta T,u,v)$, compare these solutions, etc. In Fig. 9 these results are presented for $T=2\Delta T,5\Delta T,10\Delta T,20\Delta T$.

It turns out that at those time steps the iterations converge much faster as compared with the test for the first step reported in Section 6.1. This convergence is summarized in Table 2. This partially could be explained by the order of approximation in time, since at the first step the scheme was of the order one, while at the next steps it is of the order 2.

As this can be seen from the presented results, the convergence of the method slightly improves with time $T$, the elapsed time remains the same, while the relative error between the traditional method and the PDE solver slightly increases with $T$. The later can obviously be explained by accumulation of the approximation errors at every time step with the increase of $T$. However, for instance, for $T=0.2$ this error is still acceptable (to be 20-50 bps). On the other hand, the biggest deviation from the traditional results corresponds to high $u$ and $v$. In particular, in our numerical example the point $u=v=0.8$ maps to the strike K=260\$$ which, in turn, corresponds to the moneyness S/M=0.385$and the option price 20$. Hence, this is a deep OTM option, and the corresponding implied volatility is very high, 2.61. Therefore, this case is almost impractical.

7 Discussion

In this paper we derive a backward and forward quasilinear parabolic PDEs that govern the implied volatility of a contingent claim whenever the latter is well-defined. This would include at least any contingent claim written on a positive stock price whose payoff at a possibly random time is convex. Alternative to that we also derived a forward nonlinear hyperbolic PDE of the first order which also governs evolution of the implied volatility in $(K,T,Z)$ space. We also discuss suitable initial and boundary conditions for those PDEs. Finally, we demonstrate how to solve them numerically by using an iterative method. All these results are new.

We compare performance of two methods. The first one uses a traditional scheme of getting the implied volatilities by using a root solver and solving Eq.(2). The second method relies on a numerical solution of the derived PDEs. We show that based on our numerical experiments, if the implied volatility is needed at every point on the rather dense grid, the PDE solver significantly outperforms the root solver. As follows from Table 2, the PDE solver is almost 40 times faster.

However, a few comments shoud be made with this regard. First, our method utilizes the traditional method in two ways. We use it to compute the boundary conditions, and we also use it at the first time step $T=\Delta T$ as this was discussed in Section 6.1 since we don’t have a good start value of ${\mathfrak{S}}$ for iterations to kick off. Also this drops the order of approximation in time from $O((\Delta T)^{2})$ to $O(\Delta T)$ at this step. Nevertheless, all those issues are relatively minor, and for the longer time horizons our method still outperforms the traditional one.

The second obstacle is that, most likely, in practice we rarely need implied volatilities on so dense grid in both space and time. Usually, we need them only at those points in $(K,T)$ which correspond to the tradable market strikes and maturities, and they are not so dense. However, the FD method still requires a grid. Therefore, the number of points for the traditional method could be much smaller that that for the FD method. This could compensate a slower performance of the traditional method and make both methods to be comparable in the elapsed time.

Third, as was already mentioned, functions $d_{1}({\mathfrak{S}}),d_{2}({\mathfrak{S}})$ in the Black-Scholes formula (see Eq.(19)) weakly depend on $T$ (since $T$ is now embedded into ${\mathfrak{S}}$), namely only via $Q_{T}$ and $D_{T}$. Therefore, the solution of Eq.(2) also depends on $T$ only via $Q_{T},D_{T}$. Therefore, since we already computed ${\mathfrak{S}}(\Delta T,K,Z)$ by using the traditional way, then for $T=2\Delta T$ Eq.(2) could be solved by expanding the left hands side into series on $\Delta T\ll 1$. In more detail, we write

[TABLE]

and define ${\mathfrak{S}}(T+\Delta T,K,Z)={\mathfrak{S}}(T,K,Z)+\Delta{\mathfrak{S}}(T,K,Z)$. Now expanding the second line of Eq.(70) into series on $\Delta T\ll 1,\Delta{\mathfrak{S}}\ll 1$, and using the first line in Eq.(70), in the first order we obtain

[TABLE]

Since ${\mathfrak{S}}(T,K,Z)$ is already known, the expression in Eq.(71) can be easily computed.

Once all ${\mathfrak{S}}(T+\Delta T,K,Z)$ are computed in such a way, we can proceed to $T+2\Delta T$. The problem with this approach, however, is that its accuracy decreases with every time step. As an alternative, periodically at some time steps the accurate value of ${\mathfrak{S}}(T,K,Z)$ can be recomputed by using the traditional approach. However, same could be done for the PDE approach as well.

It should also be mentioned that in this paper we present results for European plain vanilla options using the forward equation, while for e.g., for American options one has to solve the backward PDE, which computationally is more expensive.

Acknowledgments

We thank Archil Gulisashvili and Daniel Duffy for useful discussions. Any errors are our own.

Appendix A FD approximation of the mixed derivative term in Eq.(60)

In this Appendix we show how to solve the third line in Eq.(62) with the accuracy $O(({\Delta T})^{2})$ in time and $O(h_{u}^{2}+h_{v}^{2}+h_{u}h_{v})$ in space folloiwng the method first proposed in [Itkin, 2017a]. To recall, the equation under consideration reads

[TABLE]

To achieve this accuracy in time we use the Padé approximation $(1,1)$ to obtain a Crank-Nicholson scheme (wrtten with some loose of notation to make it better readable)

[TABLE]

This numerical scheme is, of course, implicit due to the fact that, despite we know $\bar{\mathfrak{S}}(T)$, to find ${\mathfrak{S}}(T+{\Delta T})$ we need to solve Eq.(A.2). It is known that, after some discretization on the FD grid is applied to Eq.(A.2), the matrix in the left hands side of Eq.(A.2) must be an M-matrix, to ensure that this scheme is stable and preserves the positivity of the solution. However, this is impossible by using the standard FD discretizations of the second order, in more detail see [Itkin, 2017b].

Therefore, we rewrite Eq.(A.2) by dropping the dependence on $T$ and using a trick proposed in [Itkin, 2017a]444The trick is motivated by the desire to build an ADI scheme which consists of two one-dimensional steps, because for the 1D equations we know how to make the rhs matrix to be an EM-matrix, see [Itkin, 2016, Itkin, 2017b] and references therein.

[TABLE]

Here $\rho_{u,v}=-1/4,$555We introduce this notation to make it easy to compare the description of the method in this paper with that in [Itkin, 2017a]. $P,Q,$ are some positive numbers which have to be chosen based on some conditions, e.g., to provide diagonal dominance of the matrices in the parentheses in the lhs of Eq.(A.3), see below. Also the intuition behind this representation is discussed in detail in [Itkin, 2017a].

Note, that in [Itkin, 2017a] we considered a mixed derivative term of the form $\rho_{u,v}f(v)g(u)\partial_{u}\partial_{v}$, where $f(u)$ is some function of $u$ only, and $g(v)$ is some function of $v$ only. Here we are dealing with a more general term $\rho_{u,v}f(u,v)\partial_{u}\partial_{v}$, so the presented approach is a generalization of that in [Itkin, 2017a]. Also our method in this paper is of the second order in time while the method in [Itkin, 2017a] is of the first order.

The Eq.(A.3) can be solved using fixed-point Picard iterations. One can start with setting ${\mathfrak{S}}^{0}=\bar{\mathfrak{S}}$ in the right hands side of Eq.(A.3), then solve sequentially two systems of equations

[TABLE]

Here ${\mathfrak{S}}^{j}$ is the value of ${\mathfrak{S}}$ at the $j$-th iteration666The reader shouldn’t miss these iterations (marked by $j$) used to solve Eq.(A.1) with those (marked by $k$) for the solution of the nonlinear equation which are discussed in Section 5..

To solve Eq.(A) we propose two FD schemes. The first one (Scheme A) is introduced by the following Propositions777For the sake of clearness we formulate this Proposition for the uniform grid, but it is transparent how to extend it for the non-uniform grid.:

Proposition A.1**.**

Let us approximate the left hands side of Eq.(A) using the following finite-difference scheme:

[TABLE]

Then this scheme is unconditionally stable in time step ${\Delta T}$, approximates Eq.(A) with $O(\sqrt{{\Delta T}}\max(h_{u},h_{v}))$ and preserves positivity of the matrix ${\mathfrak{S}}(u,v,T)$ if $Q=\beta\sqrt{\Delta T}/h_{v},\ P=\beta\sqrt{\Delta T}/h_{u}$, where $h_{v},h_{u}$ are the grid space steps correspondingly in $v$ and $u$ directions, and the coefficient $\beta$ must be chosen to obey the condition:

[TABLE]

Proof.

The proof can be found in [Itkin, 2017a, Itkin, 2017b]. ∎

The computational scheme in Eq.(A.5) should be understood in the following way. At the first line of Eq.(A.5) we begin with computing the product ${\mathfrak{S}}_{1}=\alpha^{+}\bar{\mathfrak{S}}$. This can be done in three steps. First, the product ${\mathfrak{S}}_{1}=Q\sqrt{\Delta T}\rho_{u,v}U(u,v)A^{1F}_{1,u}\bar{\mathfrak{S}}$ is computed in a loop on $v_{i},i=1,...,N_{v}$. In other words, if $\bar{\mathfrak{S}}$ is a $N_{u}\times N_{v}$ matrix where the rows represent the $u$ coordinate and the columns represent the $v$ coordinate, each $j$-th column of ${\mathfrak{S}}_{1}$ is a product of matrix $Q\sqrt{\Delta T}\rho_{u,v}U(u,v)A^{1F}_{1,u}$ and the $j$-th column of $\bar{\mathfrak{S}}$. The second step is to compute the product ${\mathfrak{S}}_{2}=P\sqrt{\Delta T}V(u,v)A^{1F}_{1,v}\bar{\mathfrak{S}}$ which can be done in a loop on $u_{i},\ i=1,...,N_{u}$. Finally, the right hand side of the first line in Eq.(A.5) is $(PQ+1)\bar{\mathfrak{S}}-{\mathfrak{S}}_{1}+{\mathfrak{S}}_{2}-{\mathfrak{S}}^{j}$. Then in a loop on $u_{i},\ i=1,...,N_{u}$, $N_{u}$ systems of linear equations have to be solved, each giving a row vector of ${\mathfrak{S}}^{*}$. The advantage of the representation Eq.(A.5) is that the product $\alpha^{+}\bar{\mathfrak{S}}$ can be precomputed.

If $\sqrt{\Delta T}\approx\max(h_{u},h_{v})$, then the whole scheme becomes of the second order in space. However, this would be a serious restriction inherent to the explicit schemes. Therefore, we don’t rely on it. Note, that in practice the time step is usually chosen such that $\sqrt{\Delta T}\ll 1$, and hence the whole scheme is expected to be closer to the second, rather than to the first order in $h$. Note, that this approach is similar to [Toivanen, 2010, Chiarella et al., 2008] for negative correlations, where a seven point stencil breaks a rigorous second order of approximation in space.

As shown in [Itkin, 2017a], the rate of convergence of the Picard iterations is linear.

The above results can be further improved by making the whole scheme to be of the second order of approximation in $h_{u}$ and $h_{v}$. We call this FD scheme as Scheme B.

Proposition A.2**.**

Let us approximate the lhs of Eq.(A) using the following finite-difference scheme:

[TABLE]

Then this scheme is unconditionally stable in time step ${\Delta T}$, approximates Eq.(A) with $O(\max(h^{2}_{u},h^{2}_{v}))$ and preserves positivity of the vector ${\mathfrak{S}}(u,x,T)$ if $Q=\beta\sqrt{\Delta T}/h_{v},\ P=\beta\sqrt{\Delta T}/h_{u}$, where $h_{v},h_{u}$ are the grid space steps correspondingly in $v$ and $u$ directions, and the coefficient $\beta$ must be chosen to obey the condition:

[TABLE]

The scheme Eq.(A.6) has a linear complexity in each direction.

Proof.

See [Itkin, 2017a]. ∎

The coefficient $\beta$ should be chosen experimentally. In our experiments described in the following sections we used

[TABLE]

Appendix B Norm of the Frèchet derivative for operators in Eq.(60)

To compute the norm of the Frèchet derivative recall that by definition

[TABLE]

If $\nu=(\nu_{1},...,\nu_{m})\in[L^{\infty}(-\infty,0)]^{m}$, then

[TABLE]

In other words, in this case the Frèchet derivative coincides with the Gâteaux derivative.

To find $\mathbb{D}({\mathcal{E}}^{\mathrm{NL}})({\mathfrak{S}})$ we need some results from functional analysis, namely:

The product and quotient rules for the Frèchet derivatives, [Schwartz and Karcher, 1969]. 2. 2.

The property that the Frèchet derivative of a linear operator is the operator itself, [Schwartz and Karcher, 1969]. 3. 3.

Let $A\to A^{r}$ be the map that takes a positive definite matrix to its $r$-th power, and let $\mathbb{D}A^{r}$ be the Frèchet derivative of this map. Then it is known that $\|\mathbb{D}A^{r}\|=\|rA^{r-1}\|$ if $r\leq 1$ or $r\geq 2$, [Bhatia and Holbrook, 2003]. 4. 4.

Triangle inequality.

Let us denote the non-linear part of $\tilde{\mathcal{L}}$ as $\tilde{\mathcal{L}}^{\mathrm{NL}}$, and the linear part - as $\tilde{\mathcal{L}}^{L}$. Accordingly, $\tilde{\mathcal{L}}=\tilde{\mathcal{L}}^{\mathrm{NL}}+\tilde{\mathcal{L}}^{L}$. Using the properties of the Frèchet derivatives, we obtain

[TABLE]

It can be seen that the first norm $\|\cdot\|_{1}$ in the right hands part of Eq.(B.1) is bounded, and moreover, $\|\cdot\|_{1}\leq 2$. The second norm $\|\cdot\|_{2}$ is also bounded since we work on a truncated grid with $K<K_{\max}$. Therefore, the values of ${\mathfrak{S}}$ are also bounded. For instance, in our numerical experiments ${\mathfrak{S}}<8$. For the third norm $\|\triangledown_{u}^{2}\|$ it is well-known, [Starzak, 1989], that at the uniform grid

[TABLE]

where $N_{u}$ is the size of the matrix $A_{2,u}^{C}$, and $h_{u}$ is the grid step. For the uniform grid $h_{u}=1/(N_{u}-1)$. Therefore,

[TABLE]

where $\|\Sigma\|$ is the norm of matrix $A(\Sigma)$ on the grid.

A similar analysis can be provided for the norm $\left\|\mathbb{D}({\Delta T}L^{\mathrm{NL}}_{v})({\mathfrak{S}})\right\|$ to obtain

[TABLE]

since

[TABLE]

Finally for the mixed derivative term we have

[TABLE]

Similarly to Eq.(B.2),

[TABLE]

and also

[TABLE]

Therefore,

[TABLE]

Finally, using these expressions, and properties of the Frèchet derivatives, we obtain

[TABLE]

where $\|L^{L}\|$ is the norm of the linear part of the operator $\tilde{\mathcal{L}}_{u}^{(k-1)}+\tilde{\mathcal{L}}_{v}^{(k-1)}+\tilde{\mathcal{L}}_{uv}^{(k-1)}$ in Eq.(60). It is easy to show that

[TABLE]

where $|\alpha^{L}|\leq 1$.

As we mentioned already, on the FD grid $\|{\mathcal{E}}\|$ is a spectral norm of the matrix exponential $e^{\mathcal{M}}$ where

[TABLE]

Here $\alpha^{\mathrm{NL}},\alpha^{L}$ are some terms that don’t depend on ${\Delta T},h_{u},h_{v}$, i.e. on the FD steps in all spatial and temporal directions. Also, by construction $M=A({\mathcal{M}})$ is the Metzler matrix with all eigenvalues being negative. Therefore, it can be seen that by decreasing the steps $h_{u},h_{v}$, or increasing the step ${\Delta T}$ we decrease the norm $\|\mathbb{D}({\mathcal{E}})({\mathfrak{S}})\|$ (as the exponential term dominates). Therefore, to decrease the norm $\|\mathbb{D}({\mathcal{E}})({\mathfrak{S}})\|$ we can, e.g., fix the time step ${\Delta T}$, and then decrease the step $h_{u}$ or $h_{v}$ (or both) until we get $\|\mathbb{D}({\mathcal{E}})({\mathfrak{S}})\|<1$.

This could be compared with the familiar stability condition for the Euler explicit FD scheme, which reads, [Itkin, 2017b]

[TABLE]

with $D$ being some diffusion coefficient. The latter condition is restrictive, because it sets the upper limit on the time step given the space step $h$. In contrast, our condition, using a loose notation, reads

[TABLE]

Since this is a stability condition, the convergence of our iterative method is conditional. However, our condition is far less restrictive then that in Eq.(B.9). Indeed, given the space step $h$ we need the time step to exceed the value $h^{2}/(2D)$, while Eq.(B.9) requires it to be less than $h^{2}/(2D)$. Since $h^{2}$ is normally small, and $D$ is proportional to $\|\Sigma\|^{2}>1$ on the grid, it could be satisfied with no serious restrictions on ${\Delta T}$. Therefore, the iterative scheme is "almost" unconditionally stable, where "almost" means - at the values of the model parameters reasonable from the practical point of view.