An Expanded Local Variance Gamma model
Peter Carr Andrey Itkin







petercarr@nyu.edu


aitkin@nyu.edu

Tandon School of Engineering, New York University, 

12 Metro Tech Center, RH 517E, Brooklyn NY 11201, USA


The paper proposes an expanded version of the Local Variance Gamma model of Carr and Nadtochiy by adding drift to the governing underlying process. Still in this new model it is possible to derive an ordinary differential equation for the option price which plays a role of Dupire's equation for the standard local volatility model. It is shown how calibration of multiple smiles (the whole local volatility surface) can be done in such a case. Further, assuming the local variance to be a piecewise linear function of strike and piecewise constant function of time this ODE is solved in closed form in terms of  Confluent hypergeometric functions. Calibration of the model to market smiles does not require solving any optimization problem and, in contrast, can be done term-by-term by solving a system of non-linear algebraic equations for each maturity, which is faster.



local volatility, stochastic clock, Gamma distribution, piecewise linear variance, Variance Gamma process, closed form solution, fast calibration, no-arbitrage.







§ INTRODUCTION


Local volatility model was introduced by <cit.> and <cit.> as a natural extension of the celebrating Black-Scholes model to take into account an existence of option smile. It is able to exactly replicate the local volatility function σ(T,K) where K,T are the option strike and time to maturity, at given pairs (T,K)  where the European options prices or their implied volatilities are known. This process is called calibration of the local volatility (or, alternatively, implied volatility) surface. Various approaches to solving this important problem were proposed, see, eg, survey in  <cit.> and references therein.

As mentioned in <cit.>, there are two main approaches to solving the calibration problem. The first approach relies on some parametric or non-parametric regression to construct a continuous implied volatility (IV) surface matching the given market quotes. Then the corresponding local volatility surface can be found via the well-known Dupire's formula, see, e.g., <cit.> and references therein.

The second approach relies on the direct solution of the Dupire equation using either analytical or numerical methods. The advantage of this approach is that it guarantees no-arbitrage [But only if an analytical or numerical method in use does preserve no-arbitrage. This includes various interpolations, etc.]. However, the problem of the direct solution can be ill-posed, <cit.>, and is rather computationally expensive. For instance, in <cit.> the Dupire equation (a partial differential equation (PDE) of the parabolic type) is solved by i) first using the Laplace-Carson transform, and ii) then applying various transformations to obtain a closed form solution of the transformed equation in terms of Kummer hypergeometric functions. Still, it requires an inverse Laplace transform to obtain the final solution.

With the second approach in use one also has to make an assumption about the behavior of the local/implied volatility surface at strikes and maturities where the market quotes are not known. Usually, by a tractability argument the corresponding local variance is seen either piecewise constant, <cit.>, or piecewise linear <cit.> in the log-strike space, and piecewise constant in the time to maturity space [See, however, comments in <cit.> about their assumptions.].

To improve computational efficiency of calibration, an important step is made in <cit.> where Local Variance Gamma (LVG) model has been introduced (the first version refers to 2014 and can be found in <cit.>).  This model assumes that the risk-neutral process for the underlying futures price is a pure jump Markov martingale, and that European option prices are given at a continuum of strikes and at one or more maturities. The authors construct a time-homogeneous process which meets a single smile and a piecewise time-homogeneous process, which can meet multiple smiles. However, in contrast to eg, <cit.>, their construction leads not to a PDE, but to a partial differential difference equation (PDDE), which permits both explicit calibration and fast numerical valuation. In particular, it does not require application of any optimization methods, rather just a root solver. In <cit.> this model is used to calibrate the local volatility surface assuming its piecewise constant structure in the strike space.

One of the potential criticism of this calibration method is the fact that the resulting local volatility function has a finite number of discontinuities. So it would be advantaged to relax the piecewise constant behavior of the surface. This is similar to how <cit.> was developed to overcome the same problem as compared with <cit.>.

On this way, recently <cit.> applied the LVG model to the FX options market where usually option prices are quoted only at five strikes. They assumed that the local volatility function is continuous, piecewise linear in the four inner strike subintervals and constant in the outer subintervals. A closed form solution of the PDDE derived in <cit.>) is obtained with this parametrization, and calibration of some volatility smiles is provided. Still, to calibrate the model the authors rely on a residual minimization by using a least-square approach. So, despite an improved version of the LVG model is used, computational efficiency of this method is not perfect.

Another remark of <cit.> is about the limitation that the risk-neutral price process of the underlying is assumed to be a martingale, i.e. the main driving process in (<ref>) doesn't have a drift. However, the drift may not be negligible. If the drift is deterministic, e.g when the interest rate and dividends are deterministic, and the drift is a deterministic function of them,  the calibration problem can be reduced to the driftless case by discounting, but this assumption might be inconsistent with the market. Therefore, an expansion of the proposed model that allows for a non-zero and stochastic drift is very desirable. In particular, it would be interesting to expand the LVG model to a risk-neutral price process obtained by stochastic time change of a drifted diffusion. In this way, similar to local Variance Gamma model, <cit.>, we introduce both stochastic volatility and stochastic drift.

With this in mind, our ultimate goals in this paper are as follows. First, we propose an expanded version of the LVG model by adding drift to the governing underlying process. It turns out that to proceed we need to re-derive and re-think every step in construction proposed in <cit.>. We show that still it is possible to find an ordinary differential equation (ODE) for the option price which plays a role of Dupire's equation for the standard local volatility model, and how calibration of multiple smiles (the whole local volatility surface) can be done in such a case.

Further, assuming the local variance to be a piecewise linear function of strike and piecewise constant function of time we solve this ODE in closed form in terms of Confluent hypergeometric functions. Calibration of the model to market smiles does not require solving any optimization problem. In contrast, it can be done term-by-term by solving a system of non-linear algebraic equations for each maturity, and thus is much faster.

The rest of the paper is organized as follows. In Section <ref> the Expanded Local Variance Gamma model is formulated. In Section <ref> we derive a forward equation (which is an ordinary differential equation (ODE)) for Put option prices using a homogeneous Bochner subordination approach. Section <ref> generalizes this approach by considering the local variance being piece-wise constant in time. In Section <ref> a closed form solution of the derived ODE is given in terms of Confluent hypergeometric functions. The next Section discusses computation of a source term of this ODE which requires a no-arbitrage interpolation. Using the idea of <cit.>), we show how to construct non-linear interpolation which provides both no-arbitrage, and a nice tractable representation of the source term, so that all integrals in the source term can be computed in closed form. In Section <ref> calibration of multiple smiles in our model is discussed in detail. To calibrate a single smile we derive a system of nonlinear algebraic equations for the model parameters, and explain how to obtain a smart guess for their initial values.  In Section <ref>
asymptotic solutions of our ODE at extreme values of the model parameters are derived which improve computational accuracy and speed of the numerical solution. Section <ref> presents the results of some numerical experiments where calibration of the model to the given market smiles is done term-by-term. The last Section concludes.




§ PROCESS
 

Below where possible we follow the notation of <cit.>.

Let W_t be a ℚ standard Brownian motion with time index t ≥ 0. Consider a stochastic process D_t to be a time-homogeneous diffusion

    d D_t = μ D_t d t + σ(D_t) d W_t,

where the volatility function σ is local and time-homogeneous, and  μ is deterministic.

A unique solution to (<ref>) exists if σ(D) : ℝ→ℝ is
Lipschitz continuous in D and satisfies growth conditions at infinity. According to (<ref>) we have D_t  ∈ (-∞,∞) while t ∈ [0,∞). Since D is a time-homogeneous Markov process, its infinitesimal generator  is given by

    ϕ(D) ≡[μ D _D  + 1/2σ^2(D) ^2_D ] ϕ(D)

for all twice differentiable functions ϕ. Here _x is a first order
differential operator on x. The semigroup of the D process is

    ^D_t ϕ(D_t) = e^t ϕ(D_t)  = [ϕ(D_t)|D_0 = D],   ∀ t ≥ 0.


In the spirit of Variance Gamma model, <cit.> and similar to <cit.>, introduce a new process D_Γ_t which is D_t subordinated by the unbiased Gamma clock Γ_t. The density of the unbiased Gamma clock Γ_t at time t ≥ 0 is

    ℚ{Γ_t ∈ dν} = ν^m-1 e^-ν m /t(t^*)^m Γ(m) dν,   ν > 0,    m ≡ t/t^*.

Here t^* > 0 is a free parameter of the process, Γ(x) is the Gamma function. It is easy to check that

    [Γ_t] = t.

Thus, on average the stochastic gamma clock Γ_t runs synchronously with the calendar time t.

As applied to the option pricing problem, we introduce a more complex construction.
Namely, consider options written on the underlying process S_t. Without loss of generality and for the sake of clearness let us treat below S_t as the stock price process. Here, in contrast to <cit.>, we don't ignore interest rates r and continuous dividends q assuming them to be deterministic (below for simplicity of presentation we treat them as constants, but this can be easily relaxed). Then, let us define S_t as

    S_t = D_Γ_X(t),      X(t) = 1 - e^-(r-q) tr-q.

It is clear that in the limit r → 0,  q → 0 we have X(t) = t, i.e., in this limit our construction coincides with that in <cit.> who considered a driftless diffusion and assumed S_t = D_Γ_t. Also based on (<ref>)

    [Γ_X(t)] = X(t).

Function X(t) starts at zero, ie, X(0) = 0, and is a continuous increasing function of time t. Indeed, if r - q > 0, then X(t) is increasing in t on t ∈ [0,∞), and at t →∞ it tends to constant. The infinite time horizon is not practically important, but for any finite time t function X(t) can be treated as an increasing function in t. If r - q < 0, function X(t) is strictly increasing ∀ t ∈ [0,∞). Thus, X(t) has all properties of a good clock. Accordingly, Γ_X(t) has all properties of a random time.

