An arbitrage-free interpolation of class C 2 for option prices
AArticle
An arbitrage-free interpolation of class C for optionprices Fabien Le Floc’h
Delft Institute of Applied Mathematics, TU Delft, Delft, The Netherlands * Correspondence: [email protected]: date; Accepted: date; Published: date
Abstract:
This paper presents simple formulae for the local variance gamma model of Carr and Nadtochiy,extended with a piecewise-linear local variance function. The new formulae allow to calibrate the modelefficiently to market option quotes. On a small set of quotes, exact calibration is achieved under onemillisecond. This effectively results in an arbitrage-free interpolation of class C . The paper proposesa good regularization when the quotes are noisy. Finally, it puts in evidence an issue of the modelat-the-money, which is also present in the related one-step finite difference technique of Andreasen andHuge, and gives two solutions for it. Keywords: volatility surface; interpolation; arbitrage; quantitative finance; European options
1. Introduction
The financial markets provide option prices for a discrete set of strike prices and maturity dates. Inorder to price over-the-counter vanilla options with different strikes, or to hedge complex derivativeswith vanilla options, it is useful to have a continuous arbitrage-free representation of the option prices, orequivalently of their implied volatilities. For example, the variance swap replication of Carr and Madanconsists in integrating a specific function over a continuum of vanilla put and call option prices [1,2]. Moregenerally, Breeden and Litzenberger [3] have shown that any path-independent claim can be valued byintegrating over the probability density implied by market option prices. An arbitrage-free representationis also particularly important for the Dupire local volatility model [4], where arbitrage will translate to anegative local variance. A option price representation of class C is also key to guarantee the second-orderconvergence of numerical schemes applied to the Dupire partial differential equation, commonly used toprice exotic financial derivative contracts. In this paper, we describe a new technique to interpolate themarket option prices in an arbitrage-free manner, with very high accuracy.A rudimentary, but popular representation is to interpolate market implied volatilities with a cubicspline across option strikes. Unfortunately this may not be arbitrage-free as it does not preserve theconvexity of option prices in general. A typical convex interpolation of the call option prices by quadraticsplines or rational splines is also not satisfactory in general since it may generate unrealistic oscillationsin the corresponding implied volatilities, as evidenced in [5]. Kahalé [6] designs an arbitrage-freeinterpolation of the call option prices, which however requires convex input quotes, employs twoembedded non-linear minimizations, and it is not proven that the algorithm for the interpolation functionof class C converges.More recently, Andreasen and Huge [7] have proposed to calibrate the discrete piecewise constantlocal volatility corresponding to a single-step finite difference discretization of the forward Dupire equation.In their representation of the local volatility, the authors use as many constants as the number of marketoption strikes for an optimal fit. It is thus sometimes considered to be "non-parametric". Their techniqueworks well in general but requires some care around the choice of discretization grid: it must be sufficiently a r X i v : . [ q -f i n . P R ] M a y of 20 dense so that two market strikes do not fall in between the same consecutive grid nodes, and sufficientlywide to properly model the boundary behaviour. Those two requirements complicate, and slow downthe non-linear optimization involved in the technique. Furthermore the output is a discrete set of optionprices, which, while relatively dense, must still be interpolated carefully to obtain the price of optionswhose strike falls in between grid nodes.Le Floc’h and Oosterlee [8] derived a specific B-spline collocation to fit the market option prices, whileensuring the arbitrage-free property at the same time. While the fit is quite good in general, it may not beapplicable to interpolate the original quotes with high accuracy. For example, input quotes may already besmoothed out if they stem from a prior model, or from a market data broker, or from another system in thebank. In those cases, it is desirable to use a nearly exact interpolation.Here, we extend the local variance gamma arbitrage-free interpolation of class C presented in [9].Instead of a piecewise-constant parametrization for the variance function, we use a piecewise-linearrepresentation. The resulting partial differential difference equations can still be solved explicitly and theoutcome is an interpolation of class C . Being based on the local variance gamma model, we find that itbehaves properly, even on challenging examples. Furthermore, we explain how the technique may beapplied to handle multiple option maturities, in a similar fashion as the method of Andreasen and Huge[7]. After writing this paper, we found out that the idea of using a piecewise-linear parametrization in thelocal variance gamma model was also developed in [10,11], with an additional stochastic drift term. Asour focus is on building an arbitrage-free interpolation of market option quotes, and not on using the localvariance gamma model to simulate the value of financial derivatives, we do not add a stochastic drift, andas a consequence, the equations we present here are much simpler. This allows for a significantly fastercalibration.We then propose a good regularization to calibrate the model against noisy option quotes. Finally,we put in evidence a flaw of the local variance gamma model, shared by the one-step finite differencetechnique of Andreasen and Huge, when the local variance function is at least of class C at the money,and propose two solutions for it.
2. Dupire’s PDDE in the local variance gamma model
We recall Dupire’s partial difference differential equation (PDDE) for a call option price C ( T , x ) ofstrike x and maturity T [9]: C ( T , x ) − max ( X ( ) − x , 0 ) T = a ( x ) ∂ C ( T , x ) ∂ x , (1)for a Martingale asset price process X ( t ) of expectation E Q [ X ( t )] = X ( ) .Let { x , x , ..., x m , x m + } be a increasing set of the strike prices, such that x = L , x m + = U with theinterval ( L , U ) being the spatial interval where the asset X lives. Furthermore, we require the following tohold ∃ s ∈ [ m ] | x s = X ( ) .In practice, ( x , ..., x m ) will correspond to the strike prices of the options of maturity T we want to calibrateagainst, along with the forward price.We consider a to be a continuous piecewise-linear function with values ( α i ) i = m at the knots ( x , ..., x m ) . The local variance gamma volatility function reads a ( x ) = α i + − α i x i + − x i ( x − x i ) + α i for x i ≤ x < x i + , i =
0, ..., m . (2) of 20 Let V be the function defined by V ( x ) = C ( x , T ) − max ( X ( ) − x , 0 ) . V is effectively the price of anout-of-the-money option (the price of a call option for x > X ( ) and of a put option for x < X ( ) ). TheDupire PDDE leads to V ( x ) = a ( x ) TV (cid:48)(cid:48) ( x ) , (3)on the intervals ( L , X ( )) and ( X ( ) , U ) . the function V is not a solution of the above equation on thewhole interval ( L , U ) since V (cid:48) ( x ) jumps at x = X ( ) . Indeed, the continuity of ∂ C ∂ x at x = X ( ) implieslim x → X ( ) − V (cid:48) ( x ) = + lim x → X ( ) + V (cid:48) ( x ) . (4)In order to define a unique V , we also impose the absorbing boundary conditions V ( L ) = = V ( U ) . (5)
3. Explicit solution
When α i + = α i , the local variance gamma functional a is constant on the interval [ x i , x i + ] . Thesolution to Equation 3 on this interval is given in Carr and Nadtochiy [9] and reads V ( x ) = χ i ( x ) [ Θ ci cosh ( ω i ( z i ( x ) − z i ( x i ))) + Θ si sinh ( ω i ( z i ( x ) − z i ( x i )))] , (6)with z i ( x ) = x , ω i = α i (cid:114) T , χ i = V on eachinterval is similar to a tension spline and we may use the fast and stable algorithm of Renka [12] for thefunctions coshm ( x ) = cosh ( x ) − ( x ) = sinh ( x ) − x , which gracefully handles the case whenthe argument | x | < α i + (cid:54) = α i , the solution is of the same form, but with z i ( x ) = ln | x + r i q i | , ω i = (cid:115) + q i T , χ i = (cid:118)(cid:117)(cid:117)(cid:116) x + r i q i x i + r i q i .with q i = α i + − α i x i + − x i , r i = α i − q i x i . The absolute value in z i ( x ) handles the case where the slope q i is negative. Proof.
On the interval [ x i , x i + ] , we have a ( x ) = q i x + r i = α i + − α i x i + − x i ( x − x i ) + α i . In particular, a ( x ) willtake values in the interval [ α i , α i + ] . As α i and α i + are assumed positive, a ( x ) is also positive.When the slope q i is positive, x + r i q i = a ( x ) q i is positive. We may apply the change of variable z = ln ( x + r i / q i ) , to Equation 3. This then leads to the following ODE V ( z ) = q i T ∂ V ∂ z − q i T ∂ V ∂ z .This is a second-order linear ODE with constant coefficients and the general solution is known to be of theform V ( z ) = A i e λ + i ( z ( x ) − z ( x i )) + B i e λ − i ( z ( x ) − z ( x i )) , of 20 where λ + i , λ − i are the roots of the quadratic λ − λ − Tq i = λ + i = + (cid:113) + q i T λ − i = − (cid:113) + q i T q i is negative, − x − r i / q i = a ( x ) − q i is positive and we may apply the change of variable z ( x ) = ln ( − x − r i / q i ) . This then leads to the same ODE as above.The two cases may be merged together by defining transform z ( x ) = ln (cid:12)(cid:12)(cid:12) x + r i q i (cid:12)(cid:12)(cid:12) , with z (cid:48) ( x ) = x + riqi .The derivative of V reads V (cid:48) ( x ) = χ i ( x ) z (cid:48) i ( x ) [( κ i Θ ci + ω i Θ si ) cosh ( ω i ( z i ( x ) − z i ( x i ))) + ( κ i Θ si + ω i Θ ci ) sinh ( ω i ( z i ( x ) − z i ( x i )))] ,(7)with κ i = (cid:40) α i = α i + , otherwise , (8) z (cid:48) i ( x ) = α i = α i + , x + riqi otherwise . (9)The boundary condition at x = x = L translates to θ c =
0. One may choose θ s arbitrarily to start thealgorithm. and its relation with V (cid:48) is θ s = V (cid:48) ( L ) ω . The final value will result from the continuity conditionsat i = s .The conditions to impose continuity of V and its derivative at x = x i + results in the following linearsystem cosh i Θ ci + sinh i Θ si = Θ ci + χ i ( x i + ) , (10) ( κ i cosh i + ω i sinh i ) Θ ci + ( ω i cosh i + κ i sinh i ) Θ si = (cid:0) κ i + Θ ci + + ω i + Θ si + (cid:1) z (cid:48) i + ( x i + ) χ i ( x i + ) z (cid:48) i ( x i + ) (11)for i =
0, ..., s −
2, withcosh i = cosh ( ω i ( z i ( x i + ) − z i ( x i ))) , sinh i = sinh ( ω i ( z i ( x i + ) − z i ( x i ))) .This system is solved from i =
0, starting with the calculation of θ ci + as given by Equation 10. The value θ ci + is then used to compute θ si + through Equation 11.At x = x m + = U , the boundary condition translates to θ sm = V (cid:48) ( U ) cosh m ω m z (cid:48) m ( x m + ) , θ cm = − θ sm sinh m cosh m .Similarly to the lower boundary, we may start with an arbitrary θ sm and then solve Equations 10, 11,downwards from i + i , for i = m − m −
2, ..., s . . When cosh i and sinh i are very large (for example,in the case of a small slope q i ), the system may be ill-defined numerically and extra care must be taken. In of 20 this case, a simple solution is to solve approximately the system, by adding a small regularization constant:when sinh i cosh i = ± ± ∓ (cid:101) .The free parameters θ s and θ sm are determined by the jump condition at i = s defined by V s − ( x s ) = V s ( x s ) and V (cid:48) s − ( x s ) = + V (cid:48) s ( x s ) . If θ s is multiplied by a factor ρ L , the system above implies that θ ci , θ si will also be multiplied by ρ L for i < s , and as a consequence, V s − will be multiplied by ρ L . Similarly, if θ sm is multiplied by a factor ρ R , V s will be multiplied by ρ R . We want to adjust θ s and θ sm by so that θ s ρ L and θ sm ρ R verify exactly the jump condition. This leads to the system V s − ( x s ) ρ L = V s ( x s ) ρ R , V (cid:48) s − ( x s ) ρ L = + V (cid:48) s ( x s ) ρ R ,with V s − ( x s ) = χ s − ( x s )( θ cs − cosh s − + θ ss − sinh s − ) , V s ( x s ) = θ cs , V (cid:48) s − ( x s ) = χ s − ( x s ) z (cid:48) s − ( x s ) (cid:2) ( κ s − θ cs − + ω s − θ ss − ) cosh s − +( ω s − θ cs − + κ s − θ ss − ) sinh s − (cid:3) , V (cid:48) s ( x s ) = ( κ s θ cs + ω s θ ss ) z (cid:48) s ( x s ) .And the solution is ρ R = V s ( x s ) V s − ( x s ) V (cid:48) s − ( x s ) − V (cid:48) s ( x s ) , (12) ρ L = V s ( x s ) V s − ( x s ) ρ R . (13)In order to make the solution verify the jump condition, we then multiply ( θ ci ) i = s − , ( θ si ) i = s − by ρ L ,and ( θ ci ) i = s ,..., m , ( θ si ) i = s ,..., m by ρ R .So far, we have ensured that V is a solution to Equation 3 of class C on the intervals ( L , X ( )) and ( X ( ) , U ) . By definition, C ( x , T ) = V ( x ) + max ( X ( ) − x , 0 ) . For i ∈ {
1, ..., s − s +
1, ..., m } , we have V (cid:48)(cid:48) i ( x i ) = a ( x i ) V i ( x i ) , and V (cid:48)(cid:48) i − ( x i ) = a ( x i ) V i − ( x i ) . By construction, we know that V i ( x i ) = V i − ( x i ) and a is continuous, thus V (cid:48)(cid:48) is continuous at x i , and we deduce that C ∈ C ( L , X ( )) , C ∈ C ( X ( ) , U ) .At x = X ( ) , the jump condition and the continuity of a ensures that ∂ C ∂ x is also continuous. We have thusobtained a solution of class C to the Dupire PDDE.
4. Calibration
The calibration of a single maturity consists in finding the parameters α , ..., α m + such that the marketoption prices ( ˆ C i ) i = m of respective strikes ( x i ) i = m match exactly the function C ( x , T ) , solution of theDupire PDDE of class C , at each market option strike.The problem is under-determined as we have m + m reference prices. In order toresolve this discrepancy, the parameters α and α m + may be set arbitrarily to control the extrapolation. Wechoose a flat extrapolation in terms of the local variance gamma function a , that is, α = α and α m + = α m .In general, the market strikes will not include X ( ) . In this case, X ( ) must be added to the knots { x i } i = m used in the local variance gamma representation. This adds one more parameter α s to the of 20 representation, where s is the index corresponding to X ( ) in the set of knots. We take this parameter to bethe linear interpolation of the enclosing parameters: α s = x s − x s − x s + − x s − σ s + + x s + − x s x s + − x s − σ s − .The calibration problem may be expressed as the least-squares minimization of the error measure E defined by E = m ∑ i = µ i ( σ ( α , x i ) − ˆ σ i ) , (14)with α i > i =
1, ..., m and where σ ( α , x ) is the implied volatility corresponding to the option pricesobtained with the piecewise-linear local gamma variance model and ˆ σ i is the market implied volatility atstrike x i , ( µ i ) i = m are weights associated to the accuracy of the fit at each point.In order to solve this non-linear least-squares problem, we will use the Levenberg-Marquardtalgorithm as implemented by Klare and Miller [13]. The box constraints α i > R to a subset of R + (for example throughthe function x → x + (cid:101) with some small positive (cid:101) ).The implied volatility for a given option price may be found efficiently and accurately through thealgorithm of Jäckel [15]. Alternatively, we may directly solve an almost equivalent formulation in terms ofoption prices, using the error measure E V defined by E V = m ∑ i = w i (cid:0) C ( α , x i ) − ˆ C i (cid:1) , (15)with C ( α , x ) being the local variance gamma option price with parameter α and strike x , and the cappedinverse Vega weights w i given by w i = min (cid:18) ν i , 10 X ( ) (cid:19) µ i , (16)where ν i = ∂ ˆ C i ∂σ is the Black-Scholes Vega corresponding the market option price ˆ C i , and 10 is a cap appliedto avoid numerical issues related to the limited machine accuracy (see Appendix A for an explanation ofthis relation between w and µ ). Another approach would be to use a fixed-point method, similar to the one described in [16,17].One could also explore solving exactly each α i successively using a one-dimensional non-linear solver,given an initial guess for A , and then either solve for A , or use a fixed point iteration on A . Eachiteration would require to solve again the ( α i ) i = m .We found the performance of a straightforward Levenberg-Marquardt minimization acceptable inpractice, and the least-squares approach to be more flexible if some additional regularization is needed. Asinitial guess, we found that various simple choices were almost equally effective such as: σ atm X ( ) , or σ atm x i , corresponding to an approximately flat Bachelier guess or an approximately flat lognormal guess,with σ atm being the (approximate) at-the-money Black-Scholes implied volatility in the set of options to fit. In the real world, the quoted options are likely to be on an asset with a non-zero time-dependent drift.For example, an equity will involve the interest rate and dividend yield evolution up to the maturity ofthe option, and a foreign exchange rate will involve the domestic and foreign interest rates evolution. Asdescribed in [18], it is always possible to translate the problem towards a problem on a driftless process X ,with call prices on a scaled strike. For example, for an asset S with forward price to time t , F ( t ) = E [ S ( t )] , of 20 we may consider X ( t ) = S ( t ) F ( t ) . The process X is a martingale and E [ X ( t )] =
1. Then, the price C of a calloption of maturity T and strike x on X is related to the price C of a call option of maturity T on S by C ( T , x ) = B ( T ) F ( T ) C ( T , F ( T ) x ) , (17)and thus the problem with strikes K i on the asset S is equivalent to a problem with strikes x i = K i F ( T ) on X . If a second maturity T is involved, with the same strikes K i , the original quotes will be converted toquotes on X at different set of strikes x i = K i F ( T ) .Let us consider two option maturities T , T with T > T >
0. Let § = (cid:8) x , ..., x m (cid:9) be the setof strikes for the first maturity and § = (cid:8) x , ..., x m (cid:9) the strikes for the second maturity. We start bycalibrating the first maturity as described in Section 4.1 using the knots { L , X ( ) , U } ∪ § . This leads tothe parameters ( α i ) . Then we add § − § ∩ § to the knots for the second maturity, and use a linearinterpolation to define the extra parameters α i corresponding to the set § − § ∩ § .We would like then to solve Dupire’s PDDE from T to T : C ( T , x ) − C ( T , x ) T − T = a ( x ) ∂ C ( T , x ) ∂ x , (18)If C ( T , x ) is obtained from a calibrated LLVG model between 0 and T and is thus given by Equation 6,then Equation 18 dictates that C ( T , x ) will not have the same analytical form, unless the local variancefunction a ( x ) is the same at both time-steps, or is a piecewise-constant function. We would however like tofind a piecewise-linear function a ( x ) at T , different from the one used to calibrate the options of maturity T , in order to fit the quotes at T , and then it is not obvious if there is an analytical solution at all.Instead of using a calibrated LLVG model for C ( T , x ) , we may consider a continuous piecewise linearinterpolation of the LLVG option prices at T . Typically we would use the market strikes as knots of thisinterpolation, but we may also use a finer grained representation, based on a calibrated LLVG model at T .Then, C ( T , x ) will be of the form given by Equation 6. The solution for θ c , θ s is however slightly different,since the first derivative must jump at each knot, instead of only at X ( ) , but the same technique is fullyapplicable. It is then guaranteed that at each knot of the piecewise linear interpolation, C ( T , x ) > C ( T , x ) ,which means that there is no calendar spread arbitrage at those points. It is not guaranteed everywhere,but we can always increase the number of knots if there is a concern.A much more straightforward approach, which practitioners sometimes apply, would be to calibratethe two maturities independently. As the calibration is almost exact, and option quotes are arbitrage-free,we know that this approach leads to no arbitrage at the discrete set of points corresponding to the union ofall the options we use in the calibration. Of course, there may however be some spurious calendar spreadarbitrage at a strike in between two quotes strikes. The hope is that it is a rare occurrence and we may findother ways to deal with those cases then.Clearly, this is an area where the technique of Andreasen and Huge [7] has an edge over the LLVGmodel, as the former will guarantee the absence of calendar spread arbitrages on a denser grid, and allowsfor any function a ( x ) (not only piecewise-constant or piecewise-linear). The input market prices, or equivalently, their implied volatilities, may not be arbitrage-free. In thiscase, the local variance α i will either move towards zero (in case of a calendar spread arbitrage) or towardsinfinity (in case of a butterfly spread arbitrage) and the minimization will not be exact. It will be thusimportant to cap and floor the local variance during the minimization. of 20 An alternative approach is to de-arbitrage the input quotes via the quadratic programmingrepresentation described in Appendix B. The quadratic programming problem is fast to solve usinga standard optimization library such as CVXOPT [19], OSQP [20] or quadprog [21], typically much fasterthan the calibration. A side-effect of this method is to remove any outlier automatically: they will effectivelybe smoothed out by the de-arbitraging process. In the direct approach described previously, outliers mayskew significantly the result, although one may remedy this by assigning smaller weights to the identifiedoutliers, which may be related to a large bid-ask spread. For example, an inverse bid-ask spread weightingscheme may be effective.Even if the quotes are de-arbitraged, they may still be noisy. As a consequence, the implied probabilitydensity would present many spikes, and the option gamma sensitivity may be of poor quality.In those cases, it may be useful to add regularization to the minimization as well. An interestingcandidate for the regularization is to minimize the strain energy of the beam that is forced to pass throughthe given data points [22]: E = m ∑ i = µ i ( σ ( α , x i ) − σ i ) + λ m − ∑ i = µ i σ (cid:48)(cid:48) ( α , x i )[ + σ (cid:48) ( α , x i ) ] ( x i + − x i ) ,The first term of the objective E corresponds to the square of the RMSE, while the second term is theregularization. The regularization parameter λ controls the smoothness of the spline interpolation. Wemay wish to use a simpler version, where we regularize directly the local variance function, by penalizingdifferences in the three-points estimate of its second derivative: E = m ∑ i = µ i ( σ ( α , x i ) − σ i ) + ˜ λ m − ∑ i = α i + − α i x i + − x i − α i − α i − x i − x i − , (19)with ˜ λ = X ( ) λ ∑ mi = µ i . We will see however that, while it leads to a smooth density, it also resultsin some spurious spike in the density, located at-the-money, at x = X ( ) . Our goal is to obtain a smoothprobability density, and we may directly perform the regularization on the density via E = m ∑ i = µ i ( σ ( α , x i ) − σ i ) + ˜ λ m − ∑ i = ln V (cid:48)(cid:48) ( x i + ) − ln V (cid:48)(cid:48) ( x i ) x i + − x i − ln V (cid:48)(cid:48) ( x i ) − ln V (cid:48)(cid:48) ( x i − ) x i − x i − , (20)where V (cid:48)(cid:48) ( x i ) = θ ci a ( x i ) = θ ci α i .
5. Numerical examples
Jäckel [5] shows that undesired oscillations can appear in the graph of the implied volatility againstthe option strikes when the option prices are interpolated by a monotonic and convex spline. Table A1 inappendix C presents a concrete example . Here, the option quotes are not direct market quotes, but thesolution of a sparse finite difference discretization of a local stochastic volatility model: the market neverquotes so far out-of-the-money option prices. His data has a few interesting properties: • some of the option prices are extremely small: the interpolation must be very accurate numerically. We are grateful to Peter Jäckel for kindly providing this data. of 20 • the option prices are free of arbitrage. In theory, an arbitrage-free interpolation can be exact. • a cubic spline interpolation on the volatilities or the variances, often used by practitioners, is notarbitrage-free. • a convexity preserving C -quadratic, or C -rational spline results in strong oscillations in the impliedvolatility.The interpolation proposed in [5] possesses unnatural spikes at the points of clamping, in particular, theimplied density is not continuous.We apply the quadratic convex spline interpolation of Schumaker [23], the rational convex spline ofSchaback [24], the one-step finite difference method of [7], and our linear local variance gamma (LLVG)interpolation on those quotes. As already evidenced in [5], the convex spline interpolations of the optionprices lead to unnatural oscillations in the implied volatility. The quadratic spline results in much strongeroscillations. The Andreasen-Huge technique is not exempt of oscillations, if a piecewise constant discretelocal volatility representation is used, as in Figure 1(b). The discontinuities in the representation leads to Log−moneyness I m p li ed v o l a t ili t y i n % Method
Andreasen−HugeLLVGReferenceSchabackSchumaker (a) Volatility smile for log-moneyness log xX ( ) larger than zero. Andreasen−Huge LLVG Schaback Schumaker−2 0 2 −2 0 2 −2 0 2 −2 0 20.000.250.500.751.00
Log−moneyness I m p li ed p r obab ili t y den s i t y (b) Probability density. Figure 1.
Implied volatility and probability density for Case I (Table A1 in appendix C). discontinuities in the implied probability density. If, instead of a piecewise constant function, we use a piecewise linear function, the oscillations in the implied volatilities with the technique of Andreasen andHuge disappear. The LLVG interpolation does not present any oscillations in the implied volatilities andresults in a smooth implied probability density.Furthermore, it is nearly exact on this example, and slightly faster than Andreasen and Huge techniqueapplied on a grid of 400 points (Table 1). It is however still more than an order of magnitude slower thanthe convex spline interpolations.
Table 1.
Root mean square error (RMSE) of the interpolated implied volatilities against the impliedvolatilities of Table A1 in Appendix C. For Andreasen-Huge and LLVG, the solver error tolerance is set to10 − . Method Case I Case IIRMSE Time (ms) RMSE Time (ms)LLVG 2 · − · − · − · − · − · − · − · − · − · − The second example (Case II of Table A1) is more challenging numerically, since some option pricesare very close to an arbitrage: the difference of consecutive option price slopes c i + − c i x i + − x i − c i − c i − x i − x i − , at strike x i = − , close to machine epsilon accuracy. As a consequence, the Levenberg-Marquardtsolver takes more iterations to find the solution for the LLVG model, and finds only an approximatesolution for Andreasen and Huge method (even on finer grids). All of the methods leads to a similarimplied volatility interpolation, and there is no oscillation on this example (Figure 2(a)).The optimal implied probability density corresponding to Case I is relatively smooth everywhere,and especially for log moneyness larger than zero. Figure 3 shows however that the probability densityimplied by technique of Andreasen and Huge [7] exhibits a staircase shape when zoomed-in, with bothpiecewise-constant and piecewise-linear representations. This is due to the interpolation in between thefinite difference grid nodes. In contrast, the probability density implied by the LLVG model stays verysmooth, and is continuous by construction. We consider option quotes on the SPX500 index of maturity 1 month, as of March 18, 2020, taken fromthe Chicago Board Options Exchange (CBOE). For each strike and maturity, at a given time, the marketquotes two prices for an option contract: the bid price and the ask price. In order to calibrate directlythe LLVG model to the market quotes, we need a single estimate of the implied cumulative probabilitydensity. It is common practice to use the average of the bid and ask prices, the mid price for this purpose.Alternatively, we could also build two distinct representations: one for the bid prices and one or the askprices.When taken separately, the bid, ask or mid prices are not guaranteed to be arbitrage-free in theory:there can be theoretical arbitrages within the bid-ask spread that can not be taken advantage of in practice.If we calibrate the LLVG model towards the closest arbitrage-free quotes, using the algorithmdescribed in Appendix B, the fit is not exact, mainly for two reasons: the problem is high-dimensionalas there are 344 quotes, and the problem is not well conditioned, since the optimal implied density isextremely jagged as evidenced in Figure 4(b).
Log−moneyness I m p li ed v o l a t ili t y i n % (a) Volatility smile for the LLVG model. Andreasen−Huge LLVG Schaback Schumaker0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 31e−191e−131e−071e−01
Log−moneyness I m p li ed p r obab ili t y den s i t y (b) Probability density in log-scale. Figure 2.
Implied volatility and probability density for Case II (Table A1 in appendix C).
The addition of a regularization during the calibration is thus a necessity. The regularization on theparameters α (Equation 19) is not appropriate here. It results in a spurious spike located at the forwardprice X ( ) = α is deeply linkedto the local variance gamma model itself. Indeed, if we set the parameters to a constant, for example α = X ( ) . This is relatively intuitive fromEquation 1 since the left hand side C ( T , x ) − max ( X ( ) − x , 0 ) is composed of a smooth function minus afunction with a discontinuous derivative at X ( ) .With noisy quotes, the regularization helps to significantly reduce the calibration time as the problemis better conditioned (Table 2). It takes around 1 second to calibrate the 344 option prices of maturity 1month, and 0.20 second to calibrate the 165 option prices of maturity 2 months, using an Intel Core i77600U processor. Strike I m p li ed p r obab ili t y den s i t y Method
Andreasen−HugeAndreasen−Huge (Linear)LLVG
Figure 3.
Implied probability density for the LLVG and Andreasen-Huge methods, for strike moneyness K ∈ [ ] , calibrated to the market data of Table A1 in Appendix C. Table 2.
Root mean square error (RMSE) of the interpolated implied volatilities against the impliedvolatilities of Table A1 in Appendix C.Method Case III Case IVRMSE Time (s) RMSE Time (s)LLVG regularized 0.00137 0.99 0.00027 0.20LLVG unregularized 0.00126 5.86 0.00002 4.00 ( ) The regularization on the density may not always fully solve the issue: when the number of knots(which correspond to the number of quotes) is small, and the peak density is not close to the forwardprice, the spike may appear in between the two knots which encompass the forward price. The Andreasenand Huge technique suffers from the same issue. When fitting to market data, it is not easy to find anexample, since, with a few points, the implied distribution will not tend to be very smooth in general. Asan illustration, we thus consider a manufactured example: we fit the LLVG model to 10 option prices ofstrikes (0.85, 0.90, 0.95, 1, 1.05, 1.1, 1.15, 1.2, 1.3, 1.4), obtained by the Black-Scholes model with constantvolatility σ B = T = α s is computed by linear interpolation, a large spike appearsin the probability density implied from a calibrated LLVG model, even with regularization (Figure 5).When we add a fictitious point at x = α s is left as a free parameter), and use aregularization constant λ = − during the calibration, the spike disappears and the implied probabilitydensity is much closer to the lognormal distribution. On this example, the RMSE in implied volatilities is3 · − with regularization, the fit is almost exact. An alternative which does not require regularization, isto choose α s such that our piecewise representation V is of class C in the interval ( x s − , x s + ) . By applyingthe derivative to V (cid:48)(cid:48) , along with Equation 3, the C continuity relation at x = x s reads (cid:18) Va (cid:19) (cid:48) ( x − s ) = (cid:18) Va (cid:19) (cid:48) ( x + s ) , Strike V o l R e f * Method
Mid QuotesRegularizedUnregularized (a) Volatility smile for the LLVG model.
Strike D en s i t y Method
Alpha RegularizationDensity RegularizationUnregularized (b) Probability density.
Figure 4.
Implied volatility and probability density of 1m SPX options as of March 18, 2020. where the notation x − s , x − s denotes the value of the limit towards x s respectively from the left and from theright. Using the continuity of V ( x s ) = θ cs , the jump condition of V (cid:48) at x s and the continuity of a ( x s ) = α s ,this leads to V (cid:48) ( x + s ) + α s − θ cs α s − α s − ( x s − x s − ) α s = V (cid:48) ( x + s ) α s − θ cs α s + − α s ( x s + − x s ) α s , Strike D en s i t y Method
Fictitious PointLinearLognormal
Figure 5.
Implied probability density for the LLVG model, using a fictitious point with regularization, or alinear interpolation for the parameter α s , fitted to a Black-Scholes model with constant volatility σ B = T = or equivalently α s = θ cs [ α s − ( x s + − x s ) + α s + ( x s − x s − )] θ cs ( x s + − x s − ) − ( x s + − x s )( x s − x s − ) . (21)This is not a linear problem, as θ cs depends on α s through θ cs + , θ ss + in a non-linear way (Equations 10 and11). Starting with the algorithm described in Section 2 to compute θ c , θ s , adopting a linear interpolation asinitial guess for α s , we may however apply the following iteration • Update α s through Equation 21. • Update θ cs , θ ss from θ cs + , θ ss + using this new α s via Equations 10 and 11. • Update θ cs − , θ ss − from θ cs − , θ ss − using this new α s via Equations 10 and 11. • Recalculate ρ R , ρ L to ensure the jump condition at x s , via Equations 12 and 13. Scale θ c , θ s by the new ρ R , ρ L .One iteration is good enough for practical purposes, three iterations is nearly exact (Figure 6). In order to illustrate the difference between the two calibration strategies outlined in Section 4.3, weconsider the market data of [6] for options on the SPX500 index as of October 1995. The quality of themarket date is not great, but it is a good illustration of what may happen in not-so-liquid markets.If we calibrate the LLVG model to each option maturity T i independently, from T = T = T i , andplot the total variance σ ( T i , x ) T i as a function of the log-moneyness y = ln xF ( T i ) , where σ is the calibratedLLVG model implied volatility and F ( T i ) = e ( r − q ) T i is the forward to maturity, with r , q respectively theinterest rate and dividend yield. It took around 0.4 ms to calibrate a single maturity and 2.2 ms to calibratethe full volatility surface. This is several orders of magnitude faster than the calibration time reported in[10, Table 4]. We notice that the lines for different maturities cross in the extrapolation part, around y = T i − to T i , using a linear interpolation of options prices at T i − , the lines do not Strike D en s i t y Method
Fictitious PointIteration 1Iteration 2Iteration 3Lognormal
Figure 6.
Zoom of Figure 5, showing probability density implied by the calibration of the LLVG modelwith one, two, or three iterations to update σ s . Log−moneyness V a r i an c e Maturity (a) Independent calibration.
Log−moneyness V a r i an c e Maturity (b) Bootstrap calibration.
Figure 7.
Implied variance against log-moneyness with two different calibration strategies for the LLVGmodel, on SPX500 options as of October 1995. cross anymore (Figure 7(b)). We however needed to be careful to add many option prices for the linearinterpolation at T i − , otherwise, the resulting smile was not necessarily as smooth.The implied volatility smiles are nearly the same within the interpolation range, even afterbootstrapping many maturities (see Figure 8).
6. Conclusion
We have presented simple formulae, which allow to efficiently calibrate the local variance gammamodel with a piecewise-linear representation of the local variance (LLVG model). In particular, we found We chose 50 equidistant option prices, 20 may be sufficient.6 of 20
Log−moneyness V o l a t ili t y i n % Calibration
BootstrapOneReference
Figure 8.
Implied volatility smile for the maturity T = the calibration to be orders of magnitude faster than what has been reported in [10]. An exact interpolationof a small set of option quotes typically requires less than one millisecond.For the calibration of a single maturity, the LLVG model possesses many advantages over theone-step finite difference technique of Andreasen and Huge [7]: it is faster, offers a continuous, smooth,interpolation, is more robust on the challenging examples of Jäckel [5], and there is no need to choose aproper discretization grid. We also have presented how to calibrate the LLVG model to multiple maturities,with a focus on avoiding calendar spread arbitrages.We proposed a specific regularization for the LLVG model, which leads to a smooth impliedprobability density, even when the options quotes used are noisy. We also put in evidence a flaw inthe LLVG model if the variance function a ( x ) is linear across the payoff discontinuity, or more generallyof class C . The one-step finite difference method suffers from the same issue. It may be solved by theproposed regularization, or via a simple iterative algorithm without regularization.We leave for further research an extension of the local variance gamma model with each piece of thetype a ( x ) = x q x + r . It is still analytically tractable as the corresponding Dupire PDDE can be seen as amodified Bessel equation.Another area of research would be to improve the calibration technique for an exact interpolation.Although we have found augmented Gauss-Newton solvers like Levenberg-Marquardt to work well, canthe calibration be reduced to a sequence of one-dimensional non-linear problems? Would it be more robustor more efficient? Funding:
This research received no external funding.
Conflicts of Interest:
The authors declare no conflict of interest.1. Carr, P.; Madan, D. Towards a theory of volatility trading.
Option Pricing, Interest Rates and Risk Management,Handbooks in Mathematical Finance , pp. 458–476.2. Carr, P.; Lee, R. Robust replication of volatility derivatives. Prmia award for best paper in derivatives, mfa2008 annual meeting, 2008.
3. Breeden, D.T.; Litzenberger, R.H. Prices of state-contingent claims implicit in option prices.
Journal of business , pp. 621–651.4. Dupire, B. Pricing with a smile.
Risk , , 18–20.5. Jäckel, P. Clamping Down on Arbitrage. Wilmott , , 54–69.6. Kahalé, N. An arbitrage-free interpolation of volatilities. Risk , , 102–106.7. Andreasen, J.; Huge, B. Volatility interpolation. Risk , , 76.8. Le Floc’h, F.; Oosterlee, C.W. Model-Free Stochastic Collocation for an Arbitrage-Free Implied Volatility, Part II. Risks , , 30.9. Carr, P.; Nadtochiy, S. Local variance gamma and explicit calibration to option prices. Mathematical Finance , , 151–193.10. Carr, P.; Itkin, A. An expanded local variance gamma model. arXiv preprint arXiv:1802.09611 .11. Carr, P.; Itkin, A. Geometric Local Variance Gamma Model. The Journal of Derivatives , , 7–30.12. Renka, R.J. Algorithm 716: TSPACK: Tension spline curve-fitting package. ACM Transactions on MathematicalSoftware (TOMS) , , 81–94.13. Klare, K.; Miller, G. GN–a Simple and Effective Nonlinear Least-Squares Algorithm for the Open SourceLiterature. .14. Kanzow, C.; Yamashita, N.; Fukushima, M. Levenberg-Marquardt methods for constrained nonlinear equationswith strong local convergence properties. J. Computational and Applied Mathematics , , 375–397.15. Jäckel, P. Let’s be rational, 2015.16. Reghai, A. The hybrid most likely path. Risk , , 34–35.17. Reghai, A.; Boya, G.; Vong, G. Local volatility: Smooth calibration and fast Usage. Available at SSRN 2008215 .18. Bühler, H. Volatility and dividends-volatility modelling with cash dividends and simple credit risk.
SSRNWorking Paper Series .19. Andersen, M.; Dahl, J.; Vandenberghe, L. CVXOPT: A Python package for convex optimization. abel. ee. ucla.edu/cvxopt .20. Stellato, B.; Banjac, G.; Goulart, P.; Bemporad, A.; Boyd, S. OSQP: An Operator Splitting Solver for QuadraticPrograms.
ArXiv e-prints , [arXiv:math.OC/1711.08013].21. Turlach, B.A.; Weingessel, A. quadprog: Functions to solve quadratic programming problems.
CRAN-Packagequadprog .22. Glass, J. Smooth-curve interpolation: A generalized spline-fit procedure.
BIT Numerical Mathematics , , 277–293.23. Schumaker, L.I. On shape preserving quadratic spline interpolation. SIAM Journal on Numerical Analysis , , 854–864.24. Schaback, R. Spezielle rationale splinefunktionen. Journal of Approximation Theory , , 281–292.25. Gatheral, J. The volatility surface: a practitioner’s guide ; Vol. 357, Wiley. com, 2006.
Appendix A Relation between the weights in the least squares minimization of option prices andimplied volatilities
We can find a weight w i that makes the solution similar to the one under the measures E and E V bymatching the gradients of each problem. We compare m ∑ i = w i ∂ C ∂α ( α , x i ) (cid:0) ˆ C i − C ( α , x i ) (cid:1) ,with m ∑ i = µ i ∂σ∂α ( α , x i ) ( ˆ σ i − σ ( α , x i )) . As we know that ∂ C ∂α = ∂σ∂α ∂ C ∂σ , we approximate ∂ C ∂σ by the market Black-Scholes Vega, the term (cid:0) ˆ C i − C ( α , x i ) (cid:1) by ∂ C ∂α ( α opt − α ) , and ( ˆ σ i − σ ( α , x i )) by ∂σ∂α ( α opt − α ) to obtain w i ≈ ∂ ˆ C i ∂ ˆ σ i µ i . (A1)In practice the inverse Vega needs to be capped to avoid taking into account too far out-of-the moneyprices, which won’t be all that reliable numerically and we take w i = min (cid:18) ν i , 10 F (cid:19) µ i , (A2)where ν i = ∂ c i ∂σ is the Black-Scholes Vega corresponding the market option price c i . Appendix B De-arbitraging the market option prices: quadratic programming problem
Let ( c i ) i = m be the undiscounted market call option prices on an asset S , with strike y i and forwardprice to maturity T given by f = E [ S ( T )] . The closest arbitrage-free option prices ˜ c are the solution thefollowing quadratic programming problem [8]:˜ c = argmin z ∈ R n + (cid:107) W · ( z − c ) (cid:107) (A3)subject to − < z i − z i − y i − y i − < z i + − z i y i + − y i < z i > max ( f − y i , 0 ) , for i =
2, ..., m − W is a diagonal matrix of weights. For equal weights, W is the identity matrix I m . We can includeinformation on the bid-ask spread, for example by taking w i to be the inverse of the bid-ask spread atstrike y i .We have (cid:107) W · ( z − c ) (cid:107) = z T W T Wz − ( W T Wc ) T z + ( Wc ) T Wc .The minimization problem can thus be formulated as a quadratic programming problem:˜ c = argmin z ∈ R m , Gz ≤ h z T Qz + q T z (A5)with Q = W T W , q = − W T Wc ,and the elements G i , j of the matrix G , that specifies the linear constraints in (A5), are G i , i − = − y i − y i − , G i , i = y i − y i − + y i + − y i , G i , i + = − y i + − y i , for i =
2, ..., m −
1, and G = y − y , G = − y − y , G m , m = y m − y m − , G m , m − = − y m − y m − .and the vector h by h = − (cid:101) , h i = − (cid:101) for 1 < i ≤ m . The lower bound constraint translates to G m + i + i = − h m + i + = − max ( f − y i , 0 ) − (cid:101) , for i =
1, ..., m .The constant (cid:101) defines a maximum acceptable slope and ensures that the call prices are strictly convex. Appendix C Market data
Table A1.
Black-Scholes implied volatilities against moneyness xX ( ) for an option of maturity T = Table A2.
Black-Scholes implied volatilities for SPX500 options in October 1995 from Nabil Kahale in [6].The spot price was S =
590 and the interest and dividend rates were r =