Model-free Stochastic Collocation For An Arbitrage-free .

Model-free stochastic collocation for an arbitrage-free 683Table 1 Collocation polynomials at the Gauss–Hermite nodes for the Black–Scholes model with volatility25%, and for the SVI parameters corresponding to a least-squares fit of SPX500 options of maturity 10 yearsxipiBlack–Scholes 20 yearsyiciSVI 10 yearsyici 3.3242574335521193 0.99955672840809971.302 56.569 1.8891758777537109 0.97056586707072936.475 60.2961.19690.181 0.6167065901925942 0.731285863197276726.861 24.38539.39212.0400.000608 88.6160.61670659019259420.26871413680272327106.662 11.119146.921 6.348255.167 991.210378.0420.457The polynomial is expressed with the coefficients ci as g N (x) 5ii 0 ci x3 Calibration of the stochastic collocation to market option pricesA Lagrange polynomial g N cannot always be used to interpolate directly on the collocation points implied by the market option strikes (yi )i 0,.,N , because on one sideN might be too large for the method to be practical (there are typically more thanhundred market option prices on the SPX500 equity index for a given maturity), andon the other side, there is no guarantee that the Lagrange polynomial will be monotonic, even for a small number of strikes. Grzelak and Oosterlee (2017) propose to relyon a set of collocation points (xi )i 0,.,N determined in an optimal manner from thezeros of an orthogonal polynomial. It corresponds to the set of the Hermite quadrature points in the case of the Gaussian distribution. This presupposes that we knowthe survival distribution function values at strikes which do not correspond to anyquoted market strike. In Grzelak and Oosterlee (2017), those values are given by theSABR model. Even with known survival distribution function values at the Hermitecollocation points, the resulting polynomial is not guaranteed to be monotonic. Forexample, we consider options expiring in 20 years on an asset with spot S 100 thatfollows the Black–Scholes model with a constant volatility σ 25%. The Lagrangecollocation polynomial of degree N 3 or N 5 implied by the Gauss–Hermitenodes is not monotonically increasing; we have g5 ( 2.34) 15.2 (see Table 1for the polynomial details). Another simple example we encountered comes fromfitting SPX500 options of maturity 10 years, with the Gatheral SVI parameterization(Gatheral and Lynch 2004). It corresponds to the following SVI parameters a 0.004,b 0.027, s 0.72, ρ 0.99 m 1.0. The corresponding Lagrange quintic polynomial obtained at the Gauss–Hermite nodes decreases around x 2.36 as we haveg5 ( 2.36) 16.5.In this paper, we don’t want to assume a prior model. Instead of using a Lagrangepolynomial g N to interpolate on well-chosen (xi )i 0,.,N as in the step 3 of the collocation method described in Sect. 2, we will directly calibrate a monotonic polynomialg N to the market option prices at strikes (yi )i 0,.,m , with m typically much larger thanN . The monotonicity will be guaranteed through a specific isotonic parameterization.The proposed parameterization will also conserve the first moment of the distributionexactly.123

684F. Le Floc’h, C. W. OosterleeIn order to apply the stochastic collocation directly to market option prices, we thusneed to:– find an estimate of the survival density from the market option prices (corresponding to the step 1 of the collocation method described in Sect. 2),– find a good initial guess for the monotonic polynomial g N ,– optimize the polynomial coefficients so that the collocation prices are closest tothe market option prices.We will detail each step.3.1 A rough estimate of the market survival densityKahalé (2004) proposes a straightforward estimate. Let (yi )i 0,.,m be the marketstrikes and (ci )i 0,.,m the market call option prices corresponding to each strike; thecall price derivative ci toward the strike K i can be estimated byci li li 1ci ci 1where li 2yi yi 1(6) l .for i 1, . . . , m 1, and with c0 l1 , cmmIf the market prices are arbitrage-free, that is when 1 ci 1 cici ci 1 0, for i 1, . . . , m 1yi yi 1yi 1 yi(7)it is guaranteed that 1 ci 0 and the ci are increasing. A more precise estimateconsists in using the parabola that passes through the three points ci 1 , ci , ci 1 toestimate the slopes:ci li (yi 1 yi ) li 1 (yi yi 1 )ci ci 1where li yi 1 yi 1yi yi 1(8) l . It will still lead to 1 c 0 andfor i 1, . . . , m 1, and with c0 l1 , cmmi increasing ci .We can build a continuous representation of the survival density by interpolatingthe call prices (yi , ci )i 0,.,m with the C 1 quadratic spline interpolation of Schumaker(1983), where additional knots are inserted to preserve monotonicity and convexity.2By construction, at each market strike, the derivative will be equal to each ci .The survival density corresponds toG Y (y) C(y), y(9)2 A C 1 polynomial spline on a fixed set of knots cannot preserve monotonicity and convexity in the generalcase (Passow and Roulier 1977).123