Under a risk-neutral measure ℚ, the total gain process, including the underlying price appreciation and dividends, after discounting at the risk free rate should be a martingale, see, eg, <cit.>. This process obeys the following stochastic differential equation

    d ( e^-r t S_t e^q t)    =  e^(q-r)t[ (q-r) S_t dt + d S_t].

Taking an expectation of both parts we obtain

    [d ( e^(q-r)t S_t)]    =
    e^(q-r)t{ (q-r)  [S_t] d t + d [S_t] }.

Observe, that from (<ref>), (<ref>)

    [d S_t]    = [d ] = μ [ d ] + [σ(D_Γ_X(t)) d W_Γ_X(t)] = μ [ d ],

because the process W_Γ_t is a local martingale, see <cit.>, chapter 6. Accordingly, the process W_Γ_X(t) inherits this property from W_Γ_t, hence [σ(D_Γ_X(t)) d W_Γ_X(t)] = 0.

Further assume that the Gamma process Γ_t is independent of W_t (and, accordingly,  is independent of W_). Then the expectation in the RHS of (<ref>) can be computed, by first conditioning on , and then integrating over the distribution of  which can be obtained from (<ref>) by replacing t with X(t), i.e.

    [ d   | S_s]    = ∫_0^∞ [ d |  = ν] ν^m-1 e^-ν m /X(t)(t^*)^m Γ(m)
       = ∫_0^∞ [D_ν] ν^m-1 e^-ν m /X(t)(t^*)^m Γ(m) d ν,   ν > 0,    m ≡ X(t)/t^*.

The find [D_ν] we take into account (<ref>) to obtain

    d  [D_ν] =   [d D_ν]
    =  [μ D_ν dν + σ(D_ν) D_ν d W_ν] = μ [D_ν] dν.

Solving this equation with respect to y(ν) =  [D_ν  | D_s],  we obtain [D_ν  | D_s] = D_s e^μ (ν-s). Since we condition on time s, it means that D_s = D_Γ_X(s) = S_s, and thus
[D_ν  | D_s] = S_s e^μ (ν-s).

Further, we substitute this into (<ref>), set the parameter of the Gamma distribution t^* to be t^* = X(t) (so m = 1) and integrate to obtain

    d [S_t  | S_s] = [d S_t  | S_s] = μ [ d ] = S_s e^- s μμ/1 - μ X(t).

Setting now m = r-q and solving this equation we find

    [S_t  | S_s] = S_s (r-q) e^(q-r)(s-t).

Substituting <ref> and <ref> into <ref> yields d ( e^-r t S_t e^q t) = 0. Thus, if we chose μ = r-q, the right hands part of (<ref>) vanishes, and our discounted stock process with allowance for non-zero interest rates and continuous dividends becomes a martingale. So the proposed construction can be used for option pricing.

This setting can be easily generalized for time-dependent interest rates r(t) and continuous dividends q(t). We leave it for the reader.

The next step is to consider connection between the original and time-changed processes. It is known from <cit.>  that the process G_Γ_t defined as

    d G_t = σ^2(G) d W_t

is a time-homogeneous Markov process. As the deterministic process μ t is also time-homogeneous, the whole process D_t defined in (<ref>) is also a time-homogeneous Markov process. Accordingly, the semigroups T^S_t of S_t and T^D_t of D_Γ_X(t) are connected by the Bochner integral

    ^S_t U(S) = ∫_0^∞^D_ν U(S) ℚ{∈ dν},   ∀ t ≥ 0,

where U(S) is a function in the domain of both ^D_t and ^S_t.
It can be derived by exploiting the time homogeneity of the D process, conditioning on the gamma time first, and taking into account the independence of Γ_t and W_t (or Γ_ and W_ in our case).

We set parameter t^* of the gamma clock to t^* = X(t). Then (<ref>) and (<ref>) imply

    ^S_t U(S) = ∫_0^∞^D_ν U(S) e^-ν/X(t)X(t) dν.

In what follows for the sake of brevity we will call this model as an Expanded Local Variance Gamma model, or ELVG.



§ FORWARD EQUATION FOR PUT OPTION PRICES
 

Following <cit.> we interpret the index t of the semigroup ^S_t as the maturity date T of a European claim with the valuation time t = 0. Also let the test function U(S) be the payoff of this European claim, ie,

    U(S_T) = e^-r T(K - S_T)^+.

Then define

    P(S_0,T,K) = ^S_T U(S_0)

as the European Put value with maturity T at time t=0 in the ELVG model. Similarly

    P^D(S_0,ν,K) = ^D_ν U(S_0)

would be the European Put value with maturity ν at time t=0 in the model of (<ref>)[Below for simplicity of notation we drop the subscript '0'  in S_0.]. Then the Bochner integral in (<ref>) takes the form

    P(S,T,K) =  ∫_0^∞ P^D(S,ν,K) p e^- p ν d ν,     p ≡ 1/X(T).

Thus, P(S,X(T),K) is represented by a Laplace-Carson transform of P^D(S,ν,K) with p being a parameter of the transform. Note that

    P(S,0,K) = P^D(S,0,K) = U(S).

To proceed, we need an analog of the Dupire forward PDE for P^D(S,ν,K).



 §.§ Derivation of the Dupire forward PDE 


Despite this can be done in many different ways, below for the sake of compatibility we do it in the spirit of <cit.>.

First, differentiating (<ref>) by ν with allowance for (<ref>) yields

    _ν P^D(S,ν,K)    = e^-r ν e^ν[  - r] U(S)
    = e^-r ν[  - r] U(S).

We take into account the definition of the generator  in (<ref>), and also remind that at t=0 we have D_0 = S_0. Then (<ref>) transforms to

    _ν P^D(S,ν,K) = -r P^D(S,ν,K) + (r-q) S _S P^D(S,ν,K) + e^-r ν1/2[ σ^2(S) _S^2 U(S) ].

However, we need to express the forward equation using a pair of independent variables (ν,K) while (<ref>) is derived in terms of (ν,S). To do this, observe that

    e^-r ν[σ^2(S) _S^2 U(S) ]    = e^-r ν[σ^2(S) δ(K-S)] = e^-r ν[σ^2(K) δ(K-S)] 
       =  e^-r ν[σ^2(K) _K^2 U(S)] = σ^2(K) _K^2 P^D(S,ν,K).

where the sifting property of the Dirac delta function δ(S-K) has been used. Also

    -r     P^D(S,ν,K) + (r-q) S ∇_S P^D(S,ν,K) 
       = e^-r ν[ -r (K-S)^+ + (r-q)S (K-S)^+S] 
       = e^-r ν[ -r (K-S)^+ - (r-q)(K-S) (K-S)^+S + (r-q)K (K-S)^+S] 
       = e^-r ν[ -r (K-S)^+ + (r-q)(K-S)^+  - (r-q)K (K-S)^+K] 
       = -q P^D(S,ν,K) - (r-q) K ∇_K P^D(S,ν,K).

Therefore, using (<ref>) and (<ref>), (<ref>) could be transformed to

    ∇_ν P^D(S,ν,K)    = -q P^D(S,ν,K) - (r-q) K ∇_K P^D(S,ν,K) + 1/2σ^2(K) K^2 ∇_K^2 P^D(S,ν,K) 
       ≡^K P^D(S,ν,K),  
    ^K    = -q  - (r-q) K ∇_K  + 1/2σ^2(K) K^2 ∇_K^2.

This equation looks exactly like the Dupire equation with non-zero interest rates and continuous dividends, see, eg, <cit.> and references therein. Note, that ^K is also a time-homogeneous generator.



 §.§ Forward partial divided-difference equation
 

Our final step is to apply the linear differential operator ^K defined in  (<ref>) to both parts of (<ref>). Using time-homogeneity of D_t and, again, the Dupire equation (<ref>), we obtain

    -q    P(S,T,K) - (r-q) K _K P(S,T,K) + 1/2σ^2(K) _K^2
    P(S,T,K) 
       = ∫_0^∞ p e^- p ν[-q P^D(S,ν,K) - (r-q) K _K P^D(S,ν,K) + 1/2σ^2(K) _K^2 P^D(S,ν,K)] d ν
       = ∫_0^∞ p e^- p ν_ν P^D(S,ν,K) d ν  = - p P^D(S,0,K) + p ∫_0^∞ P^D(S,ν,K) p e^- p ν d ν
       = p[P(S,T,K) - P^D(S,0,K) ] = p[P(S,T,K) - P(S,0,K) ],

where in the last line (<ref>) was taken into account.

Thus, finally P(S,T,K) solves the following problem

    -q P(S,T,K)    - (r-q) K _K P(S,T,K) + 1/2σ^2(K) _K^2
    P(S,T,K) 
       = P(S,T,K) - P(S,0,K)X(T),      P(S,0,K) = (K-S)^+.

At r = q = 0 this equation translates to the corresponding equation in <cit.>. In contrast to the Dupire equation which belongs to the class of PDE, (<ref>) is an ODE, or, more precisely, a partial divided-difference equation (PDDE), since the derivative in time in the right hands part is now replaced by a divided difference. In the form of an ODE it reads

    [1/2σ^2(K) _K^2 - (r-q) K _K - (q + 1X(T)) ] P(S,T,K) =
    - P(S,0,K)X(T).

This equation could be solved analytically for some particular forms of the local volatility function σ(K) which are considered later in this paper. Also in the same way a similar equation could be derived for the Call option price C_0(S,T,K) which reads

    [1/2σ^2(K) _K^2 + (r-q) K _K    - (q + 1X(T)) ] C_0(S,T,K) = - C_0(S,0,K)X(T), 
    
    C_0(S,0,K)    = (S-K)^+.


Solving (<ref>) or (<ref>) provides the way to determine σ(K) given market quotes of Call and Put options with maturity T. However, this allows calibration of just a single term. Calibration of the whole local volatility surface, in principle, could be done term-by-term (because of the time-homogeneity assumption) if (<ref>), (<ref>) could be generalized to this case. We consider this in the following Section.



