A First Option Calibration of the GARCH Diffusion Model by a PDE Method
1 A FIRST OPTION CALIBRATION OF THE GARCH DIFFUSION MODEL BY A PDE METHOD
Yiannis A. Papadopoulos and Alan L. Lewis ABSTRACT
Time-series calibrations often suggest that the GARCH diffusion model could also be a suitable candidate for option (risk-neutral) calibration. But unlike the popular Heston model, it lacks a fast, semi-analytic solution for the pricing of vanilla options, perhaps the main reason why it is not used in this way. In this paper we show how an efficient finite difference-based PDE solver can effectively replace analytical solutions, enabling accurate option calibrations in less than a minute. The proposed pricing engine is shown to be robust under a wide range of model parameters and combines smoothly with black-box optimizers. We use this approach to produce a first PDE calibration of the GARCH diffusion model to SPX options and present some benchmark results for future reference.
1. Introduction
Stochastic volatility models are a natural generalization of the seminal Black-Scholes-Merton (BSM) option theory. In such models, the constant volatility parameter π of the BSM theory is promoted to a random process: ππ π‘ = ππ π‘ ππ‘ + π π‘ π π‘ ππ π‘π . Indeed, there is general agreement in finance that volatility (in its many forms) is best modelled as some sort of mean-reverting stochastic process. Starting from that premise, there are many possibilities. One of the simplest has the instantaneous variance rate π£ π‘ β‘ π π‘2 evolving as a positive diffusion process following the SDE: ππ£ π‘ = π (π£Μ β π£ π‘ )ππ‘ + ππ£ π‘ ππ π‘π£ . Here π π‘π£ is an additional Brownian motion, (π , π£Μ , π) > 0 are constant parameters, and the two Brownian motions (π π‘π , π π‘π£ ) are correlated with constant parameter π . Coupled with the (risk-neutral) stock price evolution above , this defines the GARCH diffusion model. The GARCH diffusion model has several nice properties. First, ignoring the drift term for a moment, π£ π‘ evolves as a geometric Brownian motion (GBM) -- a natural way in finance to achieve a positive stochastic process. GBM was originally introduced into finance by M. F. M. Osborne in the 1950βs to model stock prices under constant volatility. Indeed, time series analysis seems to favor GBM volatility over the popular Heston β93 (square-root) volatility process. Second, with π£Μ = 0 , the model nests a variant of the SABR model β very popular in interest rate modelling. The virtue of the SABR-GARCH connection is very tractable small-time behavior β due to a close connection of the small-time dynamics with hyperbolic Brownian motion. (While tractable small-time behavior facilitates time-series analysis by Maximum Likelihood, we found it not especially helpful in option chain calibration). Finally, the model name comes from the property (due to D. Nelson) that there exists a continuous-time limit of a discrete-time GARCH model (GJR-GARCH) that leads to a GARCH diffusion model. How well can the model fit option chains? Answering that is called calibration. Unfortunately, one desirable β but absent β property is an analytic solution, leaving numerics. While a simulation-based (Monte Carlo) approach doesnβt seem like the most efficient approach , in fact the model has been calibrated using Monte Carlo to a large options data set extending over several years by Christoffersen et al. in [1]. They find the GARCH diffusion a better fit than the oft- calibrated Heston β93 model and the so-called 3/2-model, their points of comparison. Thessaloniki, Greece; email: [email protected] Newport Beach, California, USA; email: [email protected] Briefly summarized (with π = 0) in Bollerslev and Rossiβs 1995 D. Nelson remembrance piece [20]. Given the nice properties, prior calibration results, and the general challenge, we were motivated to develop an efficient, accurate PDE calibrator for this model. Here, we report our methods and first results. and KBEβs . An important property that the model shares with a wide class of models is stock price level-independence, a well-known scaling relation for vanilla option prices. Specifically, at some initial time π‘ , consider a vanilla European call option price πΆ(π‘ , π , π£ ; πΎ, π) with strike price πΎ , expiration π , and state variables (π , π£ ) . Then πΆ(π‘ , π , π£ ; πΎ, π) = πΎ π(π‘ , π , π£ ; π) where the standardized option pricing function π(π‘, π , π£; π) is independent of πΎ and π = π /πΎ . Fixing and suppressing (πΎ, π) , consider the pricing function πΆ(π‘, π, π£).
It satisfies the KBE (Kolmogorov backwards equation) problem: βππΆ/ππ‘ = πΏ
π,π£ πΆ with terminal condition πΆ(π, π, π£) = (π β πΎ) + , and where πΏ π,π£ is the process generator. Then, of course, π(π‘, π , π£ ) satisfies the same PDE with c (π, π , π£) = (π β 1) + . Now fix πΎ , say πΎ = πΎ β‘ π , and solve the (continuum) KBE problem once for expiration π . This gives π(π‘ , π , π£ ) , a function of π for π β (0, β) since π(π‘ , π , π£ ) = πΆ(π‘ , π πΎ , π£ ) / πΎ and the r.h.s is known for all values of π . For any other strike then, say πΎ = πΎ , one immediately gets πΆ(π‘ , π , π£ ; πΎ , π) = πΎ π(π‘ , π /πΎ , π£ ) . The point is that a single KBE solution yields all the (vanilla) option values for different strikes at a given expiration. While obvious in hindsight, the KBE implication of the MAP property initially eluded us. Early on we thought a forward equation (Fokker-Planck) was the only way for pricing βall -options-at- onceβ at a fixed expiration . Exploiting the scaling property resulted in significant performance improvements over our original βone option at -a- timeβ approach: 3Γ β π πππ‘ππππ /π ππ₯πππππ‘ππππ . The somewhat subtle reasons are discussed in Sec. 3.1. Given our KBE approach, one must make choices on how to solve the pricing PDE. As with the Heston model, option prices under the GARCH diffusion model are governed by a 2-D convection-diffusion-reaction PDE with a mixed derivative term. Key characteristics of a suitable numerical scheme would be a) stability under practical usage, b) good accuracy to execution time ratio, and c) robustness (good oscillation-damping properties). As noted in [2], spurious oscillations in numerical computation of option prices can have three distinct causes: convection dominance, time-stepping schemes that are unable to sufficiently damp the high frequency errors stemming from the payoff discontinuity, and finally negative coefficients arising from the discretization of the diffusion terms. Here we take a closer look at the last two. For the spatial discretization we use the finite difference method on non-uniform grids. We employ standard central finite difference formulae for the diffusion and convection terms, but opt for a less common formula for the mixed derivative term, one that helps reduce oscillations that may take the solution to negative values. Although not our first choice, we also discretize the PDE cast in the natural We are well-aware of the general limitations of simple stochastic volatility diffusions. For example, they have difficulty fitting short-dated SPX option smiles and VIX options. Overcoming the limitations seems to require jump-processes. But, even if you want to include jumps into so- called βnon - affineβ models (like the GARCH diffusion), you need to start with a good PDE solver. For American options, barrier options, and other more exotic options, individual KBE solutions are needed. More generally, replace π£ in the scaling argument above by π , a (D-1) vector-valued state variable for a D-dimensional jump- diffusion or whatever. Then, scaling (and thus βall -options-at- onceβ) for Euro -style vanillas holds if: the process ( π π‘ , π π‘ ) is a MAP (Markov Additive Process), where π π‘ β‘ log π π‘ is the additive component and π π‘ is the Markov component. MAPs are defined in Γinlar [19]; modelling implications are stressed in Lewis (2016) [5]. Note this generality admits even discrete-time processes. Thus, for MAPs, it suffices to solve the backwards evolution problem once for a single strike to get all the vanilla option prices at a given expiration π . Admitting jumps, the backwards evolution problem (in continuous-time) is generally a PIDE (partial integro-differential equation) problem. logarithm of the asset price; combined with the mixed derivative scheme, this can further guard against negative values (but not preclude them altogether). With spatial discretization in place, one is left with a large system of stiff ordinary differential equations and must adopt a time-marching method. We employ two commonly used schemes plus a rather unusual one. For this type of PDE the most popular choice would be a cross-derivative-enabled ADI variant (see [3] for an overview). We opt for the Hundsdorfer-Verwer and the Modified Craig-Sneyd schemes that offer the best overall characteristics. Our alternative is the BDF3 fully implicit scheme which, as far as we know, has not been used in such a context in the financial literature. It may have already become apparent that we do not aim for one sole scheme that is necessarily monotone by design; we believe that such a scheme would likely be less accurate or slower than it needs to be. What we aim for instead is a reliable set-up, enabling as fast and accurate calibrations as possible. To this end we propose a strategy that involves occasional re-evaluations and a hybrid engine that switches from ADI to the slower but more robust BDF3 scheme in such cases. The optimization is done with commercial software. We mainly use local constrained optimization routines, but we also try a global method (Differential Evolution). The latter, while proving too slow to be the recommended option, can be used to add confidence that the local optimizer is indeed finding global minima (which we have found to be the case in all our tests). The rest of this paper is organized as follows: Sec. 2 presents the numerical methods for the solution of the pricing PDE, with non-standard implementation specifics given in more detail. Sec. 3 describes the calibration phase and proposed strategy for optimizing performance. Sec. 4 contains various numerical results. We compare the computational efficiency of the time-marching schemes and examine the effectiveness of Richardson extrapolation in both space and time. This is followed by reference calibration results to real data and comparisons with the Heston model. We conclude with a brief exploration of other non-affine models that are readily handled by our framework. We finally present our conclusions and suggestions for further development.
2. Numerical solution of the GARCH diffusion PDE
The GARCH diffusion stochastic volatility model is described (under the risk-neutral measure) by ππ π‘ = (π π β π π )π π‘ ππ‘ + βπ£ π‘ π π‘ ππ π‘π , (1) ππ£ π‘ = π (π£Μ β π£ π‘ )ππ‘ + ππ£ π‘ ππ π‘π£ . Here the Brownian noises associated to the underlying asset π π‘ (here the SPX) and its variance π£ π‘ are correlated; i.e., ππ π‘π ππ π‘π£ = π ππ‘. A compatible real-world evolution is given in Appendix A. Time-series analysis (of similar real-world models) suggests that the correlation coefficient π is negative with typical values of around -0.75 (Ait-Sahalia & Kimmel [4]). So here we assume π < 0 . The variance process π£ π‘ has volatility π > 0 and reverts to its long-run mean π£Μ > 0 with a mean-reversion rate of π > 0. π is the time of an option expiration. Generally, our model assumes an environment with deterministic interest rate and dividend yields: ( π π‘, π π‘ ). But we write ( π π, π π ) to indicate that we are using stepwise constants for each option expiration. (There will be some deterministic behavior for ( π π‘, π π‘ ) compatible with this). Let then π(π, π£, π‘) denote the price of a European option when at time
π β π‘ the underlying asset price equals π and its variance equals π£ . It is easy to verify that under the above specification π(π, π£, π‘) must satisfy the following parabolic PDE ππππ‘ = 12 π π£ π πππ + ππππ£ π πππππ£ + 12 π π£ π πππ£ + (π π β π π )π ππππ + π (π£Μ β π£) ππππ£ β π π π (2) for
We can also cast the equation in terms of the natural logarithm of the price
π = ln(π) ππππ‘ = 12 π£ π πππ + πππ£ π πππππ£ + 12 π π£ π πππ£ + (π π β π π β 12 π£) ππππ + π (π£Μ β π£) ππππ£ β π π π (3) for
Equations (2) and (3) are categorized as time-dependent convection-diffusion-reaction PDE β s on an unbounded spatial domain. While (3) has a slightly simpler form, it is harder to allocate points on the π -grid optimally. This is especially true since in many cases the grid needs to start from very small S values (to avoid loss of accuracy from grid truncation), which then means that a lot of X -points will be placed in an area of low interest. We will thus discretize and solve (2) primarily, but that the code is also (trivially) adapted for switching to solving (3) as well. Initial and boundary conditions
As initial conditions to (2) we have the vanilla call and put payoffs π ππππ (π, π£, 0) = max(π β πΎ, 0), π ππ’π‘ (π, π£, 0) = max(πΎ β π, 0) , (4) where K is the strike of the option. We impose (numerical) boundary conditions of Dirichlet type (5) - (6) and Neumann type (7) - (9) on the left and right-side boundaries respectively: π ππππ (π πππ , π£, π‘) = 0, (5) π ππ’π‘ (π πππ , π£, π‘) = πΎπ βπ π π‘ β π πππ π βπ π π‘ , (6) ππ ππππ ππ (π πππ₯ , π£, π‘) = π βπ π π‘ , (7) ππ ππ’π‘ ππ (π πππ₯ , π£, π‘) = 0, (8) ππ ππππ ππ£ (π, π£ πππ₯ , π‘) = ππ ππ’π‘ ππ£ (π, π£ πππ₯ , π‘) = 0. (9) Under this model π£ = 0 is an entrance boundary for all (π , π) > 0 , meaning that π£ = 0 is unreachable whenever the process starts at π£ π > 0 . However, the process may in principle be started at π£ π = 0 , after which it immediately enters the interior and never hits the origin again (for more details the reader is referred to Lewis [5], pg. 102). Therefore, from a mathematical standpoint no boundary condition is necessary. The PDE itself can be applied at π£ = π£ πππ = 0 (where all diffusion terms vanish due to the presence of factor π£ ) and there is no need for any extra condition from a numerical point of view either. The choice of the grid truncation boundaries π πππ , π πππ₯ (or π πππ , π πππ₯ ) and π£ πππ₯ is discussed in Sec. 2.4.1.1. The boundary conditions are set in an equivalent manner for equation (3). We discretize in space using the finite difference method and work on non-uniform grids which we consider necessary for the efficient solution of the pricing PDE. In the S -direction, allocating more points around the strike can significantly reduce the error stemming from the initial delta discontinuity there. In the v -direction, allocating more points near π£ makes sense since we want to resolve better the area where we want to obtain a price. Also, since typically we have π£ πππ₯ β« π£ π (see Sec. 2.4.1.1), a non-uniform v -grid is all but necessary to both adequately resolve the area around π£ π and at the same time reach out to π£ πππ₯ with a reasonable number of grid points. We use the standard central finite difference formulas for the first and second derivatives in (2) and (3) and a rather less standard seven-point stencil representation for the mixed derivative. All formulas give second-order accurate approximations, provided the grid step variation is sufficiently smooth (as is indeed the case for the grid construction proposed in Sec. 2.4.1.2). This means that there are indeed non-trivial option price solutions for π£ = 0. Let the grid in the S -direction be defined by ππ + 1 points, 0 β€ π πππ = π < π < β― < π ππ =π πππ₯ and the corresponding grid steps Ξπ π = π π β π πβ1 , π = 1,2, β¦ , ππ . We then define the discretized versions of the first and second derivatives ππ π,π /ππ and π π π,π /ππ at π = π π as ππ π,π ππ β ββπ π+1 βπ π (βπ π + βπ π+1 ) π πβ1,π + βπ π+1 β βπ π βπ π βπ π+1 π π,π + βπ π βπ π+1 (βπ π + βπ π+1 ) π π+1,π , (10) π π π,π ππ β 2βπ π (βπ π + βπ π+1 ) π πβ1,π β 2βπ π βπ π+1 π π,π + 2βπ π+1 (βπ π + βπ π+1 ) π π+1,π . (11) We use the equivalent expressions for the derivatives ππ π,π /ππ£ and π π π,π /ππ£ at π£ = π£ π in the v -direction where the grid is defined by ππ + 1 points, < π£ < β― < π£ ππ£ = π£ πππ₯ and Ξv π = π£ π β π£ πβ1 , π = 1,2, β¦ , ππ . An exception is the π£ = 0 boundary where we use the one-sided (upwind) second-order formula for ππ π,π=0 / ππ£ : ππ π,0 ππ£ β β (2βπ£ + βπ£ )βπ£ (βπ£ + βπ£ ) π π,0 + (βπ£ + βπ£ )βπ£ βπ£ π π,1 β βπ£ βπ£ (βπ£ + βπ£ ) π π,2 . (12) For the mixed derivative term, we opt for a custom second-order scheme based on a 7-point stencil, which is very similar but not identical to that proposed by Ikonen & Toivanen in [6]. Such a scheme can be constructed so that it contributes fewer negative off-diagonal coefficients to the resulting systemβs discretization matrix A than the standard second-order scheme based on the 9-point stencil. This in turn makes the solution less likely to produce a negative valuation. When the correlation coefficient π is negative (which we take to be the case here as discussed in Sec. 2.1), an appropriate formula for approximating π π π,π /ππππ£ at (π, π£) = (π π , π£ π ) is given by π π π,π ππππ£ = 1π· (βπ π+1,πβ1 + 2π π,π β π πβ1,π+1 + (βπ π+1 β βπ π ) ππ π,π ππ + (βπ£ π+1 β βπ£ π ) ππ π,π ππ£ +
12 (βπ π+12 +βπ π2 ) π π π,π ππ + 12 (βπ£ π+12 +βπ£ π2 ) π π π,π ππ£ ) (13 ) where π· = βπ π+1 βπ£ π + βπ π βπ£ π+1 . (14) Formula (13) is readily obtained considering Taylor expansions of the option value π π,π at the neighboring upper left and lower right grid points (π πβ1 , π£ π+1 ) and (π π+1 , π£ πβ1 ) . Such a formula can be used in conjunction with a specially constructed grid (with limitations imposed on the grid steps) and some use of first-order upwind formulas for the convection terms ππ π,π ππβ and ππ π,π ππ£β , to make A an M-matrix by design (see [6] for example). This would ensure that the solution cannot produce a negative valuation in any case, which would be a particularly useful feature for our calibration (that requires the calculation of implied volatilities). Such an approach though is not favored in the present work, since we believe that it would unnecessarily reduce the average accuracy of the solution through suboptimal grid construction: The grid points allocation should be driven by the problemβs physical characteristics (e.g. the location of the payoff discontinuity) and not be forced upon through the mathematical requirement of nonnegative coefficients . We revisit this in Sec. 3 where we explain how we handle the occasional negative values that are indeed possible under the proposed discretization. We can now replace the spatial derivatives on the right-hand side of equation (2) with their discretized versions described above to obtain its semi-discretized form Using the first-order (two-point) upwind formula would be better from a stability point of view, but would result in loss of accuracy and the overall second-order convergence of the discretization. In practice we have seen no stability issues arising from the use of (12) in extensive tests throughout numerous calibration exercises. Zvan et al. [2] provide a similar discussion, albeit in the context of finite volume/element discretization . ππ π,π ππ‘ = π π π π π,π ππ + π π£ π π π,π ππ£ + π π ππ π,π ππ + π π£ ππ π,π ππ£ + π ππ£ (βπ π+1,πβ1 + 2π π,π β π πβ1,π+1 ) β π π π π,π , (15) where the diffusion, convection and mixed derivative coefficients are given by π π = 12 π π2 π£ π + 12 π ππ£ (βπ π+12 +βπ π2 ) , (16) π π£ = 12 π π£ π2 + 12 π ππ£ (βπ£ π+12 +βπ£ π2 ) , (17) π π = (π π β π π )π π + π ππ£ (βπ π+1 β βπ π ) , (18) π π£ = π (π£Μ β π£ π ) + π ππ£ (βπ£ π+1 β βπ£ π ) , (19) and π ππ£ = 1π· πππ π π£ π3 2β . (20) Equation (15) is applied at each grid point (π π , π£ π ) for π = 1, 2, β¦ ππ and π = 0,1, β¦ , ππ . We do not need to solve for π π=0 = π πππ since the Dirichlet boundary conditions (5) and (6) specify constant values there for the option value π . At π£ π=0 = 0 the second and mixed derivative terms in (15) vanish and the upwind discretization (12) means we only use values within our grid. This is not the case though for the far-boundary grid lines π π=ππ = π πππ₯ and π£ π=ππ = π£ πππ₯ : the Neumann-type conditions (7) - (9) imply that the mixed derivative π π ππππ£β (and thus all the terms in (15) - (19) multiplied by π ππ£ ) vanish. But we still have the second derivatives, whose stencils reference a point outside the grid. Such points are treated as fictitious and their value is obtained through extrapolation based on the last actual grid point and the known value of the gradient there. With spatial discretization in place we are now left with a large system of stiff ordinary differential equations (ODEβs) in time, which we can write as π½ β² (π‘) = π(π‘, π½(π‘)), π½(0) = π½ , where π(π, π½(π)) βΆ= π¨π½(π‘) + π(π‘) for 0 β€ π‘ β€ π. (21) Here π½ β² (π‘), π½(π‘), π(π‘) and π½ are vectors of size π and π¨ is the π Γ π spatial discretization matrix, where
π = ππ Γ (ππ + 1) is the total number of unknowns. The elements of π(π‘) will depend on the boundary conditions (5) - (9) and those of π½ on the initial conditions (4). We now need to adopt a time-marching method to solve (21). Popular choices for 1-D problems, such as the Implicit Euler and Crank-Nicolson schemes, become inefficient in higher dimensions, leading to large sparse systems that are a lot more expensive to solve than the small (typically tridiagonal) ones in the 1-D case. ADI-type splitting schemes are thus the most popular choice in 2-D and 3-D. However, standard (non-splitting) schemes can still be competitive for 2-D problems if a fast, sparse direct solver is used. This is especially true for the vanilla option pricing problem as the coefficients are time-independent and the matrix factorization step only needs to be performed once (or a few times). We employ one such method not often used in finance, namely the BDF3 (or 4-Level Fully Implicit) scheme. Our main workhorses though will be two popular ADI schemes which we briefly present first. For a detailed review of
ADI methods for PDEβs with mixed derivatives in finance, the reader is referred to [3]. The first step for all such methods is to decompose π¨ in (21) into three submatrices: π¨ = π¨ + π¨ + π¨ . (22) π¨ contains all terms stemming from the discretization of the mixed derivative term in (2), (3), i.e., all terms in (15) including π ππ£ as a factor. π¨ and π¨ contain all the terms corresponding to the discretized derivatives in the S -direction and v -direction respectively. The source term π π π is evenly distributed between π¨ and π¨ . By virtue of our 3-point central discretizations for the convection and diffusion terms, π¨ and π¨ are tridiagonal matrices . We split the vector π(π‘) and function π(π‘, π½) from (21) accordingly as π(π‘) = π (π‘) + π (π‘) + π (π‘) and π(π‘, π½) = π (π‘, π½) + π (π‘, π½) + π (π‘, π½). We will use a uniform temporal grid which is defined by the points π‘ π = π β βπ‘, 0 β€ π β€ ππ, βπ‘ = πππ . Let π be a real parameter which will control the exact splitting. We now outline our two main schemes, chosen for their optimal combination of stability, accuracy and inherent oscillation-damping properties [7]. Hundsdorfer-Verwer (HV) scheme { π = π πβ1 + βπ‘πΉ(π‘ πβ1, π πβ1 ), π π‘ππ 1π π = π πβ1 + πβπ‘ (πΉ π (π‘ π, π π ) β πΉ π (π‘ πβ1, π πβ1 )) (π = 1,2), π π‘πππ 2 & 3πΜ = π + βπ‘ (πΉ(π‘ π, π ) β πΉ(π‘ πβ1, π πβ1 )) , π π‘ππ 4πΜ π = πΜ πβ1 + πβπ‘ (πΉ π (π‘ π, πΜ π ) β πΉ π (π‘ π, π )) (π = 1,2), π π‘πππ 5 & 6π π = πΜ Modified Craig-Sneyd (MCS) scheme { π = π πβ1 + βπ‘πΉ(π‘ πβ1, π πβ1 ), π π‘ππ 1π π = π πβ1 + πβπ‘ (πΉ π (π‘ π, π π ) β πΉ π (π‘ πβ1, π πβ1 )) (π = 1,2), π π‘πππ 2 & 3πΜ = π + πβπ‘ (πΉ (π‘ π, π ) β πΉ (π‘ πβ1, π πβ1 )) , π π‘ππ 4πΜ = πΜ + ( β π)βπ‘ (πΉ(π‘ π, π ) β πΉ(π‘ πβ1, π πβ1 )) , π π‘ππ 5πΜ π = πΜ πβ1 + πβπ‘ (πΉ π (π‘ π, πΜ π ) β πΉ π (π‘ πβ1, π πβ1 )) (π = 1,2), π π‘πππ 6 & 7π π = πΜ Both schemes employ multiple intermediate steps to advance the solution from π½ πβ1 to π½ π . The HV scheme starts with a forward Euler (predictor) step (1), followed by two unidirectional implicit (corrector) steps (2 & 3) which serve to stabilize the explicit first step. Then a second predictor step (4) is followed by two more implicit corrector steps (5 & 6). The MCS scheme has an identical structure except for the double second predictor step (steps 4 & 5). The implicit steps require the solution of tridiagonal systems which we solve efficiently with LU decomposition. We use the HV scheme with π = 1 β β2 2β (which we shall refer to as HV1) and π = 1 2β + β3 6β (HV2). It was conjectured in [3] that HV1 is only conditionally stable (but more accurate), and HV2 unconditionally stable (and less accurate). For the MCS scheme we use π = 1 3β , recommended in [8] as an optimal value based on stability analysis and experiments. Regardless of the value of π, both schemes are second-order. We note that despite proven unconditionally (von Neumann-) stable, these schemes do not always sufficiently damp local high-frequency errors caused by discontinuities in the initial conditions. This may result in spurious oscillations and reduced order of convergence; see for example [3, 9]. In this case a technique known as Rannacher time-stepping can be used to palliate the issue. This involves using a different scheme for the first time-step (which is divided into two equal sub-steps), one that can successfully damp oscillations and is usually first-order (typically the Euler Implicit scheme). This is not strictly true for matrix π¨ because the one-sided formula (12) used for the π£ = 0 boundary involves one more point off the diagonal . For a more robust alternative that could conceivably provide smoother inputs to the (gradient-based) optimizer, we look at the third-order BDF3 scheme. Although not a typical choice, it is nonetheless simple to implement and has good stability (it is almost A-stable in an ODE sense) and oscillation damping properties. To solve the resulting systems, we use the Eigen C++ matrix library that offers simple interfaces to several direct sparse system solvers. The fastest one for the present system structure seems to be the UMFPACK solver, which we used for our experiments here. The scheme simply amounts to replacing the time derivative π½ β² (π‘) in (21) with a one-sided, 4-level backward finite difference expression. The discretized version of (21) then looks like
116 π½ π β 3π½ πβ1 + 32 π½ πβ2 β 13 π½ πβ3 βπ‘ = π¨π½ π + π π = π(π‘ π, π½ π ), (25) and the values π½ π at time level n are calculated given the values at the previous three time-levels as
116 π½ π = 3π½ πβ1 β 32 π½ πβ2 + 13 π½ πβ3 + βπ‘(π¨π½ π + π π ). (26) Since values are required not only from the previous time-level (like the ADI methods), but also from two levels before that, we must use some alternative scheme for the first two steps of the integration. We use the first-order Implicit Euler (IE) scheme and the second-order BDF2 scheme. The IE scheme is given by π½ π = π½ πβ1 + βπ‘(π¨π½ π + π π ) and requires the factorization of π¨ πΌπΈ = (π β βπ‘π¨) . The BDF2 scheme is given by π = 2π½ πβ1 β 0.5π½ πβ2 + βπ‘(π¨π½ π + π π ) and requires the factorization of π¨ π΅π·πΉ2 = (1.5π β βπ‘π¨) . In order to improve accuracy for the first time-step, we employ Richardson extrapolation like this: we first use the IE scheme for 4 sub-steps of size βπ‘/4 to get the values π½ at the end of the first time-step. We then repeat, this time using 2 sub-steps of size βπ‘/2 to obtain π½ and get the final composite values for the first time-step as π½ = 2π½ β π½ . Note that this requires 2 matrix factorizations, corresponding to π¨ πΌπΈ with βπ‘/4 and βπ‘/2 . To get the values π½ at the end of the second time-step we use the BDF2 scheme. In total, the present implementation requires 4 expensive factorizations which add a substantial upfront computational cost. Let us loosely define computational efficiency (CE) as the accuracy achieved per unit CPU time. A PDE-based solver cannot match the CE of semi-analytical solutions, such as those available for the Heston model. We therefore need to look into ways of improving the CE of our set-up. Here we consider grid construction, smoothing of the initial conditions and Richardson extrapolation.
Our domain is semi-infinite (or infinite for X in equation (3)), so in practice the grid needs to be truncated at some point. If the grid does not extend far enough then the imposed boundary conditions will not hold exactly true and forcing them on the solution will introduce some error. If the grid extends further than it needs to then the grid step sizes will be larger for the same number of points, resulting in less accurate finite difference approximations. There is no obvious way to determine the truncation limits, so here we make the empirical choices below. Note the dependence of the limits on the model parameters, which means that the grids used will be different for each objective function evaluation (based on the parameters set by the optimizer each time).
S-direction
For the S -grid, we truncate to the right at π πππ₯ = π (ππ(πππ₯(πΎ,π ))+ ππ ππ π‘ βπ) , where we set M = 5 and π ππ π‘ = 0.5(βπ£ + βπ£ πΏ ) . We then set: πππ₯ < 20πΎ. This choice leads to good solution accuracy overall, but for extreme model parameter regimes and benchmark calculations we additionally multiply π πππ₯ by a safety factor of 2 to 3. For the left boundary and for equation (2), we set π πππ =
0. For equation (3) in
π = ln(π) , we truncate at π πππ = ln(πππ(πΎ, π )) β ππ ππ π‘ βπ , where we set M = 6. We then further require that π πππ β€ πΌπΎ , where πΌ is some constant. We normally set πΌ = 0.1 but for high accuracy we recommend πΌ β€ 0.025. v-direction We set π£ πππ = 0 , i.e., we do not truncate the left boundary. To set an appropriate right boundary, we note that for π β β , π£ π‘ follows an Inverse Gamma distribution (see Appendix B). Given the distribution we can then set π£ πππ₯ = π£ ππππ‘ (π) = πΉ β² (π) (27) and πΉ β² is the inverse cumulative (Inverse Gamma) probability function. We find that a value of between β5 and β6 is necessary for accurate valuations. For short-dated options an empirical fraction of (27) can be used whenever π β π < 1 . Alternatively, one can numerically calculate the exact distribution β and thus π£ ππππ‘ (π) β for each expiration. This is described in Appendix B and is used for our experiments in Sec. 4 with β6 . We finally note that typically it will be π£ πππ₯ β« π£ . This observation alone necessitates the use of a non-uniform grid, described next. Computational efficiency can be improved significantly and any problems due to discontinuities mitigated, with a grid that concentrates more points where theyβre needed.
We employ a well-known one-dimensional grid-generating (stretching) function based on the inverse hyperbolic sine, which satisfies certain criteria for use with finite difference methods. The interested reader is referred to Vinokur [10]. The same function but in slightly different form is often used in the financial literature, see for example Tavella & Randall [11] and
In βt Hout & Foulon [3]. The grid in the S -direction is given by: π π = π πππ + πΎ (1 + sinh (π π ( πππ β π π )) sinh(π π π π )β ) (28) for π = 0,1, β¦ , ππ , where K is the strike (or more generally the desired clustering point) and π π , π π are free parameters. π π represents the percentage of total points that lie between π πππ and K and π π controls the degree of non-uniformity. We set π π = 4.5 which corresponds to moderate non-uniformity and generally results in low error profiles across the moneyness spectrum. Given π π , π π can be set so that the grid goes up to π πππ₯ using: π π = ln((π΄ + π π π ) (π΄ + π βπ π )β ) 2π π β , where π΄ = (π πππ₯ β πΎ) (πΎ β π πππ β ) . We then make sure that the strike falls exactly on a grid point by making a further slight adjustment to π π : we find π πΎ = βπ π ππβ and then reset π π = π πΎ /ππ . Finally, we use the same approach for generating the X -grid in equation (3). For the v -direction we again use the same grid-generating function: π£ π = π£ (1 + sinh (π π£ ( πππ β π π£ )) sinh(π π£ π π£ )β ) , (29) for j = 0, 1, β¦ , NV , which clusters points around π£ . Since π£ πππ₯ β« π£ we set π π£ = 8.5 which is as non-uniform as we can get before CE starts dropping. We first set π π£ so that the grid goes up to π£ πππ₯ : π π£ = ln((π΄ + π π π£ ) (π΄ + π βπ π£ )β ) 2π π£ β , where π΄ = π£ πππ₯ π£ β 1 . If we wanted to place the strike in the middle between grid points we would use π π = (π πΎ + 0.5)/ππ instead. We then find π π£ = max(βπ π£ ππβ, β0.2ππβ, 6) and reset π π£ = π π£ ππβ to ensure that π£ lies exactly on a grid point. When the input ππ is low and/or π£ πππ₯ β« π£ , then the above π π£ adjustment results in the last grid point now falling short of π£ πππ₯ . In such cases we keep adding points using (29) until π£ π β₯ π£ πππ₯ , which means that the final grid size will be ππ β , with ππ β β₯ ππ . Typically, ππ β will be up to 50% higher than the input ππ . An alternative construction that seems to have some advantage over the one just described, is a hybrid one, having the narrow (but most important) zone around π£ uniformly-spaced and the rest non-uniform. Haentjens & In βt Hout [12] propose one way of constructing such a grid in the S -direction for the solution of the Heston PDE. Here we use our first construction above as the base, to determine π π£ and ππ β . The segment (
0, 2π£ ) is then made uniform with step Ξπ£ π = π£ π π£ β . We then use the simple stretching function π£ π = π ππ sinh ( π π£ ππ ππ ) sinh(π π£ )β , π ππ = ππ β β 2π π£ + 1 , π ππ = π£ πππ₯ β 2π£ + Ξπ£ π (30) to generate the non-uniform part, choosing π π£ so that the first step is equal to Ξπ£ π . This can be easily achieved with any one-dimensional root-finding method. The two v -grid constructions generally result in comparable performance. But when used with Richardson extrapolation (Sec. 2.4.3), the second (hybrid) variant is always preferable. We thus use the hybrid construction for all the numerical experiments of Sec. 4. Whenever there is a discontinuity at some point in the initial conditions, it is usually a good idea to apply some sort of averaging for that point using the value(s) of adjacent point(s). That is effectively to smooth out the discontinuity (in this case located at the strike K ) before solving the PDE. The reason is that such discontinuities increase the solution error. To this end, here we just replace the (zero) initial condition values along the π = πΎ line of the grid (remember we made sure that there is a grid point π π πΎ on K ) with a simple average over nearby space as proposed in [13]. For vanilla options this amounts to setting: ππππ‘πΆπππ π πΎ ,π = 0.25βπ π π πΎ +1 β π π πΎ β1 , for j = 0, 1, β¦ , NV , where we have βπ = π π πΎ +1 β π π πΎ for calls and βπ = π π πΎ β π π πΎ β1 for puts. Richardson extrapolation (RE) can significantly increase accuracy for many problems adding only a small computational overhead. It simply involves calculating solutions based on two different grids (either spatial or temporal, usually with grid-step sizes ratio of
Here we apply it on the spatial level as follows: for a given resolution
ππ Γ ππ , we first generate a (
ππ/2 Γ ππ/2 ) grid and calculate an option price on it, π πππππ π . Then using the same grid parameters ( a S , b S ) and ( a V , b V ), we generate a ( ππ Γ ππ ) grid and use it to calculate π ππππ . Given that our discretization is full second-order in both S and v , we can then calculate the extrapolated price as π π πΈ = π ππππ + π πππππ π . Note that the fine grid will contain all the coarse gridβs points and add new ones in between. It is important that the relative location of the strike K is the same for the two grids, as is indeed the case (both have points exactly on the π = πΎ line). The main advantage is that while the computational cost increase is merely 25%, the accuracy is typically improved by 1-2 orders of magnitude (depending on the resolution used and the model parameters). RE works very well when sufficiently fine grids are used and not so well when the grids are too coarse (in which case it may well give worse accuracy than the single evaluation). This is because the To guarantee that the solution around π£ is always adequately resolved, we make sure that there is a minimum number of allocated grid points up to π£ , at least 20% of the total and no less than 6. premise for RE is that the two solutions are in the asymptotic range, i.e., that the observed order of convergence for the grids used is (very close to) the theoretical one. Down to the lowest resolution ( ππ Γ ππ ) = (40 Γ 20) used in our experiments, weβve found RE to clearly outperform the single evaluation in terms of CE. RE is also less effective for 2-D and 3-D problems when non-uniform grids with different stretching functions for each dimension are used (as is the case here). We find that this is effectively countered with the use of the hybrid v -grid that makes the grid in the v -direction uniform in the region of interest. This helps to regularize convergence, which in turns leads to improved RE performance. Discontinuities and/or singularities in the initial or boundary conditions will also often cause the observed order of convergence to be less than the theoretical one (and make convergence overall erratic), again reducing the effectiveness of RE. If those can be treated somehow, then convergence order is restored and RE performance improved. This is one more reason for applying the smoothing procedure described in Sec. 2.4.2.
3. Calibration
The main goal of the present work is to fit the GARCH diffusion model to a market of options. Some people choose to fit to option prices and others to the implied volatilities (IVs). We are strong proponents of the second approach.
IVβs are a natural way to regularize a set of option prices -- which can range from $0.05 to hundreds of dollars.
IVβs are the same order of magnitude across all the options . For SPX and other broad- based indices, using IVβs will also weight higher the influence of deep out-of-the money puts. Given that such options are a difficult regime for models (especially diffusions) to fit well, we like this property as well. It stresses an area where models have difficulty. Specifically, we try to fit the model to the option data by defining the following objective function to be minimized:
π πππΈ πΌπ = β 1π β(πΌπ ππππππ β πΌπ πππππππ‘ ) = π(π£ , π£Μ , π , π, π, ππ, ππ, ππ), (31) where N is the number of options we wish to include in the calibration . We calibrate to two SPX option chains, denoted Chain A and Chain B. Chain A used 246 SPX option quotes from Mar 31, 2017, filtered from quotes and IVβs calculated by the CBOE. The data and notes for that are found in Appendix C. For the optimization we use tools available in popular software . We test two local optimizers, Excelβs Solver tool which is based on the Generalized Reduced Gradient (GRG) method and
Mathematicaβs FindMinimum function which is based on an interior point method. We also use
Mathematicaβs NMinimize function , based on the global optimization Differential Evolution algorithm. All routines accept constraints which we impose on the model parameters in (31) in a way so that they encompass all plausible values. To work with Excel and Mathematica , we build a dll exporting a function that returns the RMSE IV taking just the PDE engineβs configuration as inputs. The function then reads the option chain data from a file, prices the options and evaluates (31). This is readily parallelized at the chain level, distributing the N options across all available CPU cores. We apply some basic load balancing since resolution (and thus calculation time - roughly proportional to ( ππ Γ ππ Γ ππ) ) may vary, as we discuss next. This excludes options with very small market prices which, as described in Sec. 3, will usually be priced with higher resolution than the nominal one input for the calibration. Obviously, NS, NV and NT (as well as the rest of the PDE engineβs configuration like choice of scheme, etc) are kept constant throughout a calibration (except when a negative value is detected as explained next). While a more customizable solution integrated with the PDE engine would likely be made to converge faster, we wish to keep things simple here and focus mostly on the PDE engine. Calling the function from Mathematica is trivial using the .NET/Link. For options of different expirations to be priced with similar accuracy, we need to have the number of time-steps ππ increasing with the expiration T . At the same time, the initial period of the valuation (close to the discontinuous initial conditions) always requires a minimum ππ to be resolved adequately. We roughly satisfy these requirements by taking the nominal ππ input in (31) to be the number of time steps per year for options with π > 1 , i.e., we set ππ πππ‘πππ = βππ β max(π, 1)β . We also find it is important to ensure that some minimum spatial resolution is used for the far out-of-the-money options in the chain, since those are more likely to incur higher relative pricing errors. More specifically, whenever the market value of an option is less than 0.5% of the asset spot π , (ππ Γ ππ) πππ is set to (120 Γ 70) which is then gradually increased to (
400 Γ 100) for market values of 0.01% of π or lower. W eβve found t hese empirical choices lead to better efficiency in terms of obtaining more accurate fitted parameters faster. As was explained in Sec. 2.2, the present discretization allows for negative option values by design. In general, those occur when the resolution is too coarse (and thus the accuracy too low) and the correlation coefficient π strongly negative. In practice we found such occurrences relatively rare under reasonable resolutions. With a local optimizer and when a previous result is used as the starting point, it is not unusual to complete a calibration involving tens of thousands of individual evaluations without a single negative value occurring. Such occurrences become even less frequent if one applies the transformation π = ln(π) , i.e. discretizing and solving equation (3) instead of equation (2). When negative valuations do occur during a calibration, the implied volatility cannot be calculated. In such cases we simply repeat the failed option valuation using our most robust (but least efficient) configuration, which involves switching to equation (3) and using the BDF3 scheme. We do so repeatedly, if required, using gradually increasing resolution until a positive value is returned. This βbrute - forceβ approach can occasionally slow down a calibration, mostly when a global optimizer is used. On the other hand, it βautomaticallyβ ensures that the option is priced accurately, which wouldnβt be the case if we used restrictions on the grid steps and/or added some sort of artificial diffusion aiming for an M-matrix (as was discussed in Sec. 2.2). Finally, since a valuation may just happen to be positive at (π , π£ ) (i.e., π > 0) , but still go significantly negative in the vicinity (and thus be inaccurate overall) , the naΓ―ve check of π > 0 is not sufficient. Instead, we check for negative values at all grid points within 10% of the strike in the S -direction and 50% of π£ in the v -direction and discard any positive valuation π if a negative value of magnitude more than 1% of π is detected. In general, if the model is to produce a decent f it to the market IVβs (and thus prices), then as the optimizer homes in on the optimum parameter set, the chances of a model price being that close to zero and thus susceptible to this problem are very low. Given our KBE PDE solver, the first and most obvious approach to evaluating (31) is indeed to price each of the N options separately, i.e., solve N PDEβs for each objective function evaluation.
Each PDE is solved on a different grid based on
πΎ, π and π , as described in Sec. 2.4.1. This is also the most general strategy since it can be used if we want to include options other than vanillas in the calibration exercise. The PDE engine could easily handle American or barrier options for example, at no extra cost. We shall refer to this as Approach I hereafter. Our purpose here though is to calibrate to vanilla options, in which case we can make use of the scaling (MAP) property introduced in Sec. 1.1. This means that only one PDE solution is sufficient to provide all option prices (for all different strikes) in an expiration bucket π π , π = 1, β¦ , π πΈ , where π πΈ is the number of different expirations included in the calibration. We choose to price one put option per π π and in particular the one that is furthest out-of-the-money (with the lowest strike price K ). This put necessitates an S -grid with the highest required π πππ₯ /πΎ ratio ( see Sec. 2.4.1.1) for our data . The prices of the rest of the puts are then readily found via the scaling relation and interpolation on the S -grid. The call prices are obtained from the put price at the corresponding scaled π via put-call-parity. Overall, It also ensures there are enough grid points where necessary for smooth higher order interpolation of the calculated option price at (π , π£ ) , and helps to avoid negative values by enforcing some minimum accuracy. W e include puts that are further out-of-the-money than the calls. this strategy (henceforth referred to as Approach II) requires π πΈ PDE solutions for one evaluation of (31) and the actual N option prices are extracted from those. Naively thinking, one may expect this to result in π/π πΈ times faster computation compared to Approach I (which would represent a 30-fold increase in the case of a chain with 240 options and 8 distinct expirations for example). In the previous section we described why for some options we want to use higher spatial resolution ( ππ Γ ππ) . Each π π bucket may include one or more such options, seen as weak links in terms of computational efficiency. With Approach I this means that about 80% of the N PDEβs can be solved very fast on coarse grids , while the rest (a few options in each π π bucket) will be priced on a much finer grid and take longer. Allocating 240 PDE solvers across different CPU cores amounts to reasonable medium-grain parallelism and allows for decent load-balancing. With Approach II on the other hand these advantages are lost: if we want all the extracted option prices to be of equivalent accuracy to when calculated individually (as in Approach I), we need to account for the weakest link in each case. If each π π includes at least one option that requires a fine(r) grid, it follows that all π πΈ PDE solvers need to use some overridden (high) resolution (as per previous section). This also means that now we may well have more CPU cores available than parallel tasks, say π πΈ = 8 and π πππππ = 10 , leaving some of the processing power unutilized. It may also lead to bad load balancing if for example one solver uses higher resolution that the rest. For this reason, we lower the maximum enforced (ππ Γ ππ) πππ from (
400 Γ 100) for Approach I to (2
00 Γ 80) for Approach II. The lesser accuracy that this implies for the few deep out-of-the-money options is offset by the fact that now all option prices are extracted from fine grids (as opposed to only a few with Approach I). As we will see next, this leads to calibrated model parameter accuracy as high as under Approach I. For the ADI schemes there is an alternative strategy and that is to parallelize at the PDE solver level. This is because grid lines can be updated simultaneously during both the explicit and implicit steps. Our brief tests with 8 or 10 cores/threads show that even with basic OpenMP instructions, a parallel efficiency of 80% is readily achievable this way, resulting in similar calibration times to our main, chain-level parallelization approach.
4. Numerical experiments and results
We now analyze the performance of the PDE engine and calibrator and present results using our two sample SPX option datasets, 246-option Chain A representing the 2017 low-volatility market, and 68-option Chain B from the higher volatility environment of 2010. Chain A data are given in Appendix C. The timings were taken on a 10-core Intel i9-7900X PC; the code was written in C++ and compiled in VS2013. We perform most of our tests using Approach I of the previous section (i.e. pricing each option separately). We do so since it is obviously preferable to assess the behavior/performance of the PDE engine over a sample of 246 or 68, rather than 8 or 7 individual option evaluations.
The ADI methods, as expected, prove to be more efficient than the BDF3 scheme and can be depended upon for successful and fast calibrations. At an individual option pricing level though we found they are less robust as they are not immune to spurious oscillations, mostly in the delta and gamma of the solution. The most likely offender is the HV1 scheme (which also proves to be the most efficient overall). We found this problem much rarer with the MCS scheme for the v anilla payoffs weβre dealing with here . The fully implicit BDF3 scheme demonstrates superior damping properties; we have been unable to reproduce a single case of this problem in extensive testing. In practice, with all reasonable (ππ, ππ, ππ) the ADI schemes are oscillation-free as well. Even when trying to force the issue (using some low ππ ππβ ), our tests show that any such mild oscillations present in the individual solutions do not usually prevent optimizer convergence (but they do slow it down). We note that this is a more We will see in the next section that a resolution of (
ππ Γ ππ) = (60 Γ 30) is more than enough to obtain accurate IVβs for the bulk of the options that are not deep out -of-the-money. Wyns [9] shows that such problems are more common when the MCS scheme is used to price cash-or-nothing options, where the discontinuity in the initial conditions is much more severe. general feature: the less accurate the PDE solution is, the more steps are usually required for the optimizer to converge. Using too low NT for example, will cause the optimizer to βhuntβ more ; conversely, application of (spatial) Richardson extrapolation typically means convergence in fewer steps than when simple (non-extrapolated) solutions are used. The other relevant problem is that of the occasional negative values. The use of equation (3) results in non-negative solutions under moderate spatial grid resolutions in cases where this is impossible with equation (2) and any reasonable resolution. Moreover, we find that in such cases the BDF3 scheme does a better job than the ADI schemes. As an example, we mention in passing a particularly difficult case that arose during a calibration with the global optimizer. In that case ππ = 20 was enough for the BDF3 scheme to produce a smooth, non-negative price profile near the strike, whereas the ADI schemes needed
ππ > 20000 to achieve the same result (with the same spatial discretization).
In terms of the optimizers we tested, we found Excelβs Solver tool to be the best choice. The main reason is that it benefits from a good initial guess, whereas Mathematicaβs FindMinimum does not. The Solver will typically converge 2.5 to 5 times faster when the starting vector is not too far from the optimal. Despite both being local optimizers, we are confident that they can be used to find the true (global) minimum. Extensive testing using many different starting points shows that both converge to the same vector . Using Mathematicaβs NMinimize global optimization routine furth er confirmed the solution in every case we tested. NMinimize also served as a torture test for the PDE engine, exploring all corners of the parameter space. All schemes never failed to produce a valid price, though of course βdifficultβ parameter sets (and insufficient resolution) generally trigger the repricing mechanism. Even so, the total optimization time is not significantly affected in practice. The allowed parameter ranges we used for the tests were: β€ 0.50 , , , , β0.95 β€ π β€ 0 , which cover most market scenarios. To test the convergence behavior of the time-marching schemes, we fix the spatial resolution to (ππ Γ ππ) = (60 Γ 30) and calculate (time-converged) benchmark prices using the BDF3 scheme with
ππ = 12800 . We also apply spatial Richardson extrapolation . This way the spatial discretization error is low, but not negligible compared to the temporal error. Nonetheless, the two errors are found to be only weakly dependent, allowing the comparative performance of the schemes to be properly assessed. The prices are obtained via an objective function evaluation under Approach I, which means that every option in the chain is priced individually and under the resolution overriding / repricing rules described in Sec. 3. The pricing errors for various ππ are calculated as the differences from the benchmark prices and the RMSE is used as an indicator of the overall performance for each scheme. Figure 1 shows the results for the HV1, HV2, MCS and BDF3 schemes, plus the HV1 scheme with Rannacher time-stepping (hereafter referred to as HV1D). The points on the left correspond to practical ππ (and CPU times) while those on the right are included to better illustrate the asymptotic behavior . We plot the relative (as opposed to absolute) pricing errors, since those are more closely related to the errors in the implied volatilities and consequently the calibrated model parameters. The HV1, HV2 and MCS schemes display a linear relationship between RMSE and CPU time on the logarithmic scale. This reflects (and confirms) their theoretical second-order convergence and the fact that the execution time in their case is proportional to NT . The Implicit Euler damping step (which requires an expensive factorization of the full system matrix) introduces an upfront cost that lowers the efficiency of the HV1D scheme. The irregular first two points from the left of the HV2 curve for Chain A are an example of the repricing mechanism in action: the schemeβs accuracy is too low h ere, causing some options in the chain to fail the βnegative values testβ (which are then revalued with a different configuration). Finally, for the BDF3 scheme we have an upfront cost that just like with HV1D is due to the initial matrix factorizations, resulting in significantly reduced efficiency for practical ππ . At large ππ this effect is diluted, and the scheme is seen to confirm its theoretical third-order convergence. The (composite) RE solution combines solutions on the (ππ Γ ππ) and (
ππ/2 Γ ππ/2 ) grids. The rightmost points of the error curves in Figure 1 correspond to
ππ = 1600 for the ADI schemes and
ππ = 400 for the BDF3. Figure 1 . Computational efficiency of time-marching schemes. RMS relative temporal error vs CPU time for the pricing of Chain A (246 options, left) and Chain B (68 options, right). Model parameters from Table 1.
Overall, the HV1, HV1D and MCS schemes offer comparable performance. Based on the present and many more similar tests, the HV1 scheme may well be the best choice. Under typical usage, any spurious oscillations do not seem to significantly affect its calibration performance. Adding damping comes at a mild cost overall and the HV1D scheme in fact produces lower (relative) errors for Chain B . The MCS scheme is still the most robust of the considered ADI schemes and usually about as accurate as HV1/HV1D, but it is slightly more computationally expensive per time-step. The HV2 scheme is the worst-performing and we will not consider it here further. Finally, the BDF3 scheme perhaps fares better than expected, beginning to outperform the ADI schemes for relative errors below around 0.02% in the case of Chain A. It is interesting to note that all temporal discretization errors are about an order of magnitude smaller for Chain B, making the BDF3 scheme less competitive in this case. Overall, we (unsurprisingly) find it cannot match the best ADI schemesβ efficiency for practical accuracy goals. β Calibration (Approach I)
We now test the performance of the PDE engine in terms of the end-result, the calibrated modelβs parameters. Benchmark values are given in Table 1, calculated using the BDF3 scheme and a resolution of (ππ Γ ππ Γ ππ) = (800 Γ 200 Γ 50) using spatial RE . Table 1.
Benchmark calibrated model parameters π£ π£Μ ο«ο οΈο ο²ο Option chain A (Mar 31, 2017) 0.010935 0.039139 5.3905 6.8997 -0.74579 Option chain B (Feb 1, 2010) 0.044816 0.088529 3.6695 5.0333 -0.79206
Figure 1 suggests that the HV1 scheme is more efficient than the MCS scheme for Chain A, but the two schemes perform very similar for Chain B. This translates to the calibrated model parameters. Table 2 shows that the parameters obtained for Chain A with the HV1 scheme are more converged for any ππ than those obtained with the MCS scheme and the CPU time required for the calibration is lower as well. We note that in this case the BDF3 scheme with ππ = 25 gives about the same accuracy as the The absolute errors (not shown) are slightly lower as well compared to the non-damped HV1, for both chains A and B. This may indeed suggest the presence of mild oscillations with the HV1 scheme (resulting in lowered accuracy) that are successfully suppressed by the added damping of the HV1D scheme. Weβve found cases where the Euler damping (Rannacher time -stepping) procedure of the HV1D scheme reduces but doesnβt eliminate spurious oscillations, whereas for the same cases the MCS scheme is oscillation -free (even without damping). In contrast, we have not seen any cases where the reverse is true. -7 -6 -5 -4 -3 -2 -1 R M S E secs HV1HV1DHV2MCSBDF3 -8 -7 -6 -5 -4 -3 -2 R M S E secs HV1HV1DHV2
MCS
BDF3 HV1 scheme with
ππ = 200 and requires the same time as well. Table 3 confirms the two ADI schemes produce similar accuracy for the Chain B parameters. The latter are also more converged for all ππ compared to the Chain A parameters in Table 2, reflecting the lower overall errors on the right plot of Figure 1. As in Sec. 4.2, we used a fixed spatial resolution of (ππ Γ ππ) = (60 Γ 30) with (spatial) RE and obtained time-converged calibrated model parameters using the BDF3 scheme with ππ = 400 . Table 2. M o del parameter convergence for calibrations to Chain A (246 options) as a function of temporal resolution for the HV1, MCS and BDF3 schemes. The exact (time-converged) parameter set is ( π£ , π£ Μ , π , π, π ) = ( ) . HV1 MCS BDF3 NT
25 50 100 200 25 50 100 200 25 50 π£ π£Μ ο«ο οΈο ο²ο -0.73655 -0.74338 -0.74511 -0.74554 -0.73753 -0.74197 -0.74476 -0.74545 -0.74586 -0.74570 CPU (mm:ss)
Table 3. M o del parameter convergence for calibrations to Chain B (68 options) as a function of temporal resolution for the HV1, MCS and BDF3 schemes. The exact (time-converged) parameter set is ( π£ , π£ Μ , π , π, π ) = ( ) . HV1 MCS BDF3 NT
25 50 100 200 25 50 100 200 25 50 π£ π£Μ ο«ο οΈο ο²ο -0.79120 -0.79205 -0.79226 -0.79232 -0.79121 -0.79204 -0.79226 -0.79231 -0.79227 -0.79233 CPU (mm:ss)
A careful inspection of the sequence of parameters obtained with the two ADI schemes with increasing (doubling) ππ , reveals very close to second-order convergence for ππ = 50 and above, indicating that the theoretical order of the time- discretization translates to the βfunctionalβ of the computation, i.e., the model parameter vector. This suggests the possible use of (temporal) Richardson extrapolation on the fitted parameters. The effectiveness of such an approach is demonstrated in Table 4, where the parameter vectors have been obtained as the composite (extrapolated) result of two successive calibrations using slightly different ππ . The vector from the first calibration can be used as the starting point for the second, reducing the time of the latter . Comparing for example the ππ = 100 calibrations in Table 2 with the (composite)
ππ = (40,60) calibrations in Table 4 for Chain A, one can see that the parameters obtained with the latter are significantly more converged, while the CPU times are lower as well. As always with RE, care should be taken not to use too low a resolution (in this case ππ ) . To guarantee good RE performance we recommend using ππ πππππ π β₯ 30 and ππ ππππ /ππ πππππ π β₯ 1.2. The effect of spatial RE on the calibrated parameters is shown in Tables 5 & 6. Here we use the BDF3 scheme with ππ = 50 so that the temporal discretization error is negligible. The benchmark parameters are from Table 1. Despite resolution being low, the parameters obtained with RE are not far Unless otherwise noted, all calibration CPU times listed in the present work are obtained with Excelβs Solver starting from the average values vector (π£ , π£Μ , π , π, π) = (0.05, 0.05, 5, 5, β0.7) . Also, the ADI schemes may not always sufficiently damp oscillations when
ππ/ππ is too low, possibly leading to erratic convergence and thus RE results. This is mainly an issue for the HV1 scheme. Applied on the PDE solution (option prices) as described in Sec. 2.4.3. from the benchmark values. For Chain A they are practically converged to 4 digits (the error is at most one point off in the fourth digit), while for Chain B they are somewhat less converged. In both cases the parameters obtained without the use of RE are significantly less accurate. Table 4. M o del parameter convergence using the ADI schemes and temporal Richardson extrapolation on the fitted parameters. (ππ Γ ππ) = (60 Γ 30) with spatial Richardson extrapolation. The exact (time-converged) parameter set is (π£ , π£Μ , π , π, π) = (0.010937, 0.039131, 5.3914, 6.9010, β0.74567) for Chain A and ( π£ , π£ Μ , π , π, π ) = ( ) for Chain B. Chain A (246 options) Chain B (68 options)
HV1 MCS HV1 MCS NT π£ π£Μ ο«ο οΈο ο²ο -0.74571 -0.74569 -0.74567 -0.74568 -0.79233 -0.79233 -0.79231 -0.79233 CPU (mm:ss) 05:07 09:07 07:07 11:20 00:52 01:08 01:00 01:20 Table 5.
Effect of spatial Richardson extrapolation on model parameter convergence for Chain A.
Spatial resolution π£ π£Μ ο«ο οΈο ο²ο (NS Γ NV) = 60 Γ 30 0.010914 0.039232 5.2805 6.8705 -0.74333 (NS Γ NV) = 60 Γ 30 w RE 0.010937 0.039131 5.3915 6.9009 -0.74570 Benchmark 0.010935 0.039139 5.3905 6.8997 -0.74579 Table 6.
Effect of spatial Richardson extrapolation on model parameter convergence for Chain B.
Spatial resolution π£ π£Μ ο«ο οΈο ο²ο (NS Γ NV) = 60 Γ 30 0.044811 0.088663 3.6064 5.0181 -0.79212 (NS Γ NV) = 60 Γ 30 w RE 0.044815 0.088479 3.6727 5.0330 -0.79233 Benchmark 0.044816 0.088529 3.6695 5.0333 -0.79206 So far, we have presented calibration tests where every option in the chain is priced separately. These tests demonstrated the efficiency of the PDE engine; calibrations using this Approach are already quite fast. We now test Approach II for the objective function evaluation, i.e., making use of the MAP property. Following the winning combination of Table 4, we choose the HV1 scheme, a resolution of (NS Γ NV) = (60 Γ 30) with spatial RE, as well as temporal RE on the fitted parameters (combining the results of two successive calibrations) with NT = (30,36). As expected, Table 7 confirms that the speed-up compared to Approach I is significant, especially for the largest Chain A . The parameter accuracy is at least as good. It is obvious that this approach works well in practice and effectively decouples the calibration time from the total number of options included. For either of our datasets a CPU time of less than a minute is needed to achieve a maximum While this nominal resolution is used by most of the PDE solvers within Approach I, in the case of Approach II, all (7 or 8) solvers really use higher resolution, as explained in Sec. 3.1. We note that the code was developed on an older CPU (4-core Intel i7-920, 2009) and not the new 10-core i9-7900X CPU used for the timings reported here. The MAP performance gains on the development CPU are almost double (6 Γ) those of the newer CPU, and the maximum calibration time still around 2 mins. relative (numerical) error of 0.05% for the obtained parameters. We stress that a judicious S -grid construction (low to moderate non-uniformity) is key for keeping the solution error profile low across the moneyness spectrum and make this approach work. As we already mentioned in Sec. 4.1, lower overall pricing accuracy leads not only to less accurate parameters but also (perhaps more importantly) to slower convergence of the optimizer. Table 7.
Comparison of pricing Approach I (one PDE solution per option) and Approach II (one PDE solution per expiration). The relative errors (compared to the benchmark parameters of Table 1) are shown in parentheses.
Chain A (246 options) Chain B (68 options)
Approach I Approach II Approach I Approach II π£ π£Μ ο«ο οΈο ο²ο -0.74569 (0.01%) -0.74581 (0.00%) -0.79233 (0.03%) -0.79209 (0.00%) CPU (mm:ss) 00:03:00 00:00:55 00:00:52 00:00:40 We now present detailed calibration results demonstrating the ability of the GARCH diffusion model to fit the option market and compare it to the popular Heston model. The Heston calibrations were performed using the present PDE engine and the resulting parameter vectors were then confirmed with independent calibrations using pricing via well-known Fourier integral representations. Figures 2 and 3 illustrate the market fit for each option expiration bucket in chains A and B respectively. As a first remark we can say that the model is indeed able to capture a smile behavior in the short end and achieve an overall decent fit. (By βsmile behaviourβ, we mean that the IV curve has an evident minimum). This can be seen in the first two plots in Figure 2, but not in Figure 3 (the Heston model can be seen to capture the smile in both cases). The reason is that the strike point K* (where the GARCH diffusion modelβs IV curve turns up ) for the first two expirations in Chain B lies further to the right of the last market point included in the plot (at around K = 1225). Overall, both sample calibrations seem to indicate that the Heston model is more β flexible β , managing to fit the data better overall. The GARCH diffusion model fits on the other hand look more β rigid β . This is somewhat surprising and in contrast with the findings of Christoffersen et al. [1]. While our two-chain data set is tiny compared to theirs, we suspect the contrast in findings is due to our much wider (smile) moneyness coverage β as no downside puts are used in [1]. The apparent Heston model victory here comes with known problems. The very small obtained Feller ratios, π β‘ 2π π£Μ π β , (0.12 and 0.29) are well below one. Note that under S. Hestonβs 1993 model [14], R is the same under both P (physical) and Q (risk-neutral) model evolutions. In our experience, the Heston P-model estimates (from time series, using maximum likelihood) will typically have
π > 1.
There are some caveats to complaining about π : for either model, P-model parameter estimates are not trivially obtained because the latent volatility must either be proxied or jointly estimated. In addition, P-model time series estimates are typically quite sensitive to the inclusion or not of crash days like Oct. 19, 1987. To trivially adapt the PDE engine to the Heston pricing PDE: adjust the v -diffusion and mixed derivative coefficients π π£ and π ππ£ accordingly in (17) and (20). The only other change required is the choice of π£ πππ₯ . We believe th e βhitting rangeβ, β€ 1β , should be admitted in calibrations. However, once you find the optimal Q-estimate ratio in the hitting range, you are forced to ponder the implications. Indeed, the volatility distribution then develops an integrable divergence at zero, v = 0 becomes the most probable value, and repeated volatility hits on zero become possible. If the P-evolution model has β , arbitrage opportunities develop β at least in the idealized continuous-time world in which the models are constructed. Yet, finding a smile-calibrated Feller ratio in the hitting range is well-known to be a common occurrence [15]. Figure 2 . GARCH diffusion and Heston model implied volatilities by expiry for Chain A. The GARCH diffusion model parameters are given in Table 1. The parameters for the Heston model are π£ = 0.007316, π£Μ = 0.03608,π = 6.794, π = 2.044 and π = β0.7184 . GARCH diffusion RMSE IV = 1.68%, Heston RMSE IV = 1.28%. Heston Feller ratio = 0.12. To purposely exaggerate the Feller ratio issue, we also fitted only the first two shortest expirations in Chain A. Despite seemingly achieving a decent fit (Figure 4), the Heston model fitted π£ is practically zero and the Feller ratio is 0.05. The GARCH diffusion model on the other hand achieves a closer fit with mostly reasonable parameters and a volatility process that doesnβt hit zero. Nevertheless, our calibrated volatility-of-volatility π - values for the GARCH diffusion are also likely βtoo highβ, relative to typical P-estimates. In the case of Chain B (Figure 5) the Heston model does a slightly better (if IV K T= 21 days
MarketGARCHHeston 00.1 K MarketGARCH
Heston IV K
T= 77 days
MarketGARCHHeston 0
IV K
Market
GARCH
Heston
500 1000 1500 2000 2500 3000 IV T= 259 days
Market
GARCH
Heston 00.1
400 900 1400 1900 2400 2900
IV K
T= 441 days
Market
GARCH
Heston00.10.20.30.40.5 400 900 1400 1900 2400 2900 3400 IV T= 630 days
Market
GARCH
Heston K T= 994 days
Market
GARCH
Heston oscillatory) job, following the smile correctly (with the model K * very close to the market K * at around K = 1140), while the GARCH diffusion model K * is about 1180. Figure 3 . GARCH diffusion and Heston model implied volatilities by expiry for Chain B. The GARCH diffusion model parameters are given in Table 1. The parameters for the Heston model are π£ = 0.04576,π£Μ = 0.06862, π = 4.905, π = 1.525 and π = β0.7131 . GARCH diffusion RMSE IV = 1.19%, Heston RMSE IV = 1.01%. Heston Feller ratio = 0.29. IV T= 19 days
MarketGARCHHeston0.1
650 750 850 950 1050 1150 1250 IV K Market
GARCH
Heston 0.10.15
IV K
T= 46 days
MarketGARCH
Heston K T= 228 days
MarketGARCHHeston0.1
450 550 650 750 850 950 1050 1150 1250 1350 IV T= 319 days
MarketGARCH
Heston 0.10.15
450 650 850 1050 1250 1450 1650
IV K
Market
GARCH
Heston0.10.150.20.250.30.35 450 650 850 1050 1250 1450 1650 1850
IV K
T=1054 days
MarketGARCHHeston Figure 4 . GARCH diffusion and Heston model implied volatilities when fitted only to the first two expiries in Chain A. The GARCH model parameters are π£ = 0.008046, π£Μ = 0.02981, π = 10.93, π = 15.06 and π = β0.5669 . The Heston model parameters are π£ = 10 β6 , π£Μ = 0.02, π = 40.64, π = 5.65 and π =β0.623 . GARCH diffusion RMSE IV = 0.61%, Heston RMSE IV = 0.86%. Heston Feller ratio = 0.05. Figure 5 . GARCH diffusion and Heston model implied volatilities when fitted only to the first two expiries in Chain B. The GARCH model parameters are π£ = 0.03836, π£Μ = 0.07067, π = 11.96, π = 7.685 and π = β0.6871 . The Heston model parameters are π£ = , π£Μ = 0.06011, π = 21.3, π = 2.604 and π = β0.6637 . GARCH diffusion RMSE IV = 0.64%, Heston RMSE IV = 0.59%. Heston Feller ratio = 0.38. To summarize, at the outset, we knew both models have their nice properties and their issues. In particular, neither model handles extreme moves β in either the asset price or its underlying volatility β either naturally or well. However, although our dataset is very small, we were still surprised to see the Heston model achieve the better smile fits. Both models we βve calibrated so far are subcases of the more general power-law SV model ππ π‘ = (π π β π π )π π‘ ππ‘ + βπ£ π‘ π π‘ ππ π‘π , (32) ππ£ π‘ = π (π£Μ β π£ π‘ )ππ‘ + ππ£ π‘π ππ π‘π£ . The popular affine Heston model ( π = 0.5 ) has some well-known limitations, such as: (a) inability to capture steep short-term volatility smiles, (b) instability of the fitted parameters over recalibrations and (c) incompatibility of fitted parameters with those estimated from the P-world. The GARCH diffusion model ( π = 1 ) has received less attention in the literature, but our tests here indicate that it too suffers from (a). Regarding (b) and (c) more tests would be needed; it may well be that GARCH produces more stable parameters than Heston for example. In an effort to address these issues, researchers have introduced a variety of ideas. Again, naming some: (i) moving away from affine models, so using the more general p -model above (or other two-factor diffusion variations), (ii) randomizing the (latent) spot variance π£ , (iii) adding jumps, (iv) adding more stochastic factors, (v) using fractional Brownian IV K
T= 21 days
MarketGARCH
Heston 00.10.2 K T= 49 days
MarketGARCH
Heston
850 900 950 1000 1050 1100 1150 1200 K T= 19 days
MarketGARCH
Heston 0.1
750 850 950 1050 1150 1250 IV T= 46 days
MarketGARCH
Heston motion as the volatility driver. Among these extensions, (i) and (ii) are particularly simple to add to the present framework, so we will briefly explore them here. Our PDE engine can easily solve for either the GARCH diffusion, Heston, or any model βin betweenβ, i.e., with . It is in fact straightforward to let π float and have the optimizer decide its optimal value. As Figures 2 and 3 show, a value of π = 0.5 corresponds to more β flexible β , whereas π = 1 to more β rigid β fits. Intermediate values of π predictably yield fits that look like a mix between the two. A low value of π reduces the RMSE IV by locally enabling enough curvature to emulate the short-term market smiles and on the other hand it increases it as the (overly) flexible behavior persists for the middle expirations (Heston case). Overall though, at least for our datasets, lower RMSE IV is achievable with lower values of p . Our tests found an optimum π of 0.62 for Chain A ( RMSE IV = 1.23%), while for Chain B we found π = 0.59 ( RMSE IV = 1%). As a result, the improvement in the fit was not significant over that of the Heston model (relative RMSE IV reduction of 1% - 4%). This may not be telling the whole story though . Letβs accept that both Heston and GARCH diffusion are essentially misspecified for the cause. Because of that, the optimizer is forced into unnatural parameters (e.g. unrealistically low π£ and/or high ΞΎ ) to adapt to the observed short-term smiles. If those two models are misspecified, so will be the general p -model. Our conclusion is that any benefit from choosing a particular value of p for the general p -model, would be best assessed if some sort of short-term-smile-enabling extension was in place, leaving the diffusion model more freedom to fit the rest of the expirations. Normally this would involve adding jumps, but an easier way is available through randomization; see Mechkov [16] or Jacquier & Shi [17]. While we find the dynamic rationale of this approach perhaps not entirely clear, we also found it does succeed in βturning onβ (up) the short -term smiles. Since again the changes needed in our code to add this feature were minimal, we gave it a try. Table 8 presents summary RMSE IV results for Chain A. As we noted above, varying p when the model lacks the ability to match the short-term smile makes little difference. But when the smile is better accounted for (here via the randomization), the optimal p- model provides a more substantial improvement in the calibration fit over both the Heston and the GARCH diffusion model. This case is presented in more detail in Appendix D. We finally add that other non-affine two-factor variations can also easily be accommodated. As an example, we calibrated the Inverse Gamma SV model introduced by LangrenΓ© et al. [15]. As can be seen in Table 8, it placed between the Heston and GARCH diffusion models in terms of overall quality of fit ( RMSE IV for Chain B was 1.11%). We also note is that the fitted π£ and π£Μ for this model were quite close to each other for both datasets (not shown). Table 8.
Calibration fit for Chain A under variations of pure diffusion models.
Model
πΉπ΄πΊπ¬ π°π½ Power law models
Heston 1.28% GARCH diffusion 1.68% Optimal p- model 1.23% Heston - Randomized 0.96% GARCH diffusion - Randomized 0.94% Optimal p- model - Randomized 0.79% Inverse Gamma Vol model 1.53%
5. Conclusions
In this work we present a first (to our knowledge) full option calibration of the GARCH diffusion model using a PDE approach. The calibration is very fast and accurate (less than a minute on a modern PC) ameliorating the lack of a closed-form solution. This is accomplished with the use of an efficient yet βordinaryβ second -order finite difference based PDE engine. While here we calibrate to European vanilla options, the same pricing engine can be used with only minor modifications for fast calibrations to other types of options that can easily be handled in the PDE setting, e.g. American options or barriers. Other similar models can also easily be accommodated as our brief experiments of Sec. 4.6 show. In a small test with two SPX option chains, the smile fits with the GARCH diffusion model were inferior to the fits from the Heston β93 model. This differed from some prior literature such as [1]. The Heston fits come with very low values for the so-called Feller condition ratio, which leads to other issues. Nevertheless, we were surprised. Our more general contribution is showing closed-form solutions need not be such a strong criterion for model selection. Similar PDE engines can potentially handle various related models that are being largely ignored in practice and therefore allow a more informed choice for a particular trading area. For the future, it would be quite interesting to extend the solver to one that could handle bivariate jump-diffusions with similar high efficiency. Finally, the second author (Lewis) would like to stress that the first author (Papadopoulos) has done all the heavy lifting here: developing and implementing all the C/C++ solvers and their Excel and Mathematica interfaces. References [1]
P. Christoffersen, K. Jacobs and K. Mimouni, βVolatility dynamics for the S&P500: Evidence from realized volatility, data returns, and option prices,β
Review of Financial Studies , vol. 23, no. 8, pp. 3141-3189, 2010. [2] R. Zvan, P. Forsyth and K. Vetzal, βNegative Coefficients in Two Factor Option Pricing Models,β Journal of Computational Finance, vol. 7, no. 1, pp. 37-73, 2003. [3]
K. J. in 't Hout and S. Foulon, βADI finite difference schemes for option pricing in the Heston model with correlation,β International Journal of Numerical Analysis and Modeling, vol. 7, no. 2, p. 303β
Sahalia and R. Kimmel, βMaximum likelihood estimation of stochastic volatility models,β Journal of Financial Economics, no. 83, p. 413 β S. Ikonen and J. Toivanen, βEfficient numerical methods for pricing American options under stochastic volatility,β Numerical Methods for Partial Differential Equations, vol. 24, no. 1, p. 104β
K. J. in 't Hout and J. Toivanen, βApplication of Operator Splitting Methods in Finance,β arXiv:1504.01022 [q-fin.CP], 2015. [8]
K. J. in 't Hout and C. Mishra, βStability of the modified Craigβ
Sneyd scheme for two-dimensional convection βdiffusion equations with mixed derivative term,β Mathematics and Computers in Simulation, vol. 81, no. 11, pp. 2540-2548, 2011. [9]
M. Wyns, βCo nvergence analysis of the Modified Craig-Sneyd scheme for two-dimensional convection- diffusion equations with nonsmooth initial data,β IMA Journal of Numerical Analysis, vol. 37, no. 2, pp.
M. Vinokur, βOn One -Dimensional Stretching
Functions For Finite Difference Calculations,β NASA
Contractor Report 3313, Santa Clara, California, 1980. [11] D. Tavella and C. Randall, Pricing Financial Instruments, The finite difference method, New York: Wiley, 2000. [12] T. Haentjens and K. J. in 't
Hout, βADI schemes for pricing American options under the Heston model,β
Applied Mathematical Finance, vol. 22, no. 3, pp. 207-237, 2015. [13]
D. M. Pooley, K. R.Vetzal and P. A. Forsyth, βConvergence remedies for non -smooth payoffs in option pric ing,β Journal of Computational Finance, vol. 6, no. 4, pp. 25 -40, 2003. [14]
S. L. Heston, βA Closed -form Solution for Options with Stochastic Volatility with Application to Bond and
Currency Options,β Rev. of Financial Studies, vol. 6, no. 2, pp. 327 -343, 1993. [15]
N. LangrenΓ©, G. Lee and Z. Zili, βSwitching to non -affine stochastic volatility: A closed-form expansion for the Inverse Gamma model,β arXiv:1507.02847v2 [q-fin.CP], 2016. [16] S. Mechkov, β'Hot - start' initialization of the Heston model,β Risk, November 2016. [17] A. Jacquier and F. Shi, βThe randomised Heston model,β arXiv:1608.07158v2 [q-fin.PR], London, 2017. [18] Chicago Board Options Exchange, βVIX white paper,β 2010. [19] E. Γ inlar, Markov additive processes I, II., Z. Wahrscheinlichkeitsth. verw. Geb., 24:I:85-93, II:93-121, 1972. [20] T. Bollerslev and P. E. Rossi, βDan Nelson Remembered,β J. Business & Econ. Statistics, 13(4), pp. 361 -364, 1995. [21]
D. W. Peaceman and H. H. Rachford, βThe numerical solution of parabolic and elliptic differential equations,β SIAM Journal, no. 3, p. 28β
41, 1955. [22]
J. Douglas and H. H. Rachford, βOn the numerical solution of heat conduction problems in two and thre e space variables,β Trans. Amer. Math. Soc., no. 82, p. 421β R. Rannacher, βFinite element solution of diffusion problems with irregular data,β Numer. Math., no. 43, pp. 309-327, 1984.
Appendix A β Market prices of risk and a compatible real-world evolution
We have performed option chain calibrations after postulating that the risk-neutral (aka Q-measure) evolution has the GARCH diffusion form (1). Jumping immediately to a risk-neutral model is a common finance short-cut. More carefully, even given a target Q-model, one should begin with a compatible real-world (P-measure) evolution, and then move to the desired Q-measure evolution by a Girsanov transformation. In the presence of a deterministic stock dividend yield π π‘ , a P-measure evolution compatible with a Q-measure GARCH diffusion under this procedure has the form: ππ π‘ = (πΌ π‘π β π π‘ )π π‘ ππ‘ + π π‘ π π‘ ππ π,π‘π , ππ£ π‘ = π½ π‘π ππ‘ + ππ£ π‘ ππ π,π‘.π£ . Now (π π,π‘π , π
π,π‘π£ ) are a pair of correlated P-Brownian motions with correlation π , both (π, π) are identical under P or Q, and we used π π‘ β‘ βπ£ π‘ . Indeed, βno - arbitrageβ requires that the Q-evolution must be related to the P-evolution by the Girsanov substitutions ππ π,π‘π = ππ
π,π‘π β π π‘π ππ‘ and ππ π,π‘.π£ =ππ
π,π‘π£ β π π‘π£ ππ‘ . Under this implied change-of-measure, the variance-covariance structure of the SDE is preserved but the drifts may change. Financially, (π π‘π , π π‘π£ ) represent market prices of (equity, volatility) risk. The π functions are independent of any derivative asset but generally dependent upon (π‘, π π‘ , π£ π‘ ). The Q-evolution model then becomes ππ π‘ = (π π‘ β π π‘ )π π‘ ππ‘ + π π‘ π π‘ ππ π,π‘π , ππ£ π‘ = π½ π‘π ππ‘ + ππ£ π‘ ππ π,π‘π£ , where π π‘π = (πΌ π‘π β π π‘ )/π π‘ and π½ π‘π = π½ π‘π β ππ£ π‘ π π‘π£ . Fixing π½ π‘π = π π β π π π£ π‘ (where π π β‘ π π π£Μ π ) from our postulated Q-measure GARCH diffusion (1) still leaves a lot of freedom for the P-evolution. In this generality, the only remaining compatibility requirements under βno - arbitrageβ are: β’ πΌ π (π£ π‘ = 0) = π π‘ , since a stock holding would be instantaneously riskless, presuming a deterministic short-rate π π‘ . β’ the boundaries π£ = 0 and π£ = β should be unattainable in finite time by the P-measure π£ π‘ -process, since that is true of the postulated Q-measure process. However, the spirit of the model is that the P-measure evolution is also a GARCH diffusion (recall the origin of the name). So, let π½ π‘π = π π β π π π£ π‘ , with possibly different P-parameters. For example, letβs postulate i) a volatility-dependent equity risk premium πΌ π‘π = π π‘ + ππ£ π‘ , where π is a (positive) constant, and ii) π π£ is also constant. With those choices, our associated P-model GARCH diffusion is ππ π‘ = (π π‘ β π π‘ + π π£ π‘ )π π‘ ππ‘ + π π‘ π π‘ ππ π,π‘π , ππ£ π‘ = (π π β π π π£ π‘ ) ππ‘ + ππ£ π‘ ππ π,π‘π£ , where π π = π π and π π = π π + π π£ π . Now two additional parameters (π, π π£ ) need to be estimated and β P/Q compatibility β under our choices becomes a hypothesis to be tested. All that is outside our scope in this article. However, one expects π > 0 and π π£ < 0 for, say, SPX. Finally, note that with our choices and using (π₯ π‘ , π£ π‘ ) where π₯ π‘ = log π π‘ , both real-world and risk-neutral processes are MAPs (Markov Additive Processes). As discussed in the body, the MAP property leads to an βall -options-at- onceβ KBE solution for vanilla options; in addition, it also allows dimensional reduction by Fourier methods. Appendix B β Critical points for the stand-alone volatility process
In our model, the stand-alone volatility process evolves as ππ π‘ = π (π£Μ β π π‘ )ππ‘ + ππ π‘ ππ΅ π‘ , where π΅ π‘ is a Brownian motion. The same functional form holds under either measure (P/Q), although the numerical values of the parameters may differ. Let π(π‘, π£, π£ ) denote the transition probability density for the process; i.e. π(π‘, π£, π£ )ππ£ β‘ Pr(π π‘ β ππ£|π = π£ ) . Then, π£ π = π£ π (π, π‘, π£ , π , π£Μ , π) , the π -critical point for the associated distribution, is defined by β« π(π‘, π₯, π£ π ) ππ₯ = π , where . Having π£ π (for π close to 1) is useful for setting the π£ -grid (upper) truncation points for the full 2D process PDE solvers. Now π(π‘, π£, π£ ) is not known analytically, but solves the Fokker-Planck problem: ππππ‘ = β ππ½ππ£ , where π½(π‘, π£) β β π {π£ π} + π (π£Μ β π£)π , π£ β (0, β), and subject to the initial condition π(0, π£, π£ ) = πΏ(π£ β π£ ), using the Dirac-delta. Here π½(π‘, π£) is the probability current or flux. Mathematically, since both π£ = 0 and π£ = β are inaccessible to the process (with π£ > 0) , no boundary conditions are necessary in the continuum problem. Two relations.
It is easy to find that in the limit π‘ β β , π£ π‘ follows an Inverse Gamma distribution π(π₯) = π½ πΌ π₯ βπΌβ1 π βπ½/π₯ /Ξ(πΌ) , with shape parameter πΌ = 1 + 2π /π and scale parameter π½ = 2π π£Μ /π . Another easy relation for arbitrary t is the scaling identity: π£ π (π, π‘, π£ , π , π£Μ , π) = π£Μ Γ π£ π (π, π π‘, π£ π£Μ , π π , 1,1), which reduces the effective number of parameters by two. While the scaling relation was not used in the implementations, it was checked. Mathematica implementation.
When speed is not a factor, this full problem (solving the Fokker-Planck PDE, and calculating the critical point) is readily solved in Mathematica. Our short implementation is shown in Fig. 6. Even if you are not a Mathematica user, the syntax should be largely readable. The basic idea is to convert to new coordinates π₯ = log π£ and solve the resulting PDE problem using NDSolve. A uniformly-spaced x -grid with π₯ points is centered at π₯ = log π£ . The grid is truncated at Β±π βsigmaβsβ from π₯ , where one sigma equals πβπ . Numerical boundary conditions are taken to be (i) zero flux at π₯ πππ and (ii) a zero spatial derivative at π₯ πππ₯ . The initial condition is a lattice Dirac-delta, non-zero only at π₯ , which lies exactly on a node. C/C++ implementation.
The Mathematica implementation solves the above Fokker-Planck problem using 4 th order spatial discretization and the Method of Lines (MOL) via an ODE solver in time. This yields very accurate results but is slow (though we have not tried to port it to C++). For the tests presented in this paper we have opted for a more standard approach which we only briefly outline here: The discretization is based on uniform central second order finite differences and the Crank-Nicolson scheme with Rannacher time-stepping. Boundary conditions are the same as above. This approach works, apart from the far-left region of the grid where convection may dominate and result in oscillations/ negative densities. In such cases we locally introduce the 1 st order upwind scheme for the convection term. If negative densities are still produced, we try to bump up π π₯ . If all fails, we simply return π£ π from the stationary (Inverse Gamma) distribution. We also apply double spatial Richardson extrapolation (which should in theory result in 6 th order accuracy- if the upwind scheme is only used in areas where the density takes negligible values). All this results in accuracy even higher than Mathematicaβs , requiring CPU times of about 2-3 milliseconds with π π₯ = 800 and π π = 50 . Figure 6 . Mathematica code computing the stand-alone V-distribution critical points.
Appendix C β Data description
Option Chain A . End-of-day (EOD) SPX option data on March 31, 2017 was obtained from the
CBOEβs LiveVol service : βEnd -of- Day Option Quotes with Calcsβ.
These files record option quotes and CBOE calculated option implied volatilities (IVβs) at 15:45 New York time. This time is 15 minutes prior to the regular session close in both NYC and Chicago. From the CBOE (edited for brevity): β Implied volatility and Greeks are calculated off the 1545 timestamp, considered a more accurate snapshot of market liquidity than the end of day market. LiveVol applies a unified calculation methodology across both live and historical data sets to provide maximum consistency between back-testing and real-time applications. Cost of carry inputs (interest rates, dividends) are determined by a statistical regression
CriticalPointGarchPDE q , V0 , T , NX , vbar , kappa , xi , n1 , AG :Module X0, Xmin, Xmax, mu kappa xi ^ 2, omega, A,h0, h, dX, i, grida, gridb, grid, t, soln, cdf, critx, critv, critx99 ,Off NDSolve::eerri ;X0 N Log V0 ;Xmin X0 n1 xi Sqrt T ;Xmax X0 n1 xi Sqrt T ;dX X0 Xmin NX ;omega kappa vbar ;grida Xmin N Table i dX, i, 0, NX ;grida NX NX ;grid Join grida, gridb ;h0 X ? NumericQ : If Abs X X0 0.5 dX, 1 dX, 0 ;Clear soln, h, t ; Off NDSolve::eerr ;using zero flux condition at Xminsoln h . NDSolve t h x, t 0.5 xi ^ 2 x,x h x, t x omega E ^ x mu h x, t ,h x, 0 h0 x , 0.5 xi ^ 2 D h x, t , x . x Xmin omega E ^ Xmin mu h Xmin, t 0,D h x, t , x . x Xmax 0 , h , x, Xmin, Xmax , t, T , T , AccuracyGoal AG ,Method "MethodOfLines","SpatialDiscretization" "TensorProductGrid", "Coordinates" grid 1 ;A NIntegrate soln x, T , x, Xmin, X0 , AccuracyGoal AG ;cdf y ? NumericQ : A NIntegrate soln x, T , x, X0, y , AccuracyGoal AG ;critx99 y . FindRoot cdf y 0.99, y, X0, Xmin, Xmax ;critx y . FindRoot cdf y q , y, critx99, Xmin, Xmax ;critv E ^ critx;Return critv process β¦ . The cost of carry projected from these inputs is compared against those implied by the at-the-money options from each option expiry. If the rates differ significantly β and the option spreads for this expiry are sufficiently narrow β the implied rates replace the standard inputs. β¦β . Table 9 . SPX Option Implied Volatilities (%): Chain A (15:45, March 31, 2017).
Option expiration date (mm/dd/yy format) Strike
On any given data, option data can be quite voluminous and needs to be filtered, both to reduce calculational burdens and remove irrelevant noise. Indeed, the full March 31, 2017 data file contained 8399 option line items, which we filtered to the 246 items shown in Table 9. This was done by, first, focusing on traditional βthird Fridayβ expirations and then doing some strike filtering. The 15:45 SPX index value was π = 2367.94 (midpoint quote). We selected positive bid, out-of-the-money options, so the IVβs shown in the table are from puts when the strike
πΎ < π and otherwise calls. (The CBOEβs IV methodology is somewhat of a black-box, but it appears to be essentially put-call-parity preserving. See also [18]). For the first expiration we chose (Apr 21, 2017), the implied volatilities were smooth down to a strike πΎ = 1800 , chosen as a lower limit strike cutoff. For other expirations, the data looked smooth down to
πΎ = 500 , our cutoff for the remaining expirations. We imposed no upper strike cutoff. To achieve a rough balance between puts and calls, we selected put strikes at multiples of 100 and call strikes at multiples of 25. We believe this filtering retains the important characteristics of the full data sets. For short-term interest rates, we found U.S. Treasury debt asked yields on Mar 31, 2017 from the Wall Street Journal (WSJ) and use those as stepwise constants for each of our 8 expirations. That series was {0.00728, 0.00723, 0.00716, 0.00865, 0.00939, 0.01118, 0.01203, 0.01434} in expiration order. For the SPX dividend yield, we used a constant π = 0.0197 for all expirations, the WSJ-reported trailing 12-month SPX yield on the same date. These are not necessarily the cost-of-carry parameters used by the CBOE, the latter unavailable. The difference is unlikely to change the parameter fits or our conclusions in any way that matters. But if the reader is concerned about this small point, then take our combination of IVβs and cost -of- carryβs as one β possible β market data set (largely consistent with the 3/31/17 actual data) to which we fit various models. Option Chain B.
The second author (Lewis) collected (out-of-the-money) closing option quotes at the time (Feb. 1, 2010), using only options with positive bids.
IVβs were calculated from the bid -ask midpoint option price. Interest rates and a dividend yield was found from the WSJ as per Chain A.
Appendix D β Calibration for the randomized optimal p -model A randomized version of the Heston model was presented by Mechkov [16]. The basic idea is that instead of taking the initial value of the (latent) variance process π£ to be a known fixed value, we assume it is given by some distribution. It is a simple and appealing idea making for an easy extension to the present framework: The PDE solution automatically provides the option values for the whole range of possible initial variance values (corresponding to the grid in the v -direction). To find the "randomized" option price at the asset spot π , average the solution across the π = π grid line using the assumed distribution. As suggested in [16], it is reasonable to assume the latter should be of the same type as that of the processβ equilibrium distribution. For the GARCH diffusion model this would be the Inverse Gamma distribution. Our brief testing indicates that this choice yields the best results even for the randomization of the Heston model , so we will use it here to randomize the general p -model (32). We make the parameters of the distribution (shape πΌ and scale π½ ), as well as the power p of the model part of the calibration. The total number of parameters to be fitted is now seven. As Figure 7 shows, the overall quality of fit ( RMSE IV = 0.79% ) is considerably better than either that of the GARCH diffusion ( RMSE IV = 1.68% ) or the Heston model ( RMSE IV = 1.28% ), especially for the shortest expiration (see Figure 3). As already discussed in Sec. 4.5, the Heston model ( π = 0.5) fit implies unlikely dynamics. The optimal power of the general model was calibrated to π = 0.8 , slightly closer to the GARCH diffusion model . In contrast to the latter, here we see a steep smile captured for the first (3W) expiration, which is due to the randomization (and not the change from π = 1 to π = 0.8 ). We also note how now that the model is not bound (and thus stressed) by its inability to account for the short-term smile, the fitted volatility of volatility parameter falls to arguably more realistic levels ( π = 2) . This is much lower than GARCH diffusionβs calibrated val ue of π = 6.9 (Table 1), even after accounting for the scaling from π = 1 to π = 0.8 . A potential problem is seen with the calibrated correlation coefficient ( π = β1 ). Other similar tests also indicate that the randomization procedure tends to lead the optimizer to rather extreme correlation values. Why? From [17], the small-T asymptotic smile βexplosionβ , due to randomization, is symmetric The equilibrium (stationary) distribution of the square root variance process is a Gamma distribution, but at least for our datasets we found it actually performs worse than the Inverse Gamma as the randomizing (initialization) distribution for the Heston model. The model is much closer to the GARCH diffusion model in terms of the implied dynamics. For example, one finds that the Heston modelβs fit here implies a 41% probability that the long -run (risk-neutral) volatility is less than 1%, which is not very plausible. For both GARCH diffusion and the π = 0.8 model this probability is practically zero. The mean long-run volatility is about 20% for all three models. about π₯ πΎ = 0, where π₯ πΎ = log πΎ/π is the log-moneyness. The optimizer may be trying to compensate for the β unwanted β new tendency towards symmetry by pushing π towards -1. Excluding π from the calibration and fixing it to a more reasonable value ( π = β0.8 ) yields a very similar fit with RMSE IV =0.83% . This indicates that the much-improved fit does not depend strongly on such extreme π values. Nevertheless, this behavior (as well as needing a dynamical rationale) are issues for randomization. Figure 7 . Calibration of the general power-law model (32) for Chain A, randomizing π£ with an Inverse Gamma distribution. The calibrated model parameters are π£Μ = 0.0407 , π = 3.34, π = 2.00, π = β1.00 , the initial distribution shape and scale parameters (see Appendix B) πΌ = 1.05, π½ = 0.00124 and the model power π =0.801 . RMSE IV = 0.79%. IV K T= 21 days
Market
Randomizedp = 0.8 model IV K T= 49 days
Market
Randomizedp = 0.8 model00.1 IV K T= 77 days
Market
Randomized p= 0.8 model
IV K
T= 168 days
Market
Randomizedp = 0.8 model00.1
IV K
T= 259 days
MarketRandomizedp = 0.8 model
IV K
T= 441 days
Market
Randomizedp = 0.8 model00.10.20.3
IV K
T= 630 days
Market
Randomizedp = 0.8 model
IV K
T= 994 days