Model-free stochastic collocation for an arbitrage-free 685or equivalently through the put option prices P:G Y (y) 1 P(y). y(10)In practice, one will use out-of-the-money options to compute the survival densityusing alternative Eqs. (9) and (10). While in the case of the SABR model, it is importantto integrate the probability density (or the second derivative of the call price) from yto (Grzelak and Oosterlee 2017), here we are only interested in a rough guess.3From the survival density at the strikes yi , it is then trivial to compute the normalcoordinates xi .3.2 Filtering out the market call prices quotesIn reality, it is not guaranteed that the market prices are convex, because of the bid–askspread. While the collocation calibration method we propose in this paper will stillwork well on a nonconvex set of call prices, starting the optimization from a convexset has two advantages: a better initial estimate of the survival density and thus a betterinitial guess, and the use of a monotonic interpolation of the survival density.The general problem of extracting a “good”, representative convex subset of a nearlyconvex set is not simple. By “good”, we mean, for example, that the frontier definedby joining each point of the set with a line minimizes the least-squares error on thefull set, along with possibly some criteria to reduce the total variation. In Appendix C,we propose a quadratic programming approach to build a convex set that closelyapproximates the initial set of market prices. It can, however, be relatively slow whenthe number of market quotes is large. The algorithm takes 4.8 s on a Core i7 7700U,for the 174 SPX500 option prices as of January 23, 2018, from our example in Sect. 5.A much simpler approach is to merely filter out problematic quotes, i.e., quotesthat will lead to a call price derivative estimate lower than 1 or positive. We assumethat the strikes (yi )0 i m are sorted in ascending order. The algorithm starts from aspecific index k {0, . . . , m}. We will use k 0, but we let the algorithm to be moregeneric. The forward sweep to filter out problematic quotes consists then in:(i) Start from strike yk . Let the filtered set be S {(yk , x K )}. Let j k,c c (ii) Search for the next lowest index j, such that 1 y jj y j and j jj m. Replace j by j,(iii) Add (y j , c j ) to the filtered set S . Repeat steps (ii) and (iii).In our examples, we set the tolerance 10 7 to avoid machine epsilon accuracyissues close to 1. A small error in the derivative estimate near 1 or near 0 willlead to a disproportionally large difference in the coordinate x.While the above algorithm will not produce a convex set, we will see that it can besurprisingly effective to compute a good initial guess for the collocation polynomial.3 Integration is still possible with the quadratic spline interpolation approach.123