§ LOCAL VARIANCE PIECE-WISE CONSTANT IN TIME
 

To address calibration of multiple smiles, we need to relax some assumptions about time-homogeneity of the process D_t defined in (<ref>). This includes several steps which are described below in more detail.



 §.§ Local variance


Here we assume that the local variance σ(D_t) is no more time-homogeneous, but a piece-wise constant function of time σ(D_t,t).

Let T_1, T_2, …, T_M be the time points at which the variance rate σ^2(D_t) jumps deterministically. In other words, at the interval t ∈ [T_0,T_1), the variance rate is σ^2_0(D_t), at t ∈ [T_1, T_2) it is σ^2_1(D_t), etc. This can be also represented as

    σ^2(D_t,t)    = ∑_i=0^M σ^2_i(D_t) w_i(t), 
    
    w_i(t)    ≡1_t - T_i - 1_t - T_i+1,  i=0,…,M,
       T_0 = 0,  T_M+1 = ∞,   1_x =
    
    1,     x ≥ 0 
    
    0,     x < 0.


Note, that

    ∑_i=0^M w_i(t) = 1_t - 1_t-∞ = 1,   ∀ t ≥ 0.

Therefore, in case when all σ^2_i(D_t) are equal, ie, independent on index i, (<ref>) reduces to the case considered in the previous Sections.

It is important to notice that our construction implies that the volatility σ(D_t) jumps as a function of time at the calendar times T_0, T_1,…,T_M, and not at the business times ν determined by the gamma clock. Otherwise, the volatility function would change at random (business) times which means it is stochastic. But this definitely lies out of scope of our model. Therefore, we need to change (<ref>) to

    σ^2(, )    = ∑_i=0^M σ^2_i(D_t) w̅_i(()), 
    w̅_i(t)    = 1_X^-1(t) - T_i - 1_X^-1(t) - T_i+1,  i=0,…,M, 
    
    X^-1(t)    = 1q-rlog[1 - (r-q)t ].


Hence, when using (<ref>) we have

    σ^2(D_t, t)|_t = Γ_X(t)   = ∑_i=0^M σ^2_i(D_t) w̅_i(X(t)) = ∑_i=0^M σ^2_i(D_t) w_i(t).


Accordingly, if the calendar time t belongs to the interval T_0 ≤ t < T_1, the infinitesimal generator  of the semigroup ^D_ν is a function of σ(D_t) (and not on σ(D_ν)). As at T_0 ≤ t < T_1 we assume σ(D) = σ_0(D), i.e., is constant in time, it doesn't depend of ν. Thus,  (which for this interval of time we will denote as _0) is still time-homogeneous.

Similarly, one can see, that for T_1 ≤ t < T_2 the infinitesimal generator _1 of the semigroup ^D_ν is also time-homogeneous and depends on σ_1(D), etc.



 §.§ Bochner subordination


We start with re-definition of (<ref>), (<ref>). We now define the European Put value with maturity T at the evaluation time t=X(T_1) in the ELVG model

    P(S_0,T_1 + T,K) = ^S_T [e^-r T P(S_0,T_1,K)].

And, clearly we are interesting in the value of T to be T = T_2 - T_1.

Similarly, we define the European Put value with maturity ν at the evaluation time t=T_1 in the model given by (<ref>) as

    P^D(S_0,T_1 + ν,K) = ^D_ν [e^-r ν P(S_0,T_1,K)].

By these definitions

    P(S_0,T_1 + T,K)|_T=0 = P^D(S_0,T_1 + ν,K)|_ν=0 = P(S_0,T_1,K)].


In contrast to (<ref>), in case of multiple smiles at t > T_1 we need to change the definition of t in (<ref>) from t ↦ X(t) to

    t ↦ X(T_1+t) - X(T_1) ≡Δ x(T_1,t).

This definition implies two observations.

First, function Δ x(T_1,t) starts at zero at t=0 and is an increasing function of time. Also, in case r=q=0 we have Δ x(T_1,t) = t. Therefore, Δ x(T_1,t) can be used as a good clock. Accordingly, similar to (<ref>) we have

    [Γ_Δ x(T_1,t)] = Δ x(T_1,t).


Second, a proof that in our model the discounted stock price is a martingale
given in Section <ref> could be repeated for times t:  T_1 < t ≤ T_2.
When doing so, at t > T_1 we reset the definition of S_t to

    S_T_1+t = D_Γ_Δ x(T_1,t) ,    t ≥ 0.


Then instead of (<ref>) we now have

    [d S_T_1+t]    = [d D_Γ_Δ x(T_1,t)] = μ [D_Γ_Δ x(T_1,t) d Γ_Δ x(T_1,t)] + [σ(D_Γ_Δ x(T_1,t)) d W_Γ_Δ x(T_1,t) ] 
       = μ [D_Γ_Δ x(T_1,t)] d Δ x(T_1,t) =
    μ [D_Γ_Δ x(T_1,t)] d X(T_1+t).

On the other hand,

    [d ( e^(q-r)(T_1+t) S_T_1+t)]    =
    e^(q-r)(T_1+t){ (q-r)  [S_T_1+t] d t + d [S_T_1+t] }
       = e^(q-r)t[μ + (q-r)S_T_1e^-(r-q)T_1] dt


One can check, that with μ = r-q the RHS of <ref> vanishes, therefore this construction can be used for option pricing.

The definition in (<ref>) implies that parameter t of the Gamma random clock is reset at the point T_1, i.e., at 0 ≤ t ≤ T_1 it is t ↦ X(t) = X(t) - X(0), while at T_1 < t ≤ T_2 it is t ↦ X(T_1+t) - X(T_1). Using the definition of w_i(t) in (<ref>), this could be written as

    t ↦∑_i=0^M w_i(T_i+t)[X(T_i+t) - X(T_i)]

Resetting t was also first proposed in <cit.> but in a different form.

Then, the Bochner integral in (<ref>) transforms to

    ^S_T P(S,T_1,K) = ∫_0^∞^D_ν P(S,T_1 + ν,K) ν^m-1 e^-ν m /Δ X(T_1,T)(t^*)^m Γ(m) dν.

Since for a tractability reason we still want to have m ≡Δ X(T_1,T)/t^* = 1. we need to redefine t^* in accordance with (<ref>). Based on that, the Bochner integral in (<ref>) now finally reads

    P(S,T_1 + T,K)    = ∫_0^∞ P^D(S,T_1+ν,K) p e^- p ν d ν,    
    p ≡ 1/Δ X(T_1, T).





 §.§ Forward partial divided-difference equation for the second term
 

Now we need to derive a Forward partial divided-difference equation for the second term T_2 similar to how this is done in Section <ref>. Obviously, the Put price P^D(S_0,T_1 + ν,K) solves the same Dupire equation (<ref>). Therefore, proceeding in the same way as in Section <ref>, we apply linear differential operator L defined in  (<ref>) to both parts of (<ref>). Using time-homogeneity of D_t at the interval [T_1, T_2) and again the Dupire equation (<ref>), we obtain

    -q    P(S,T_1 + T,K) - (r-q) K _K P(S,T_1 + T,K) + 1/2σ^2(K) _K^2 P(S,T_1 + T,K) 
       = ∫_0^∞ p e^- p ν[-q P^D(S,T_1 + ν,K) - (r-q) K _K P^D(S,T_1 + ν,K)  
       + 1/2σ^2(K) _K^2 P^D(S,T_1 + ν,K)] d ν
    = ∫_0^∞ p e^- p ν_ν P^D(S,T_1 + ν,K) d ν
       = - p P^D(S,T_1,K) + p ∫_0^∞ P^D(S,T_1+ν,K) p e^- p ν dν
       = p[P(S,T_1 + T,K) - P^D(S,T_1,K)]  = p[P(S,T_1 + T,K) - P(S,T_1,K) ].

Finally, taking T = T_2  - T_1 we obtain an ODE for the Put price P(S,T_2,K).


    [1/2σ^2(K) _K^2 - (r-q) K _K - (q + 1
    X(T_2) - X(T_1)) ] P(S,T_2,K) = - P(S,T_1,K)X(T_2) - X(T_1).

Here the local variance function σ^2(K) = σ^2_1(K) as it corresponds to the interval (T_1,T_2] where the above ODE is solved.

We continue in the same way to derive an ODE for the Put price P(S,T_i,K),  i=1,…,M, which finally reads

    [1/2σ^2(K) _K^2 - (r-q) K _K - (q + 1
    X(T_i) - X(T_i-1)) ] P(S,T_i,K) = - P(S,T_i-1,K)X(T_i) - X(T_i-1).

This is a recurrent equation that can be solved for all i=1,…,M sequentially starting with i=1 subject to some boundary conditions. The natural boundary conditions for the Put option price are, <cit.>

    [                 P(S,T_i,K) = 0,                                                          K → 0,; P(S,T_i,K) = _i K - _i S ≈_i K,                                                           K →∞, ]

where _i = e^-r T_i is the discount factor, and _i = e^-q T_i.

A similar equation can be obtained for the Call option prices, which reads

    [1/2σ^2(K) _K^2 + (r-q) K _K - (q + 1
    X(T_i) - X(T_i-1)) ] C(S,T_i,K) = - C(S,T_i-1,K)X(T_i) - X(T_i-1),

subject to the boundary conditions

    [ C(S,T_i,K) = _i S,                                K → 0,;    C(S,T_i,K) = 0,                                 K →∞. ]







§ SOLUTION OF THE ODE (<REF>)
 

Below we use the approach similar to <cit.> by assuming the local variance to be a piecewise linear continuous function of strike. In contrast to <cit.>, instead of a standard local volatility model in this paper we use the ELVG model. As the result, instead of a partial differential (Dupire) equation, we face a problem of solving the ODE in (<ref>).

First, it is useful to change the dependent variable from P(S,T_j,K) to

    V(S,T_j,K) = P(S,T_j,K) - _j K,

which is known as a covered Put. The advantage of the covered Put is that according to (<ref>) its price obeys homogeneous boundary conditions.

Using this definition we now re-write (<ref>) in a more convenient form (while with some loose of notation)

    - v(x) V_x,x(x) + b_1 x V_x(x) + b_0,j V(x) = c_j, 
    
    b_1    = (r - q) p_j,    b_0,j = q p_j + 1,    c_j = V(T_j-1,x) + β x, 
    
    p_j    = X(T_j) - X(T_j-1) > 0,    x = K/S,    V(x) = V(S,T_j,x),    v(x) = p_jσ^2(x)/2 S^2. 
    β   = - S[_j(1 + p_j r) - _j-1].


In (<ref>) x is the inverse moneyness. In what follows we also assume that r > q > 0, but this assumption could be easily relaxed.

Further, suppose that for each maturity T_j,  j ∈ [1,M] the market quotes are provided at a set of strikes K_i,  i=1,…,n_j where the strikes are assumed to be sorted in the increasing order. Then the corresponding continuous piecewise linear local variance function v_j(x) on the interval [x_i,x_i+1] reads

    v_j,i(x) = v^0_j,i + v^1_j,i x,

where we use the super-index 0 to denote a level v^0, and the super-index 1 to denote a slope v^1. Subindex i=0 in v^0_j,0, v^1_j,0 corresponds to the interval (0, x_1]. Since v_j(x) is continuous, we have

    v^0_j,i + v^1_j,i x_i+1 = v^0_j,i+1 + v^1_j,i+1 x_i+1,    i=0,…,n_j-1.

The first derivative of v_j(x) experiences a jump at points x_i, i1n_j. We also assume that v(x,T) is a piecewise constant function of time, i.e., v^0_j,i, v^1_j,i don't depend on T on the intervals [T_j, T_j+1),  j ∈ [0,M-1], and jump to new values at the points T_j, j1M.

With the above assumptions in mind, (<ref>) can be solved by induction. One starts with T_0 = 0, and on each time interval [T_j-1,T_j],  j1M solves the problem (<ref>) for V(x) ↦ P(S,T_j,x) - d_j S x.

Since v(x) is a piecewise linear function, the solution of (<ref>) can also be constructed separately for each interval [x_i-1,x_i]. By taking into account the explicit representation of v(x) in (<ref>), from (<ref>) for the i-th spatial interval we obtain

    -(b_2 + a_2 x) V_xx(x)    + b_1 x V_x(x) + b_0 V(x) = c 
    
    b_2    = v^0_j,i,  a_2 = v^1_j,i.

We proceed by introducing a new independent variable z = (b_2 + a_2 x)b_1/a_2^2,  z ∈ℝ^+, so that (<ref>) transforms to

    -z V_zz(z)    + (z-q_2)V_z(z) + q_1 V(z) = χ
    
    q_1    = b_0/b_1,  q_2 = b_2 b_1/a_2^2,  χ = c/b_1.


The (<ref>) is an inhomogeneous Laplace equation, <cit.>, page 155. It is well known that if y_1=y_1(z), y_2=y_2(z) are two fundamental solutions of the corresponding homogeneous equation, then the general solution of (<ref>) can be represented as

    V(z)    = C_1 y_1(z) + C_2 y_2(z) + 1/b_1 I_12(z) 
    
    I_12(z)    = -y_2(z) ∫ y_1(z) f(z)W zd z + y_1(z) ∫ y_2(z) f(z)W zd z ≡ I_1 + I_2, 
    
    f(z)    = V(T_j-1,z) - k_1 - k_2 z,   
    k_1 = βb_2/a_2,    k_2 = - βa_2/b_1,

where W = y_1 (y_2)_z - y_2 (y_1)_z is the so-called Wronskian, and β is defined in (<ref>). Then the problem is to determine suitable fundamental solutions of the homogeneous Laplace equations. Based on <cit.>, if a_2  0, they read

    y_i(z) = 𝒱_i (q_1, q_2, z),    i=1,2

Here 𝒱_i(a,b,z) is an arbitrary solution of the degenerate hypergeometric equation, i.e., Kummer's function, <cit.>. Two types of Kummer's functions are known, namely M(a,b,z) and U(a,b,z), which are Kummer's functions of the first and second kind.

It is known, that there exist several pairs of such independent solutions. Therefore, for every spatial interval in z among all possible fundamental pairs we have to determine just one which is numerically satisfactory at this interval (see <cit.> for the detailed definition of satisfactory solutions and the corresponding discussion). Since our boundary conditions are set at zero and positive infinity, we need a numerically satisfactory solution for the positive half of the real line.

Similar to <cit.>, in the vicinity of the origin we choose the numerically satisfactory pair as, <cit.>

    y_1(χ)    = M(q_1, q_2, z) = e^z M(q_2-q_1, q_2, -z), 
    
    y_2(χ)    = z^1-q_2 M(q_1 - q_2 + 1, 2-q_2, z) =
    z^1-q_2 e^z M(1-q_1, 2-q_2, -z), 
    
    W    = sin(π q_2) z^-q_2 e^z/π.


However, in the vicinity of infinity the numerically satisfactory pair is, <cit.>

    y_1(χ)    =   U(q_1, q_2, z) = z^1-q_2 U(q_1-q_2+1, 2-q_2, z), 
    
    y_2(χ)    = e^z U(q_2 - q_1, q_2, -z) = e^z z^1-q_2 U(1-q_1, 2-q_2, -z), 
    
    W    = (-1)^q_1 - q_2 e^z z^-q_2.


As two solutions J_1(q_1,q_2, z), J_2(q_1,q_2, z) are independent, (<ref>) is a general solution of (<ref>). Two constants C_1,C_2 should be determined based on the boundary conditions for the function V(z).

The boundary conditions for the ODE (<ref>) in a strike K space (or in x space) should be set at zero and infinity. Based on the usual shape of the local variance curve and its positivity, for x → 0, we expect that v^1_j,i < 0.  Similarly, for x →∞ we expect that v^1_j,i > 0. In between these two limits the local variance curve for a given maturity T_j is assumed to be continuous, but the slope of the curve could be both positive and negative, see, e.g., <cit.> and references therein. Also, by definition z = v_j,i, and Dom(z) = ℝ^+. Thus, at high strikes a_2 = v^1_j,i > 0. Therefore, the boundary conditions for (<ref>) should be set at z = b_2 (which corresponds to the boundary K=0) and at z →∞. These are the boundary conditions given in (<ref>).



§ COMPUTATION OF THE SOURCE TERM
 

Computation of the source term p I_12 in (<ref>) could be achieved in several ways. The most straightforward one is to use numerical integration as the Put price P(x,T_i-1) as a function of x is already known when we solve (<ref>) for T=T_i. However, as this is discussed in detail in <cit.>, function P(x,T_i-1) is known only for a discrete set of points in x. Therefore, some kind of interpolation is necessary to find its values at the other points.



 §.§ No-arbitrage interpolation


As shown in <cit.>, this interpolation must preserve no-arbitrage. So, for instance, a standard linear interpolation is not a good candidate, since its violates no-arbitrage conditions. Indeed, given three Put option prices P(K_1), P(K_2), P(K_3) for three strikes K_1 < K_2 < K_3, the necessary and sufficient conditions for an arbitrage-free system read, <cit.>

    P(K_3)    > 0,      P(K_2) < P(K_3), 
    
    Bs = (K_3 - K_2)P(K_1)    - (K_3 - K_1)P(K_2) + (K_2 - K_1)P(K_3) > 0.

Suppose that we want to use linear interpolation in the strike space on the interval [K_1,K_3] to find the unknown Put option price P(K_2) given the values of P(K_1), P(K_3),

    P(K_2) ≡ P_l(K_2) = P(K_1) K_3 - P(K_3) K_1K_3 - K_1 + P(K_3) - P(K_1)K_3 - K_1 K_2.

When plugging this expression into the second line of (<ref>), the left hands side of the latter vanishes, so the third no-arbitrage condition is violated.

In <cit.> it is shown that this problem could be resolved if we use linear interpolation with a modified independent variable (further on we denote it as P_F(K)),

    P(K_2)    ≡ P_F(K_2) 
       = P(K_1) f(K_3) - P(K_3) f(K_1)f(K_3) - f(K_1) + P(K_3) - P(K_1)f(K_3) - f(K_1) f(K_2),

where f(K) is a convex and increasing function in [K_1,K_3]. Indeed, if f(K) is convex, then P(K_2) = P_F(K_2) = P_l(K_2) - ε,  ε > 0 (see Fig. 2 in <cit.>). Substitution of this expression into the second line of (<ref>) gives (K_3 - K_1) ε > 0, which is true. The second condition in (<ref>) now reads

    (P(K_1) - P(K_3))(f(K_3) - f(K_2))(f(K_1) - f(K_3)) > 0,

which is also true since f(K) is an increasing function of K.

Alternatively, one can use non-linear interpolation. In <cit.>) both approaches were combined, and it was proved that the new interpolation scheme preserves no-arbitrage. Moreover, the final representation of the modified Put price (which is a dependent variable in their approach) obtains a nice tractable representation, so the integral I_12 can be computed in closed form. Here we want to exploit the same idea, thus significantly improving performance of our model as compared with the numerical integration.

Therefore, here we propose the following interpolation scheme

    P(x)    ≡ P_F(x) = γ_1 + γ_2 x^2,    x_1 ≤ x ≤ x_3, 
    γ_1    = P(x_3) x_1^2   - P(x_1) x^2_3x_1^2 - x_3^2,     γ_2 = P(x_1) - P(x_3)x_1^2 - x_3^2.