686F. Le Floc’h, C. W. OosterleeWe could also derive a similar backward sweep algorithm and combine the twoalgorithms to start at the strike yk close to the forward price F(0, T ). On our examples,this was not necessary.3.3 An initial guess for the collocation polynomialIn order to obtain an arbitrage-free price, it is not only important that the density(zeroth moment) sums up to 1, which the collocation method will obey by default, butit is also key to preserve the martingale property (the first moment), that is g N (x)φ(x) F(t, T ).(11)Using the recurrence relation [Eq. (4)], this translates toN 1a0 2 a2i (2i 1)!! F(t, T ).(12)i 1Instead of trying to find directly good collocation points, a simple idea for an initialNbk x k corresponding to the leastguess is to consider the polynomial h N (x) k 0squares fit of xi , yi :b0 , . . . , b N Nm min(a0 ,.,a N ) R N 1 i 0 ak xik yi 2 k 0 ,(13)with the additional martingality constraint. This is a linear problem and is very fastto solve, for example, by QR decomposition. Unfortunately, the resulting polynomialmight not be monotonic.As we want to impose the monotonicity constraint by a clever parameterization ofthe problem, we will only consider the least-squares (with additional martingality constraint) cubic polynomial as starting guess. The following lemma helps us determinewhether it is monotonic.Lemma 1 A cubic polynomial a0 a1 x a2 x 2 a3 x 3 is strictly monotonic andincreasing on R if and only if a22 3a1 a3 0.Proof The derivative has no roots if and only if the discriminant a22 3a1 a3 0If our first attempt for a cubic initial guess is not monotonic, we follow the ideaof Murray et al. (2016) and fit a cubic polynomial of the form A Bx C x 3 . For thisspecific case, the linear system to solve is then given by 1 001230m2i 0 x im4i 0 x i F(t, T )Am4 B m x y .i 0 i ii 0 x imm6Cxi3 yixi 0i 0 i0(14)

Model-free stochastic collocation for an arbitrage-free 687In our case, Lemma 1 reduces to B 0 and C 0. As the initial guess, we thususe the cubic polynomial with coefficients a0 A F(t, T ), a1 B , a2 0,a3 C .3.4 The measureThe goal is to minimize the error between specific model implied volatilities andthe market implied volatilities, taking into account the bid–ask spread. The impliedvolatility error measure corresponds then to the weighted root-mean-square error ofimplied volatilities: m22i 0 μi (σ (ξ, K i ) σi ) Mσ m2i 0 μi,(15)where σ (ξ, K i ) is the Black implied volatility4 obtained from the specific modelconsidered, with parameters ξ , σi is the market implied volatility and μi is the weightassociated with the implied volatility σi . In our numerical examples, we will chooseμi 1. In practice, it is typically set as the inverse of the bid–ask spread.An alternative is to use the root-mean-square error of prices: m22i 0 wi (C(ξ, K i ) ci ) MV m2i 0 wi,(16)where C(ξ, K i ) is the model5 option price and ci is the market option price at strikeK i . We can find a weight wi that makes the solution similar to the one under themeasure Mσ by matching the gradients of each problem. We comparem i 0withm 2wi2 C(ξ, K i ) (ci C(ξ, K i )) , ξ2μi 2i 0As we know that C ξ σ C ξ σ , σ(ξ, K i ) (σi σ (ξ, K i )) . ξwe approximatethe term (ci C(ξ, K i )) by C ξ (ξoptVega,to obtainwi C σby the market Black–Scholes ξ ), and (σi σ (ξ, K i )) by1 ci σiμi . σ ξ (ξopt ξ)(17)4 Fast and robust algorithms to obtain the implied volatility from an option price are given in Jäckel (2015)and Li and Lee (2011). When no implied volatility corresponds to the model option price, which can happenbecause of numerical error, we just fix the implied volatility to zero.5 In the case of the stochastic collocation, ξ corresponds to the coefficients of the collocation polynomial.123