Then Proposition similar to that in <cit.> can be proved.

 
The interpolation scheme in (<ref>) is arbitrage free at the interval [K_1,K_3].


Observe, that the no-arbitrage conditions in (<ref>) are discrete versions of the conditions

    P > 0,    P_K > 0,    P_K,K > 0,

They, in turn, correspond to the conditions

    P > 0,    P_x > 0,    P_x,x > 0,

as x'(K) = 1/S > 0. By differentiating the first line of (<ref>) one can check that the proposed interpolation obeys these conditions provided that P is an increasing function of K (or x) given the values of all other parameters to be constant. For instance, this is the case for the Black-Scholes Puts.


As by definition z is a linear function of x, a similar interpolation scheme can be used in the z space, with a similar proof of no-arbitrage.



 §.§ No-arbitrage at consecutive intervals


Proposition <ref> guarantees that the proposed interpolation doesn't introduce an arbitrage into the solution if any three strikes belong to the same interval [K_1,K_3]. However, what if we consider strikes K_2, K_3, K_4 as this is schematically  depicted in Fig. <ref>.



Here at the interval [K_1,K_3] the exact solution is depicted by the red line, while our quadratic interpolation is in blue. Accordingly, the Put prices P_1,P_3 are the market quotes, so they assumed to be the exact prices with no market arbitrage. By our construction, these prices also don't have a model arbitrage. At the consecutive interval [K_3,K_5] a similar construction applies.

We have to emphasize that this graph is pure illustrative, and no-arbitrage interpolation guarantees  that P'(K) > 0, while the blue line in Fig. <ref> doesn't support this. However, if we draw an accurate picture by using the above formulae, it would be almost impossible to distinguish the red and blue lines. Therefore, we changed convexity and skew of the blue line to make the difference visible.

By Proposition <ref> given a set of strikes K_1, K_2, K_3 the price P_2 obtained by interpolation preserves no-arbitrage. The same is true for P_5 given the Put prices P_3, P_5 at strikes K_3, K_5.
Now assume that given K_1, K_3, K_5 and P_1, P_3, P_5 we want to check the no-arbitrage conditions for the set of strikes K_2, K_3, K_4. The Proposition <ref> doesn't help in this situation, so we need a special consideration of this case.

Obviously, the first and second conditions in (<ref>) are still satisfied in this case, so we need to check that the butterfly spread is positive. Unfortunately, at the moment we don't have a general analytical solution of this problem, while some particular cases can be addressed. Thus this remains an open question. However, we checked this condition numerically. In doing so we used the Black-Scholes Put prices P_1, P_3, P_5[This is done to preserve upper bounds on the Put price that P(S, K,T,r) ≤ K e^-rT, <cit.>.] and built a 2D plot of Bs which is the left-hands side of the third line in (<ref>). The results for two cases presented in Table <ref> are presented in Fig. <ref>, <ref>.



Overall, we ran a lot of tests and didn't find any case where the butterfly spread would become negative. This partly supports our no-arbitrage interpolation. More sophisticated cases where, e.g., instead of strike K_3 in the butterfly spread at strikes K_2, K_3, K_4 we use another strike K_6 such that K_1 < K_2 < K_3 < K_4 < K_6 < K_5, could be treated in a similar way. Again, our numerical tests didn't reveal any case where a butterfly spread would become negative.






 §.§ Computing the integrals in (<ref>) far from z = 0


Using the interpolation scheme proposed in above, consider the first integral in (<ref>).  To remind, we compute it at some interval z ∈ [z_i,z_i+1],  i1n_j. Picking together the solutions in (<ref>) with the interpolation scheme for P(z,T_j-1) and Wronskians in (<ref>), and substituting them into the first integral in (<ref>) we obtain

    ∫    y_2(z) f(z,T_j-1)W zd z = A [-B_0
    + B_1 M(-2-q_1, -1-q_2, -z)
    + B_2 M(-1-q_1, -q_2, -z) 
       + B_3 M(-q_1, 1-q_2, -z)], 
    
    A    = 1b_1^2 q_1π(1-q_2) (π q_2), 
    
    B_0    = 1/a_2^2 (q_1+1) (q_1+2)[
    a_2 b_1 (q_1+2) (a_2^2 β q_2 - b_1 (q_1+1) (β  b_2-a_2 γ_1)
    ) 
       + γ_2 (2 a_2^4 q_2 (q_2+1)-2 a_2^2 b_1 b_2 (q_1+2)q_2 + b_1^2 b_2^2 (q_1+1) (q_1+2) ) ],
    
    
    B_1    = 2 a_2^2 γ_2 q_2 (q_2+1)/(1+q_1)(2+q_1), 
    
    B_2    = (a_2 b_1 β - 2 b_1 b_2 γ_2 + 2 a_2^2 γ_2 z)  q_2/1+q_1, 
    
    B_3    = 1/a_2^2[a_2 b_1(a_2^2 β  z + a_2 b_1 γ_1 - β b_1 b_2 )
     + γ_2 (b_1 b_2 - a_2^2 z)^2].


Similarly

    ∫    y_1(z) f(z,T_j-1)W zd z = A̅[
    B̅_1 M(q_2-q_1, 1+q_2, -z)
    + B̅_2  22q_2-q_1, 1+q_2; q_2,2+q_2; -z
       + B̅_3  22q_2-q_1, 2+q_2; q_2,3+q_2; -z], 
    
    A    = π z^q_2(π q_2) Γ(q_2), 
    
    B_1    =
    a_2^2 γ_1 - a_2 β  b_2 + b_2^2 γ_2/a_2^2 Γ(1+q_2), 
    
    B_2    =Γ(q_2+1) /Γ(q_2) Γ(2+q_2) b_1 (a_2 β -2 b_2 γ_2) z, 
    
    B_3    = Γ(q_2+1) /Γ(q_2) Γ(3+q_2) b_1^2
    a_2^2(1+q_2)γ_2 z^2,


where pqa_1,...,a_p; b_1,...,b_q; z is the generalized hypergeometric function, <cit.>.




 §.§ Computing the integrals in (<ref>) far from z = ±∞


Here we proceed in the same way as in the previous section. Again, we pick together the solutions in (<ref>) with the interpolation scheme for P(z,T_j-1) and Wronskians in (<ref>), and substitute them into the first integral in (<ref>) we obtain

    ∫    y_2(z) f(z,T_j-1)W zd z = (-1)^q_2 - q_1 [C_0 J_0 + C_1 J_1 + C_2 J_2], 
    
    J_i    = ∫ z^i U(1-q_1, 2-q_2, -z) dz, 
    
    C_0    = b_2^2 γ_2/a_2^2 - β b_2/a_2 + γ_1,   
    C_1 = a_2 β - 2 b_2 γ_2/b_1,   
    C_2 = a_2^2 γ_2/b_1^2.


It is known, <cit.>, that

    J_0 = - 1/q_1U(-q_1, 1-q_2, -z).

Then, J_1, J_2 can be found using integration by parts to yield

    J_1    = z J_0 + 1/q_1(1+q_1)U(-1-q_1, -q_2, -z), 
    
    J_2    = z J_1 + 1/q_1(2+3 q_1 + q_1^2)U(-2-q_1, -1-q_2, -z)
     - z1/q_1(1 + q_1)U(-1-q_1,-q_2,-z).


Similarly

    ∫    y_1(z) f(z,T_j-1)W zd z = (-1)^q_2 - q_1 [C_0
    𝒥_0 + C_1 𝒥_1 + C_2 𝒥_2], 
     J_i    = ∫ z^i e^-z U(1+q_1-q_2, 2-q_2, z) dz.

The integrals 𝒥_i have been considered in <cit.> using the approach of <cit.>. Borrowing from there the result


    𝒥_0 = ∫ e^-z U(1+q_1-q_2,2-q_2,z) d z = -e^-z U(q_1-q_2,1-q_2,z),

and using integration by parts, we obtain

    𝒥_1    = z 𝒥_0 + e^-zU(q_1-q_2-1, -1-q_2, z), 
    𝒥_2    = z 𝒥_1 - ∫𝒥_1 dz = (z-1) 𝒥_1 - ∫ e^-zU(q_1-q_2-1, -1-q_2, z) dz 
       = (z-1) 𝒥_1 + e^-zU(q_1-q_2-2, -2-q_2, z).





 §.§ Some additional notes


Based on the no-arbitrage interpolation and some analytics proposed in this Section, we managed to find the solution (<ref>) of the forward equation  (<ref>) in closed form. This solution by construction is arbitrage free at any interval where the local variance function defined in (<ref>) is linear. In other words we proved, that if we consider, say 3 strikes 0 < K_1 < K_2 < K_3 < ∞ such that, e.g., x_1  = K_1/S ∈ [x_i, x_i+1],  x_2  = K_2/S ∈ [x_i, x_i+1],  x_3  = K_3/S ∈ [x_i, x_i+1], then the solution at these 3 points obeys no-arbitrage conditions.




§ CALIBRATION OF SMILE FOR A GIVEN TERM T_I
 

Calibration problem for the local volatility model can be formulated as follows: given market quotes of Call and/or Put options corresponding to a set of N strikes {K}:= K_j,  j ∈ [1,N] and same maturity T_i, find the local variance function σ(K) such that these quotes solve equations in (<ref>), (<ref>).

As mentioned in <cit.>, there are two main approaches to solving this problem. The first approach attempts to construct a continuous implied volatility (IV) surface matching the market quotes by using either some parametric or non-parametric regression, and then generates the corresponding LV surface via the well-known relationship between the local and implied variances also known as the Dupire formula, see, e.g., <cit.> and references therein. To be practically useful, this construction should guarantee no arbitrage for all strikes and maturities, which is a serious challenge for any model based on interpolation. If the no-arbitrage condition is satisfied, then the LV surface can be calculated using the Dupire formula. The second approach relies on the direct solution of the corresponding forward equation (which is the Dupire equation in the Black-Scholes world, or (<ref>), (<ref>) in our model) using either analytical or numerical methods. The advantage of this approach is that it guarantees no-arbitrage. However, the problem of the direct solution can be ill-posed, <cit.>, and is rather computationally intensive.

In this Section we show that the second approach could be significantly simplified when using the ELVG model, so calibration of the smile could be done very fast and accurate.

Further, for the sake of certainty, suppose that all known market quotes are Puts, despite this can be easily relaxed. Also, suppose that the shape of a local variance is given by some function σ_j(K) = f_j(K, p_1,…,p_L), where p_1,…,p_L is a set of the model parameters to be determined. For instance, in <cit.> the local variance is assumed to be a piecewise constant function of strike, while in <cit.> this is a piecewise linear function of strike.

In this paper we also assume the local variance to be a piecewise linear function of strike. Moreover, for our model we obtained a closed form representation of the Put option prices via parameters of the model given in Sections <ref>,<ref>. Therefore, calibration of the model to the given set of smiles could be provided as follows. First, using the above-mentioned closed form solution for a fixed interval in x where parameters of the model are constant, we construct the combined solution valid for all x ∈ℝ^+. At the second step, the parameters of the local variance function v^0_j,i, v^1_j,i can be found together with the integration constants C_1, C_2 in (<ref>) by solving a system of non-linear algebraic equations.




 §.§ The combined solution in x ∈ℝ^+ 


Suppose that the Put prices for T=T_j are known for n_j ordered strikes. The location of these strikes on the x line is schematically depicted in Fig. <ref>.



Recall, that the Put prices are given by (<ref>), which in a more convenient form at the interval x_i-1≤ x ≤ x_i and at T=T_j can be represented as

    P_i(x)    = ij J_1(q_1, q_2, z)
    + ij J_2(q_1, q_2, z) + 1/b_1I_12(z) + _j K, 
    
    z    ≡ (b_2 + a_2 x)b_1/a_2^2 = (v_j,i^0 + v_j,i^1 x)b_1/a_2^2.

Here, for consistency we change notation of two integration constants which belong to the i-th interval in x and j-th maturity to ij,ij.

For the open interval B_1 in Fig. <ref>, since function K_ν(z) diverges when  z → 0, we have to put 1j = 0 as the boundary condition[
Actually, since x → 0 implies z = v → b_2, so b_2 should be non-negative, b_2 ≥ 0. Therefore, the only case when z → 0 at x → 0 is when b_2 = 0.]. Therefore, (<ref>) contains just one yet unknown constant 1j. For the closed intervals x ∈ [x_i-1, x_i],  i ∈ [2,n_j] the solutions in (<ref>) have two yet unknown constants ij, ij, since x is finite on the corresponding intervals, and both solutions y_1(x), y_2(x) are well-behaved. Finally, for the interval x ∈ [x_n_j, ∞), according to the boundary conditions in (<ref>) we must set C^(2)_2,n_j+1 = 0.

Rigorously speaking, we also have to show that in the limits x → 0 and x →∞ the source term I_12(z) in (<ref>) also vanishes. This could be done similar to Proposition 2 in <cit.>.

Thus, we have 2 n_j unknown constants to be determined. Since the local volatility function v_i is continuous at the points x_i,  i=1,…,n_j, so should be the Put options prices P(x,T_j). Therefore, we require that at the points x_i,  i = 1,...,n_j the solution for Puts and its first derivative in x should be a continuous function of x. Thus, if the local variance function is known, the above constants solve a system of 2 n_j algebraic equations. This system has a block-diagonal structure where each block is a 2x2 matrix. Therefore, it can be easily solved with the linear complexity O(n_j).

When computing the first derivatives, we take into account that, <cit.>

    M(a,b,z)z   = abM(a+1,b+1,z),   U(a,b,z)z = -a U(a+1,b+1,z), 
    ∂_z I_12(z)    = [y'_1(z)/y_1(z) I_1 + y'_2(z)/y_2(z) I_2 ] a_2.

Therefore, computing the derivatives of the solution doesn't cause any new technical problem.



 §.§ Additional equations for calibration


As we have already mentioned above, the standard way of doing calibration of the local volatility model would be that described, e.g., in <cit.>. Namely, given the maturity T_j and some initial guess of the local variance parameters v_j,i^0, v_j,i^1,  ∀ i ∈ [1,n_j], the following steps represented in Panel <ref> have to be achieved, e.g., in the standard least-square method,








Here V^market(z_i) are the market Put quotes at the given strikes and maturity. Obviously, when the number of calibration parameters (strikes) is high, this algorithm is slow even if the closed form solution is known and can be used at Step 2. Things become even worse when a numerical solution at Step 2 has to be used if the closed form solution is not available.

However, in our case this tedious algorithm can be fully eliminated. Indeed, at every point i in strike space, i ∈ [1,n_j] we have four unknown variables v_j,i^0, v_j,i^1, ij, ij. We also have four equations which contain these variables, namely

    P_i(x)|_x=x_i   = P_i+1(x)|_x = x_i, 
    
    P_i(x)|_x=x_i   = P_market(x_i), 
    P_i+1(x)x|_x = x_i   =  P_i(x)x|_x = x_i, 
    
    v^0_j,i + v^1_j,i x_i   = v^0_j,i+1 + v^1_j,i+1 x_i,    i=1,…,n_j.

Also, based on (<ref>), the last line in (<ref>) could be re-written as a recurrent expression

    v^0_j,i = v^0_j,n_j + ∑_k=i+1^n_j x_k(v^1_j,k - v^1_j,k-1) ,    i=0,…,n_j-1.


The (<ref>) is a system of 4 n_j nonlinear equations with respect to 4 (n_j+1) variables v_j,i^0, v_j,i^1, ij, ij. We remind that according to the boundary conditions 1j = n_jj = 0. Therefore, we need two additional conditions to provide a unique solution. For instance, often traders have an intuition about the asymptotic behavior of the volatility surface at infinity, which, according to our construction, is determined by v^1_j,n_j and v^1_j,0.

Overall, solving the nonlinear system of equations (<ref>) provides the final solution of our problem. This can be done by using standard methods, and, thus, no any optimization procedure is necessary. However, a good initial guess still would be helpful for a better (and faster) convergence.



 §.§ Smart initial guess


The initial guess of the solution of (<ref>) can be constructed, for instance, as follows. We take advantage of the fact that according to (<ref>) the local variance function v(x) could be explicitly expressed as

    v(x) = b_1 x V_x(x) + b_0 V(x) - cV_x,x(x).


Given maturity T_j and approximating derivatives by central finite differences with the second order of approximation in step h in the strike space (see, e.g. <cit.>), (<ref>) can be represented in the form

    v^0_j,i + v^1_j,i x_i   = b_1 x V_x(x_i) + b_0,j V(x_i) - c_jV_x,x(x_i), 
    
    V_x(x_i)    = α_-1 V(x_i-1) + α_0 V(x_i) + α_1 V(x_i+1), , 
    
    V_x,x(x_i)    = δ_-1 V(x_i-1) + δ_0 V(x_i) + δ_1 V(x_i+1), 
    α_-1   = - h_i+1h_i(h_i+1 + h_i),   α_0 = h_i+1-h_ih_i+1h_i,   α_1 = h_ih_i+1(h_i+1 + h_i). 
    δ_-1   =  2h_i(h_i+1 + h_i),   δ_0 = -2h_i+1h_i,   δ_1 = 2h_i+1(h_i+1 + h_i). 
    
    h_i    = x_i-x_i-1,    i ∈ [1,n_j].

Further, associating Put prices P(S,T_j,x_i) with the given market quotes, the right hands side of the first line in (<ref>) can be found explicitly. This then can be combined with the last line of (<ref>) to produce a system of 2(n_j-1) equations for v^1_j,i and v^1_j,i,  i ∈ [1,n_j]. Finally, we take into account the asymptotic behavior of the volatility surface in x at zero and infinity, which, according to our construction, is determined by v^1_j,n_j and v^1_j,0 and is assumed to be known. Thus, we obtain a closed system of 2(n_j-1) linear equations with a banded matrix which can be easily solved with a linear complexity. This provides an explicit representation of the local variance function over the whole set of intervals in the strike space determined according to our approximation where the continuous derivatives are replace by finite differences.

Note, that at the first and last strike intervals the approximation of the first and second derivatives by central finite differences should be replaced by one-sided approximations, in more detail see <cit.>, chapter 2.

It could also happen that at some strikes this solution (the smart guess) gives rise to a negative local variance. In such a case we do another step which is a kind of smoothing. Namely, we exclude from the initial guess all values where the local variance is negative and using the remaining points create a spline. Then the negative values in the initial guess are replaced by those given by the constructed spline.

The final step utilizes the exact representation (<ref>) of the Put price in the ELVG model. As the variance function is already known from the previous step, this equation contains two yet unknown constants ij,ij. Accordingly, they can be found by solving the system of 2 linear equations represented by the first and third  lines of (<ref>). Then, after this last step is complete, all unknown variables are determined, and thus found solution could be used as an educated initial guess for solving (<ref>) numerically.



§ ASYMPTOTIC SOLUTIONS
 

In many practical situations either some coefficients a_2 = v^1_j,i, or both
b_2 = v^0_j,i,  a_2 = v^1_j,i in (<ref>) are small. Of course, in that case the general solution (<ref>) remains valid. However, in this case when computing the values of Kummer functions numerically, numerical errors significantly grow. This is especially pronounced when computing the integral I_12. The main point is that either the Kummer function takes a very small value, and then the constants ij,ij should be big to compensate, or vice versa. Resolution of this requires a high-precision arithmetics, and, which is more important, taking many terms in a series representation of the Kummer functions, which significantly slows down the total performance of the method.