688F. Le Floc’h, C. W. OosterleeIn practice, the inverse Vega needs to be capped to avoid taking into account toofar out-of-the-money prices, which won’t be all that reliable numerically and we take wi minwhere νi ci . ci σ1 106,νi F μi ,(18)is the Black–Scholes Vega corresponding to the market option price3.5 Optimization under monotonicity constraintsWe wish to minimize the error measure MV while taking into account the martingalityand the monotonicity constraints (Lemma 1) at the same time. The polynomial g Nis strictly monotonically increasing if its derivative polynomial is strictly positive.We follow the central idea of Murray et al. (2016) and express g N in an isotonicparameterization: xg N (x) a0 p(x)dx,(19)0where p(x) is a strictly positive polynomial of degree N 1 2Q. It can thus beexpressed as a sum of two squared polynomials of respective degrees at most Q andat most Q 1 (Reznick 2000):p(x) p1 (x)2 p2 (x)2 .(20)As in the case of the cubic polynomial, we can refine the initial guess by first findingthe optimal positive least-squares polynomial with the sum of squares parameterization. Let (β1,0 , . . . , β1,q ) Rq 1 be the coefficients of the polynomial p1 and(β2,0 , . . . , β2,q 1 ) Rq be the coefficients of the polynomial p2 . The coefficients(γk )k 0,.,N 1 of p can be computed by adding the convolution of β1 with itself tothe convolution of β2 with itself, that isγk k β1,l β1,k l l 0k β2,l β2,k l ,(21)l 0with β1,l 0 for l q and β2,l 0 for l q 1. The martingality condition leadstoN 1N2 γ2k 1γk 1 k(2k 1)!! x .g N (x) F(t, T ) (22)2kkk 1123k 1

Model-free stochastic collocation for an arbitrage-free 689Lemma 2 The gradient of g N toward (β1,0 , . . . , β1,q , β2,0 , . . . , β2,q 1 ) can be computed analytically, and we have βl,k x (2k 1)!! gNi(xi ) 2 βl,2k j 1 , βl, jk j 1kqqj k 1k 0(23)k 1with βl,k 0 for k 0 and β1,k 0 for k q and β2,k 0 for k q 1.Proof 2q kx g N (x) a0 β1,l β1,k l x k 0 k 0 l 0 a0 2q k 0We thus have ak 1 1k 11k 1 1k 1k β1,l β1,k l x2q 2k k 0 l 02q 2 k 1 l 0kl 0 β1,l β1,k lβ2,l β2,k l x k dxk 0kl 0 β2,l β2,k l 1 β2,l β2,k l x k 1 .k 1kl 0 kl 0 β1,l β1,k lfor 0 k 2q 2,for k 2q 1, 2q.We recall that the martingality condition impliesa0 F(t, T ) q a2k (2k 1)!!.k 1We have a0 a2k (2k 1)!!, βl, j βl, jqk 1and ak 112β1,k j for j k 2q, β1, jk 1 ak 112β2,k j for j k 2q 2, β2, jk 1 ak 1 0 for k j and l 1, 2. βl, jThus, a0 βl, jq(2k 1)!!βl,2k j 1 ,k 1k123

690F. Le Floc’h, C. W. Oosterleewith βl,k 0 for k 0 and β1,k 0 for k q and β2,k 0 for k q 1.In particular, for a cubic polynomial, we have g3 g3(xi ) 2β1,0 xi β1,1 xi2 β1,1 ,(xi ) 2β2,0 xi , β1,0 β2,0 g32β1,1 3x β1,0 ,(xi ) β1,0 xi2 β1,13 iand for a quintic polynomial, g52β1,2 3(xi ) 2β1,0 xi β1,1 xi2 x β1,1 , β1,03 i g5(xi ) 2β2,0 xi β2,1 xi2 β2,1 , β2,0 g52β1,1 3 β1,2 43xi xi β1,0 β1,2 ,(xi ) β1,0 xi2 β1,1322 g52β2,1 3(xi ) β2,0 xi2 x β2,0 , β2,13 i g52β1,0 3 β1,1 4 2β1,2 5 3x x x β1,1 .(xi ) β1,23 i2 i5 i2The cubic polynomial initial guess can be rewritten in the isotonic form as follows, a0 a1 x a2 x 2 a3 x 3 a0 x a1 2a2 t 3a3 t 2 dt0 a0 0x a23a3 t 3a3 2 2a22 dt. a1 3a3(24)Based on the initial guess (refined or cubic), we can use a standard unconstrainedLevenberg–Marquardt algorithm to minimize the measure MV , based on the isotonic parameterization. This results in the optimal coefficients (β1,0 , . . . , β1,q ) and(β2,0 , . . . , β2,q 1 ), which we then convert back to a standard polynomial representation, as described above.The gradient of the call prices toward the isotonic parameters can also be computedanalytically from Eq. (2), as we have C xK(K ) (K )C(K ) βl, j βl, j123 xK gN(x)φ(x)dx, βl, j