On the other hand, to eliminate these problems we can look for asymptotic solutions of (<ref>) taking into account the existence of small parameters from the very beginning. This approach was successfully elaborated on in <cit.>, and below we proceed in a similar spirit.



 §.§ Small a_2


We can build the solution of (<ref>) directly using an independent variable x (so not switching to the variable z). We represent it as a series on the small parameter a_2, i.e.

    V(x) = ∑_i=0^∞ a_2^i V_i(x).

In the zero-order approximation by plugging (<ref>) into (<ref>) and neglecting by terms proportional to a_2 ≪ 1 we obtain the following equation for V_0(x)

    -b_2 V_xx(x) + b_1 x V_x(x) + b_0 V(x) = c.

This equation is simpler than (<ref>). Still, its solution is given by a general formula

    V(x) = C_1 y_1(x) + C_2 y_2(x) + I_12(x),

but the fundamental solutions y_1(x), y_2(x) now read

    y_1(x) = ℋ(-b_0/b_1, √(b_1/2 b_2) x),     
    y_2(x) = M(b_0/2 b_1, 1/2, b_1/2 b_2x^2),

where ℋ(a, x),  a,x ∈ℝ is the generalized Hermite polynomial H_a(x), <cit.>.



 §.§ Small |z|


Based on the definition of z = (b_2 + a_2 x)b_1/a_2^2, this could occur in two cases: either at some finite interval in the strike space |a_2| ≫ |b_1 x|,  |a_2| ≫ |b_2|, or just z is small, so b_2 and a_2 have the opposite signs.  In any case we have a small parameter under the high-order derivative. This equation belongs to the class of singularly perturbed differential equations, <cit.>. It can be solved by using either the method of matching asymptotic expansions, <cit.>, or the method of boundary functions, <cit.>. The latter was used in <cit.> in a similar situation, so for further details we refer a reader to that paper.

However, we can partly eliminate this by constructing solutions of (<ref>) using the original variable x. Then we have to consider various cases where instead of a small parameter z some other combinations of parameters could be small or large. But if so, a general solution as a function of the original independent variable x could be represented as regular series on the new small parameter. Then, truncating the series, one gets a simplified solution.

To make it more transparent let us represent the general solution of (<ref>) expressed in variable x, rather than in z, as this was done in (<ref>)

    V(x)    = ij y_1(x) + ij y_2 (x) + I_12(x), 
    
    y_i(x)    = a_2^k (b_2+a_2 x)^k 𝒱_i(
    -1 - b_0/b_1 + b_1 b_2/a_2^2, 2 - b_1 b_2/a_2^2,
    b_1/a_2^2(b_2 + a_2 x)),    i=1,2.
    
    k    = 1 - b_1 b_2/a_2^2.

Observe, that based on the definition of b_1 in <ref>, b_1 ≈ (r-q)Δ T, so usually small. Therefore, small z doesn't mean that
w is necessarily small. Below we consider two cases.



  §.§.§ w ≪ 1


As |z| ≪ 1 and w ≪ 1 we have w ≪ |a_2^2/b_1|. So a_2 ≥√(b_1).
In this case w ≪ 1 is an actual small argument. Therefore, the general solution (<ref>) can be expanded into series on small w. The condition 0 < w ≪ 1 implies that a_2 and b_2 have the opposite signs.  If a_2 > 0 (and so b_2 < 0), then in the zero-order approximation we obtain

    y_1(w)    = (a_2 w)^k-1[Γ(-k)Γ(b_0/b_1) a_2 w + O(w^2) ] - (b_1/a_2^3)^1-kΓ(k-1)/b_1Γ(k + b_0/b1)(a_2 b_1 b_2 - b_0 a_2 w) + O(w^2) , 
    
    y_2(w)    = (a_2 w)^k-1[a_2 w + O(w^2)].

As b_1 > 0 we have k-1 > 0.

If a_2 < 0 and b_2 > 0, then both RHS in (<ref>) should be multiplied by a factor
exp(-2 i π  b_1 b_2/a_2^2).



  §.§.§ a^2_2 ≫ |b_1 w|


In this case we can also expand the solution in (<ref>) into series on small z to obtain

    y_1(z)    = 1/Γ(1 + q_1 - q_2)[Γ(1 - q_2) - q_1 Γ(- q_2) z] + z^-q_2[Γ(q_2-1)/Γ(q_1)z + O(z^2)] + O(z^2), 
    
    y_2(w)    = 1 + q_1/q_2z + O(z^2).

Note, that based on the definition q_2 = b_2 b_1/a_2^2, at large a_2 the coefficient q_2 could also be small. But z/q_2 = 1 + a_2 x/b_2 = w/b_2 = O(1).




§ NUMERICAL EXPERIMENTS
 

In our numerical test we use the same data set as in <cit.>. This is done first, to compare performance and a quality of the fit for all those models. Also, we already know that these smiles are difficult to fit precisely, see discussions in <cit.>.

To remind, we take data from <http://www.optionseducation.org> on XLF traded at NYSEArca on March 25, 2014. The spot price of the index is S = 22.64, and r = 0.0148,  q=0.01. The option implied volatilities (I_call, I_put) are given in Tables <ref>,<ref>. We take all OTM quotes and some ITM quotes which are
very close to the at-the-money (ATM). When strikes for Calls and Puts coincide, we take an average of I_call and I_put with weights proportional to 1 - |Δ|_c and 1-|Δ|_p respectively, where Δ_c, Δ_p are option Call and Put deltas [By doing so we do take into account effects reported in <cit.>, who pointed out that the IVs calculated from  Call and  Put option prices corresponding to the same strike  do not coincide,  although  they  should  be  equal  in theory. Our weights are chosen according to a pure empirical rule of thumb, and a more detailed investigation of this effect is required.].

[!htb]




XLF implied volatilities for the Put options.



1|c|2[4]*T     17c|K,Put 

2-18              10        11        12        13        14        15        16        17        18        19        19.5      20        20.5      21        21.5      22        23 

4/4/2014                                                                                                                   39.53               23.77     19.73     16.67      

4/11/2014                                                                                                         35.89     30.33     26.62     22.06     18.49     16.11      

4/19/2014                                                                                               32.90               26.79               20.14               15.19     12.93 

5/17/2014                                                                           37.66     33.27     26.88               23.08               18.94               16.12     13.86 

6/21/2014                                                       40.51     37.21     31.41     27.84     23.90               21.07               18.88               16.95     15.82 

7/19/2014                                                       36.71     33.35     29.96     26.09     22.81               20.29               18.13               16.30     14.93 

12/20/2014                                             31.98     29.38     27.21     25.30     23.75     22.09               20.67               19.44               18.36     17.60 

1/17/2015     42.75     38.79     35.60     33.26     30.94     28.82     26.52     24.96     23.12     21.67               20.29               19.10               17.90     18.07 





XLF implied volatilities for the Put options.



1|c|2[4]*T     13c|K,Call 

2-14              21        21.5      22        22.5      23        23.5      24        25        26        27        28        29        30 

4/4/2014               16.60     14.69     14.40     14.86                                                                            

4/11/2014               16.89     14.96     14.52     14.77     14.98                                                                  

4/19/2014                         15.79               13.38               15.39                                                        

5/17/2014     16.71               14.48                                   13.75                                                        

6/21/2014     16.31               14.78                                   13.92     14.28     16.58                                    

7/19/2014     16.82               15.24                                   14.36     14.19     15.20                                    

12/20/2014     17.63               16.61                                   15.86     15.47     15.12     15.18     15.03                

1/17/2015                                             16.95               17.25     16.23     15.73     15.50     15.58     15.86     16.47 







We have already mentioned that in our model for each term the slopes of the smile at plus and minus infinity, v^1_j,n_j and v^1_j,0, are free parameters. So often traders have an intuition about these values. However, in our numerical experiments we take for them just some plausible values. In more detail, for a normalized variance v(x) defined in (<ref>), for all smiles we use v^1_j,0 = -0.1, and v^1_j,n_j = 0.1. Accordingly, for the instantaneous variance σ^2(x) = 2 S^2 v(x)/p_j the slopes at both zero and plus infinity are time-dependent and can be computed by using the above formula.

When calibrating the model to market data, we use the standard Matlab fsolve function, and utilize a "trust-region-dogleg" algorithm (see Matlab documentation on fsolve). Parameter "TypicalX" has to be chosen carefully to speedup calculations.

The results of this calibration which is done term-by term, are given in Fig. <ref>. Here each subplot corresponds to a single maturity T (marked in the legend) and shows market data (discrete points) and computed values (solid line). It can be seen that this simple local calibration algorithm provides a very accurate fit for all terms[Note, that in <cit.> in the last subplot the fit is not perfect in the vicinity of X = -0.5, where X = log K/F and F = S e^(r-q)T.].


We constructed the calibration algorithm to be smart enough in a sense that based on the values of parameters at each iteration it decides itself which particular solution (full or asymptotic) should be used at this iteration. We also observed that all full and asymptotic solutions are utilized by the algorithm when calibrating these market smiles.

Table <ref> presents some performance measures of our algorithm. It can be seen that the elapsed time depends on the number of iterations and function evaluations necessary to converge to the given tolerance (we use a relative tolerance ε = 10^-4). This, in turn, depends on the number of evaluated Kummer functions (for the full solution), or number of exponential and Gamma functions (for the asymptotic solutions). Of course, the asymptotic solutions are much faster to evaluate, therefore an average time to calibrate a typical term is less than a second. For the last term 8 in Tab. <ref> calibration is slow for two reasons: i) full solution is used based on the values of parameters, and 2) the number of strikes is higher than for the other terms. But the main reason is that the market data for this term is quite irregular. In any case, performance of this model is much better than that reported in both <cit.> and <cit.>.




The local variance curves obtained as a result of this fitting are given term-by-term in Fig. <ref>. The corresponding local variance surface is represented in Fig. <ref>





By comparing the surface with that given in <cit.>, one can notice that the shape is quite different while for calibration we use the same market smiles. This is because  in <cit.> the standard local volatility model is used, where the underlying price follows a Geometric Brownian motion equipped with an instantaneous local volatility function, while in this paper the model is quite different.


To look at a more regular surface, we proceed with another example which is taken from <cit.>. In that paper an implied volatility surface of S&P500 is presented, and the local volatility surface is constructed using the Dupire formula. In our test we take data for the first 12 maturities and all strikes as they are given in <cit.>, and apply our model to calibrate the local variance surface as this is described in above. When doing so we set v^1_j,0 = -0.3, and v^1_j,n_j = 0.1 for all smiles.

The results of this calibration are presented in Fig. <ref>,<ref>,<ref>. By construction, our surface preserves no-arbitrage, while for the approach in <cit.> they have to solve some additional problems[As this is mentioned in <cit.>, the correct pricing of local volatility surface requires an arbitrage free implied volatility surface. If the input implied volatility surface is not arbitrage free, this can lead to negative transition probabilities and/or negative local volatilities and can give rise to mispricing.].








In Table <ref> we present the performance of our algorithm in this experiment.
It can be seen that here the elapsed time is similar or shorter as compared with the previous test presented in Table <ref>.







§ CONCLUSIONS


In this paper we propose an expanded version of the Local Variance Gamma model of
<cit.> which we refer as an Expanded Local Variance Gamma model, or ELVG. Two main improvements are introduced as compared with the LVG model. First, we add drift to the governing underlying process. It turns out that this a relatively minor (at the first glance) improvement requires a interesting trick to preserve tractability of the model, which is a non-trivial time-change. We show that still in this new model it is possible to derive an ordinary differential equation for the option price which plays a role of Dupire's equation for the standard local volatility model.

The second novelty of the paper as compared with the LVG model is that in contrast to <cit.> we consider a local variance to be a piecewise linear function of strike, while in <cit.> it was piecewise constant. We proceed in the spirit of <cit.> by describing a no-arbitrage interpolation, and then construct a closed-form solution of our ODE in terms of hypergeometric and generalized hypergeometric functions. An important advantage of this approach is that calibration of the model to market smiles does not require solving any optimization problem, and can be done term-by-term by solving a system of non-linear algebraic equations for each maturity, which, in general, is significantly faster, especially since we provide an algorithm for constructing a smart initial guess. We also provide various asymptotic solutions which allow a significant acceleration of the numerical solution and improvement of its accuracy in the corresponding cases (i.e, when parameters of the model at some iteration obey the conditions to apply the corresponding asymptotic).

In principle, somebody could claim that solving a system of nonlinear equations with a generic solver is not much different from solving a nonlinear optimization problem. Obviously, when our ODE is used as an alternative to the Dupire equation, the difference comes from the fact that calibration based on the Dupire equation requires solving this PDE at every iteration by either numerically, or semi-analytically by using a Laplace transform, which is obviously slower. As was mentioned in Introduction there exist many other calibration algorithms which reduce to a nonlinear optimization problem (e.g., taking a sufficiently large parametric family of local volatility functions and choosing the parameters that provide the best fit of observed prices). For the latter computation of the objective function is fast, but optimization must be constrained to preserve no-arbitrage, and, thus, slow.

In our numerical test we use same market data as in <cit.>. The results of the test demonstrate robustness of the proposed approach from both the speed and accuracy point of view, especially in cases where the above referred papers experienced some difficulties with achieving a perfect fit. An additional test performed for the S&P500 data taken from <cit.> gives rise to the same conclusion.







§ REFERENCES






 
   
27
natexlab#1#1









[#1],#1
[Abramowitz and Stegun(1964)]as64
authorAbramowitz, M., authorStegun, I.,
  year1964.
titleHandbook of Mathematical Functions.
publisherDover Publications, Inc.

[Ahoniemi(2009)]callPutIV
authorAhoniemi, K., year2009.
titleModeling and forecasting implied volatility.
Ph.D. thesis. Helsinki School of Economics.

[Balaraman(2016)]Balaraman2016
authorBalaraman, G., year2016.
titleModeling Volatility Smile and Heston Model
  Calibration Using QuantLib Python.
noteAvailable at
  <http://gouthamanbalaraman.com/blog/volatility-smile-heston-model-calibration-quantlib-python.html>.

[Bochner(1949)]BochnerPDE1949
authorBochner, S., year1949.
titleDiffusion equation and stochastic processes, in:
  booktitleProceedings of the National Academy of Sciences, USA,
  pp. pages368–370.

[Carr and Nadtochiy(2014)]CarrLVGOrig2014
authorCarr, P., authorNadtochiy, S.,
  year2014.
titleLocal variance gamma and explicit calibration to
  option prices.
noteAvailable at <https://arxiv.org/abs/1308.2326>.

[Carr and Nadtochiy(2017)]CarrNadtochiy2017
authorCarr, P., authorNadtochiy, S.,
  year2017.
titleLocal Variance Gamma and explicit calibration to
  option prices.
journalMathematical Finance volume27,
  pages151–193.

[Coleman et al.(2001)Coleman, Kim, Li and Verma]Coleman2001
authorColeman, T., authorKim, Y., authorLi,
  Y., authorVerma, A., year2001.
titleDynamic hedging with a deterministic local volatility
  function model.
journalThe Journal of Risk volume4,
  pages63–89.

[Cox and Rubinstein(1985)]CoxRubinstein1985
authorCox, J., authorRubinstein, M.,
  year1985.
titleOptions Markets.
publisherPrentice-Hall.

[Derman and Kani(1994)]derman/kani:94
authorDerman, E., authorKani, I., year1994.
titleRiding on a smile.
journalRISK , pages32–39.

[Dupire(1994)]Dupire:94
authorDupire, B., year1994.
titlePricing with a smile.
journalRisk volume7, pages18–20.

[Ekström and Tysk(2012)]Tysk2012
authorEkström, E., authorTysk, J.,
  year2012.
titleDupire's equation for bubles.
journalInternational Journal of Theoretical and Applied
  Finance volume15, pages1250041–1250053.

[Falck and Deryabin(2017)]FalckDeryabin2017
authorFalck, M., authorDeryabin, M.,
  year2017.
titleLocal variance gamma revisited.
noteAvailable at
  <https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2659728>.

[Hull(1997)]hull:97
authorHull, J.C., year1997.
titleOptions, Futures, and other Derivative Securities.
editionthird ed., publisherPrentice-Hall,
  Inc., addressUpper Saddle River, NJ.

[Itkin(2015)]ItkinSigmoid2015
authorItkin, A., year2015.
titleTo sigmoid-based functional description of the
  volatility smile.
journalNorth American Journal of Economics and Finance
  volume31, pages264–291.

[Itkin(2017)]ItkinBook
authorItkin, A., year2017.
titlePricing derivatives under Lévy models.
Number number12 in seriesPseudo-Differential
  Operators. edition1 ed., publisherBirkhauser,
  addressBasel.

[Itkin and Lipton(2018)]ItkinLipton2017
authorItkin, A., authorLipton, A.,
  year2018.
titleFilling the gaps smoothly.
journalJournal of Computational Sciences
  volume24, pages195–208.

[Lipton and Sepp(2011)]LiptonSepp2011
authorLipton, A., authorSepp, A., year2011.
titleCredit value adjustment in the extended structural
  default model, in: booktitleThe Oxford Handbook of Credit
  Derivatives. publisherOxford University, pp.
  pages406–463.

[Madan et al.(1998)Madan, Carr and Chang]MadCarrChang
authorMadan, D., authorCarr, P.,
  authorChang, E., year1998.
titleThe variance gamma process and option pricing.
journalEuropean Finance Review volume2,
  pages79–105.

[Madan and Seneta(1990)]MadanSeneta:90
authorMadan, D., authorSeneta, E.,
  year1990.
titleThe variance gamma (V.G.) model for share
  market returns.
journalJournal of Business volume63,
  pages511–524.

[Nayfeh(2000)]Nayfeh2000
authorNayfeh, A.H., year2000.
titlePerturbation methods.
publisherJohn Wiley & Sons.

[Ng and Geller(1970)]kummerInt1970
authorNg, E., authorGeller, M., year1970.
titleOn some indefinite integrals of confluent
  hypergeometric functions.
journalJournal of reserach of the Notional Bureau of
  Standards - B. Mathematica l Sciences volume748,
  pages85–98.

[Olver(1997)]Olver1997
authorOlver, F., year1997.
titleAsymptotics and Special Functions.
publisherAKP Classics.

[Polyanin and Zaitsev(2003)]PolyaninSaitsevODE2003
authorPolyanin, A., authorZaitsev, V.,
  year2003.
titleHandbook of exact solutions for ordinary differential
  equations.
edition2nd ed., publisherCRC Press Company,
  addressBoca Raton, London, New York, Washington, D.C.

[Revuz and Yor(1999)]RevuzYor1999
authorRevuz, D., authorYor, M., year1999.
titleContinuous Martingales and Brownian Motion.
edition3rd ed., publisherSpringer,
  addressBerlin, Germany.

[Shreve(1992)]Shreve:1992
authorShreve, S., year1992.
titleMartingales and the theory of capital-asset pricing.
journalLecture Notes in Control and Information SCIENCES
  volume180, pages809–823.

[Vasil'eva et al.(1995)Vasil'eva, Butuzov and Kalachov]VasBut1995
authorVasil'eva, A., authorButuzov, V.,
  authorKalachov, L., year1995.
titleThe boundary function method for singular
  perturbation problems.
Studies in Applied Mathematics, publisherSIAM,
  addressPhiladelphia.

[Wasow(1987)]Wasow1987
authorWasow, W., year1987.
titleAsymptotic Expansions for Ordinary Differential
  Equations.
publisherDover Pubns.