Model-free stochastic collocation for an arbitrage-free 691where x K is the integration cutoff point defined by Eq. (3). As gxN (β, x K ) 0, we xKcan use the implicit function theorem to compute the partial derivatives β(K ):l, j x K (β) g1N xwhere x K xK xK β1,0 , . . . , β2,q 1 (β, x K )and g N g N (β, x K ) gN gN β1,0 , . . . , β2,q 1 .4 Examples of equity index smilesWe consider a set of vanilla option prices on the same underlying asset, with the samematurity date. As an illustrating example, we will use SPX500 option quotes expiringon March 7, 2018, as of February 5, 2018, from Appendix D. The options’ maturityis thus nearly one month. The day before this specific valuation date, a big jump involatility across the whole stock market occurred. One consequence is a slightly moreextreme (but not exceptional) volatility smile.4.1 A short review of implied volatility interpolationsLet us recall shortly some of the different approaches to build an arbitrage-free impliedvolatility interpolation, or equivalently, to extract the risk-neutral probability density.We can choose to represent the asset dynamics by a stochastic volatility modelsuch as Heston (1993), Bates (1996), Double-Heston (Christoffersen et al. 2009).This implies a relatively high computational cost to obtain vanilla option prices andthus to calibrate the model, especially when time-dependent parameters are allowed.Furthermore, those models are known to not fit adequately the market of vanilla optionswith short maturities. Their implied volatility smile is typically too flat.Many practitioners revert to a parameterization-based or inspired from a stochasticvolatility model, such as the Hagan SABR expansion (Hagan et al. 2002), or theGatheral SVI model (Gatheral and Lynch 2004; Flint and Maré 2017). These aremuch faster to calibrate. SVI is one of the most popular parameterizations to representthe equity option volatility smile, because of its simplicity, its relation to stochasticvolatility models asymptotically, and its almost arbitrage-free property. However, aswe shall see, the fit for options on equities can still be poor (Fig. 1a). SVI managesto fit only a part of the left wing and fails to represent well the market quotes in theregion of high implied volatility curvature. SVI and SABR are usually much better atfitting longer option maturities.Another approach is not to assume any underlying model and use an exact interpolation. A cubic spline interpolation of the implied volatilities is not arbitrage-free,although it is advocated6 by M. Malz (2014). Kahalé (2004) proposes an arbitragefree spline interpolation of the option prices. Unfortunately, it is not guaranteed that6 Malz precises that the challenge of his approach is to find a good filter for the quotes, which he does notdescribe at all.123

692F. Le Floc’h, C. W. Oosterlee70Implied volatility in %60MethodAndreasen Huge50ReferenceSVI403020210024002700Strike(a) Implied volatility.0.020Probability densityMethodAndreasen Huge0.015SVI0.0100.0050.000210024002700Strike(b) Probability density.Fig. 1 SVI and Andreasen–Huge calibrations on 1m SPX500 options as of February 05, 2018his algorithm for C 2 interpolation, necessary for a continuous probability density,converges. Furthermore, it assumes that the input call option quotes are convex anddecreasing by strike. But the market quotes are not convex in general, mainly becauseof the bid–ask spread. While we propose a quadratic programming-based algorithm tobuild a convex set that closely approximates the market prices in Appendix C, it canbe relatively slow when the number of quotes is large. Finally, the resulting impliedprobability density will be noisy, as evidenced by Syrdal (2002).A smoothing spline or a least-squares cubic spline will allow to avoid overfittingthe market quotes. For example, Syrdal

maturity swaps (CMS) consistently with the swaption implied volatility smile. In the standard approximation of Hagan (2003), the CMS convexity adjustment consists in evaluating the second moment of the distribution of the forward swap rate. It can be computed in closed form with the stochastic collocation. This allows for an efficient