Arbitrary Number of Stencil Points

hpcquantlib1 Comment

The standard finite difference implementations of derivative pricing algorithms based on partial differential equations have a spatial order of convergence of two. Reason is that these implementations are using a three point central stencil for the first and second order derivatives. The three point stencil leads to a tridiagonal matrix. Such linear systems can be solved efficiently with help of the Thomas algorithm.

Higher order of spatial accuracy relies on stencils with more points. The corresponding linear systems can no longer be solved by the Thomas algorithm but by the BiCGStab iterative solver with the classical three point stencil as a very efficient preconditioner. The calculation algorithm for the coefficients of larger stencils on arbitrary grids is outlined in [1].

Before using higher order stencils one should first make sure, that the two standard error reduction techniques are in place, namely adaptive grid refinement around important points [2] and cell averaging around special points of the payoff at maturity or – equally effective for vanilla options – put the strike in the middle between two grid points. Let’s focus on the Black-Scholes-Merton PDE

\displaystyle \frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2\frac{\partial^2 V}{\partial x^2} + \left(r-q-\frac{\sigma^2}{2}\right)\frac{\partial V}{\partial x} -rV = 0 \ \wedge \ x= \ln S

after transforming it to the heat equation

\displaystyle \frac{\partial u}{\partial \tau}=\frac{1}{2}\sigma^2\frac{\partial^2 u}{\partial y^2} \ \wedge \ y = x + \left(r-q-\frac{1}{2}\sigma^2\right)\tau \ \wedge \ t=T-\tau \ \wedge \ u=e^{r\tau}V

to demonstrate the convergence improvements. The spatial pricing error is defined as the RMSE for a set of benchmark options. In this example the benchmark portfolio consists of OTM call or put options with strikes

k\in \{50, 75, 90, 100, 110, 125, 150, 200\}.

The other parameters in this example are

S_0=100, r=5\%, q=2.5\%, \sigma=20\%, T=1.

The Crank-Nicolson scheme with “more than enough” time steps was used to integrate the PDE in time direction such that only the spatial error remains.

convergence

As can be seen in the diagram above the spatial order of convergence increases from second order to fourth order when moving from the standard three point to the five point stencil. On the other hand solving the resulting linear systems takes longer now. All in all only if the target RMSE is below 10^{-5}, the five point stencil will be faster than the standard three point discretization.

The diagram below shoes the effective combination of a five point stencil and the Richardson extrapolation.

richardsonDiagrams are based on the test cases testHigerOrderBSOptionPricing and
testHigerOrderAndRichardsonExtrapolationg in the PR #483.

 

[1] B. Fornberg, [1988], Generation of Finite Difference Formulas on Arbitrarily Spaced Grids.
[2] Tavella, D. and C.Randall [2000], Pricing Financial Instruments The Finite Dierence Method, Wiley Series In Financial Engineering, John Wiley & Sons, New York

Advertisement

Andreasen-Huge Volatility Interpolation

hpcquantlib10 Comments

A few years ago Andreasen and Huge have introduced an efficient and arbitrage free volatility interpolation method [1] based on a one step finite difference implicit Euler scheme applied to a local volatility parametrization. Probably the most notable use case is the generation of a local volatility surface from a set of option quotes.

Starting point is Dupire’s forward equation for the prices of European call options C_t(T, K) at time t with strike K and maturity T is given by

\displaystyle \frac{\partial C_t(T, K)}{\partial T} = \frac{1}{2}\sigma_{LV}^2(T,K)K^2\frac{\partial^2 C_t(T, K)}{\partial K^2} - \left(r(T) -q(T)\right) K\frac{\partial C_t(T, K)}{\partial K}-q(T)C_t(T, K)

Define the normalized call price \hat{C}_t(T, \kappa) in terms of the discount factor P(t,T), the forward price F(t,T) and the moneyness \kappa as

\begin{array}{rcl} \displaystyle P(t,T) &=& \displaystyle  e^{-\int_t^T r(\tau)d\tau} \nonumber \\ \displaystyle  F(t,T) &=&  S_t  \displaystyle  e^{\int_t^T \left( r(\tau)-q(\tau)\right)d\tau} \nonumber \\ \displaystyle \kappa &=& \displaystyle \frac{K}{F(t,T)} \nonumber \\ \displaystyle \hat{C}_t(T, \kappa) &=& \displaystyle \frac{C_t \left( T, \kappa F(t,T) \right) }{P(t,T)F(t,T)} \nonumber \end{array}.

The Dupire forward equation for the normalized prices \hat{C}_t(T,\kappa) is then given by

\displaystyle \frac{\partial \hat{C}_t(T, \kappa)}{\partial T} = \frac{1}{2} \sigma_{LV}(T, \kappa F(t,T))\kappa^2 \frac{\partial^2 \hat{C}_t(T, \kappa)}{\partial \kappa^2} .

Rewriting this equation in terms of u=\ln \kappa amd \tilde{C}_t(T, u) = \hat{C}_t(T, \kappa) yields

\displaystyle  \frac{\partial \tilde{C}_t(T,u)}{\partial T} = \frac{1}{2} \sigma_{LV}^2 \left( T, e^u F(t,T)\right) \left( \frac{\partial^2 \tilde{C}_t(T, u)}{\partial u^2} - \frac{\partial \tilde{C}_t(T, u)}{\partial u}\right)

The normalized put prices \tilde{P}_t(T, u) are fulfilling the same equation, which can easily been shown by inserting the call-put parity into the equation above

\displaystyle \tilde{C}_t(T, u) = \tilde{P}_t(T,u) + 1 - e^u.

The numerical stability of the original algorithm [1] can be enhanced for deep ITM options by calibrating to calls and puts instantaneously. Also the interpolation scheme has a significant impact on the stability. This topic has been discussed in [2][3]. Using concentrated meshes along the current spot level for the finite difference scheme is of advantage for the stability and accuracy of the algorithm.

In order to stabilize the calculation of the local volatility function

\displaystyle \sigma_{LV}\left( T, e^uF(t,T)\right) = \sqrt{\frac{2 \frac{\partial \tilde{C}_t(T,u)}{\partial T}}{  \frac{\partial^2 \tilde{C}_t(T, u)}{\partial u^2} - \frac{\partial \tilde{C}_t(T, u)}{\partial u} }}

one should evaluate the first order derivative of \tilde{C}_t(T,u) w.r.t. time T using the fact that the derivative of the inverse of the matrix \bf{A}(t) is given by

\displaystyle \frac{\partial \bf{A}(t)^{-1}}{\partial t} = -\bf{A}(t)^{-1} \frac{\partial \bf{A}(t)}{\partial t} \bf{A}(t)^{-1}

As an example, the diagram below shows different calibrations of the Andreasen-Huge volatility interpolation to a SABR volatility skew at discrete strike sets

x \in \{ 0.02, 0.025, 0.03, 0.035, 0.04, 0.05, 0.06\} \vee x\in \{0.03, 0.035, 0.04\}

for the SABR parameter

\alpha = 0.15, \beta = 0.8, \nu = 0.5, \rho = -0.48, \text{fwd} = 0.03, T=20
sabr

The QuantLib implementation is part of the release 1.12.

[1] J. Andreasen, B. Huge, Volatility Interpolation
[2] F. Le Floc’h, Andreasen-Huge interpolation – Don’t stay flat
[3] J. Healy, A spline to fill the gaps with Andreasen-Huge one-step method

Adaptive SOR Method for Implied Volatility Calculation

hpcquantlibLeave a comment

In a recent blog contribution Fabien Le Floc’h [1] suggests to combine the adaptive successive over-relexation method [2] with an improved explicit approximate implied volatility formula [3] to calculate the initial guess. The implementation of both algorithms is straight forward.

A large set of OTM and ITM options together with different displacement factors has been identified to serve as a benchmark portfolio to compare the performance with the original QuantLib implementation. As can be seen in the diagram below the new method is – depending on the accuracy – two to three times faster than the current QuantLib implementation.

sorThe implementation is available here, the diagram above is derived from the test case testImpliedVolAdaptiveSuccessiveOverRelaxation in the class BlackFormulaTest.

[1] Le Floc_h, F (2017) Implied Volatility from Black Scholes Price

[2] Li, M. (2008) An Adaptive Successive Over-relaxation Method for Computing the Black-Scholes Implied Volatility

[3] J. Gatheral, I. Matic, R. Radoicic, D. Stefanica (2017), Tighter Bounds for Implied Volatility

Andersen-Piterbarg Integration Limits for the Time Dependent Heston Model

hpcquantlibLeave a comment

The normalized characteristic function \phi_t(z) of a piecewise constant time dependent Heston model

\begin{array}{rcl} d \ln S_t&=& \left(r_t - q_t - \frac{\nu_t}{2}\right)dt + \sqrt\nu_t dW^{S}_t \nonumber \\ d\nu_t&=& \kappa_t\left(\theta_t-\nu_t \right ) dt + \sigma_t\sqrt\nu_t dW^{\nu}_t \nonumber \\ \rho_t dt &=& dW^{S}_tdW^{\nu}_t \end{array}

with n time intervals [t_0=0, t_1], ... ,[t_{n-1}, t_n] and constant parameters within these intervals

\kappa_t = \kappa_j \wedge \theta_t = \theta_j \wedge \sigma_t=\sigma_j \wedge \rho_t=\rho_j \forall j\in [1, n] \wedge t\in [t_{j-1}, t_j]

is given by the recurrence relation

\begin{array}{rcl}  k_j &=& \kappa_j - iz\rho_j\sigma_j \nonumber \\ d_j &=& \sqrt{k_j^2 +\sigma_j^2(z^2+iz)} \nonumber \\  g_j &=&\displaystyle \frac{k_j- d_j}{k_j + d_j} \nonumber \\ \tilde{g_j} &=& \displaystyle \frac{k_j- d_j - D_{j+1}\sigma_j^2}{k_j + d_j - D_{j+1}\sigma_j^2} \nonumber \\ D_j &=&  \displaystyle \frac{k_j + d_j}{\sigma_j^2}\frac{g_j-\tilde{g_j} e^{-d_j \tau_j}}{1-\tilde{g_j} e^{-d_j\tau_j}} \nonumber \\ C_j &=& \displaystyle \frac{\kappa_j\theta_j}{\sigma_j^2} \left( (k_j- d_j )\tau_j - 2\ln\left(\frac{1-\tilde{g_j}e^{-d_j\tau_j}}{1-\tilde{g_j}}\right)\right) + C_{j+1} \nonumber \\ \tau_j &=& t_j - t_{j-1} \nonumber \\ \phi_{t_n}(z) &=& \exp{\left(C_1(z)+D_1(z)\nu_0\right)}\end{array}

and the initial condition

C_{n+1} = D_{n+1} = 0

It should be noted that the complex logarithm in this formulation can be restricted to the principal branch without introducing any discontinuities [1]. It is important to know the asymptotic behaviour of C_1(z)+D_1(z)\nu_0 in order to calculate the truncation point for the integral over the characteristic function when using the Andersen-Piterbarg approach with control variate [2]. A Mathematica script gives

\begin{array}{rcl} \displaystyle \lim_{u\to\infty}\frac{D_1(u)}{u} &=& \displaystyle-\frac{i\rho_1 + \sqrt{1-\rho_1^2} }{\sigma_1} \nonumber \\ \displaystyle\lim_{u\to\infty}\frac{C_1(u)}{u} &=& \displaystyle -\sum_{j=1}^n \frac{\kappa_j\theta_j}{\sigma_j}\left(\sqrt{1-\rho_j^2} + i\rho_j\right)\tau_j\end{array}

The implementation of the Andersen-Piterbarg method for the piecewise constant time dependent Heston model is part of the pull request #251.

[1] Afshani, S. (2010) Complex logarithms and the piecewise constant extension of the Heston model.

[2] Andersen, L. and Piterbarg, V. (2010) Interest Rate Modeling, Volume I: Foundations and Vanilla Models, (Atlantic Financial Press London).

New Semi-Analytic Heston Pricing Algorithms

hpcquantlibLeave a comment

The pricing engines for the Heston model in QuantLib have aged a bit. Meanwhile newer and better algorithms have been developed and discussed in the literature. For a comprehensive review please see [1][2]. Time to refurbish the existing engines.

The Heston model is defined by the following stochastic differential equation of the log spot

\begin{array}{rcl} d \ln S_t&=& \left(r_t - q_t - \frac{\nu_t}{2}\right)dt + \sqrt\nu_t dW^{S}_t \nonumber \\ d\nu_t&=& \kappa\left(\theta-\nu_t \right ) dt + \sigma\sqrt\nu_t dW^{\nu}_t \nonumber \\ \rho dt &=& dW^{S}_tdW^{\nu}_t \end{array}

The normalized characteristic function in the Gatheral formulation is given by

\phi_t(z) = \exp\left\{\frac{v_0}{\sigma^2}\frac{1-e^{-Dt}}{1-Ge^{-Dt}}\left(\kappa-i\rho\sigma z-D\right) + \frac{\kappa\theta}{\sigma^2}\left(\left(\kappa-i\rho\sigma z-D\right)t-2\ln\frac{1-Ge^{-Dt}}{1-G}\right) \right\}

\begin{array}{rcl} D&=&\sqrt{\left(\kappa - i\rho\sigma z\right)^2+\left(z^2+iz\right)\sigma^2} \nonumber \\ G&=&\displaystyle\frac{\kappa -i\rho\sigma z-D}{\kappa -i\rho\sigma z + D}\end{array}

Andersen and Piterbarg [3] introduced a Black-Scholes control variate to improve the numerical stability of Lewis’s formula (2001) for the price of a vanilla call option

\begin{array}{rcl} C(S_0, K, T)&=&BS\left(S_0, K, T, \sigma_{BS}\right) \nonumber \\ &+& \frac{\sqrt{Se^{(r-q)t}K}e^{-rt}}{\pi}\displaystyle\int\limits_{0}^{\infty}{\Re\left( \frac{\phi^{BS}_T\left(u-\frac{i}{2}\right) - \phi_T\left(u-\frac{i}{2}\right)}{u^2+\frac{1}{4}} e^{i u \left(\ln\frac{S}{K}+(r-q)T\right) } \right)  \mathrm{d}u}\end{array}

with the Black-Scholes price BS\left(S_0, K, T, \sigma_{BS}\right) of a vanilla call option with volatility \sigma_{BS} and the characteristic function

\phi^{BS}_T(z)=e^{-\frac{1}{2}\sigma_{BS}^2 T (z^2 + iz)}

Different proposals for the volatility of the vanilla option have been brought up in order to achieve an optimal control variate:

  • \sigma_{BS}^2 = \nu_0
  • \sigma_{BS}^2 = \frac{1}{\kappa t}\left(1-e^{-\kappa t}\right)\left(\nu_0-\theta\right) + \theta
  • \sigma_{BS}^2 = c_2 with c_2 the second cumulant of \ln \frac{F}{K} [1].

The diagram below shows the resulting integrand for the different control variate volatilities, the Heston model parameters

r=7.5\%, q=5\%, S_0=100, \nu_0=0.08, \rho=-80\%,\sigma=0.5, \kappa=4, \theta=0.05

and for a vanilla call option with strike K=200 and maturity t=2.0.

control_variate

The best choice depends on the Heston- and option parameters but it seems that

\sigma_{BS}^2 = \frac{1}{\kappa t}\left(1-e^{-\kappa t}\right)\left(\nu_0-\theta\right) + \theta

is a good option for a large variety of parameters. A direct comparison of the integrand with and without control variate shows how effective the control variate is.

control_comparison

The central trick of Andersen-Piterbarg is now to truncate Lewis’s infinite integral at a finite u_{max} such that the remaining part is smaller than a given threshold based on the following inequation

\displaystyle\int_{u_{max}}^\infty \left\vert \frac{\phi_B(u-\frac{i}{2}) - \phi(u-\frac{i}{2})}{u^2+\frac{1}{4}} \right\vert du \le e^{-C_\infty u_{max}}\displaystyle\int_{u_{max}}^\infty \frac{1}{u^2} du

Please see the original paper or [2] for further details on how to get from here to an algorithm for u_{max}, especially for short dated options.

As Andersen and Piterbarg have pointed out, the simple trapezoidal rule works surprisingly well to carry out the resulting integral if a control variate is used. Alan Lewis’s test cases with high precision Heston reference prices should serve as a test bed here, in particular the call with strike 100. The Gauss-Laguerre method is very thankful for this particular test case but the point here is that the trapezoidal rule converges much faster than the higher order Simpson rule or the even more complex adaptive Gauss-Lobatto method.

convergence_speed

The diagram below compares the recently added COS engine with the Andersen-Piterbarg method using the trapezoid rule for this test-case. The Andersen-Piterbarg method is more accurate for a similar number of points.

convergence_comparison

The implementation is part of the pull request #251.

[1] M. Schmelzle, Option Pricing Formulae using Fourier Transform: Theory and Application.

[2] F. Le Floc’h, Fourier Integration and Stochastic Volatility Calibration.

[3] L. Andersen, and V. Piterbarg, 2010,  Interest Rate Modeling, Volume I: Foundations and Vanilla Models,  Atlantic Financial Press London.

The CLV Model with a Square Root Kernel Process

hpcquantlibLeave a comment

The kernel process of the Collocating Local Volatility (CLV) model [1] defines the forward skew dynamics of the model. Although the Ornstein-Uhlenbeck process has some nice analytical features it offers sometimes too little control over the forward skew dynamics. The square root process

d\nu_t=\kappa(\theta-\nu_t) dt + \sigma \sqrt{\nu_t}dW_t

is a promising kernel to get more control over the forward skew. First, exact sampling is pretty easy for the square root process. The probability density function of \nu_t given \nu_0 is

\nu_t = \frac{\sigma^2 (1- e^{-\kappa t})}{4\kappa} \chi_d^{'2}\left(\frac{4\kappa e ^{-\kappa t}}{\sigma^2(1-e^{-\kappa t})}\nu_0\right)

where \chi_d^{'2} denotes the noncentral chi-squared distribution with

d=\frac{4\theta\kappa}{\sigma^2}

degrees of freedom and the noncentrality parameter

\lambda = \frac{4\kappa e^{-\kappa t}}{\sigma^2(1-e^{-\kappa t})}\nu_0

The boost library provides an efficient and accurate implementation of the inverse of the cumulative noncentral chi-squared distribution function, which can be used for an exact Monte-Carlo sampling scheme.

The optimal collocation points are given by the Gaussian quadrature points, which are defined by the zeros of the corresponding orthogonal polynomials. The orthogonal polynomials are defined by a recurrence relation and the collocation points are given by the eigenvalues of a symmetric, tri-diagonal matrix with the diagonal \{\alpha_i\} and the minor diagonal \{\sqrt{\beta_i}\}. Again these vectors are defined by a recurrence relation [2]

\begin{array}{rcl} z_{k, i} &=& z_{k-1, i+1} + \alpha_k z_{k-1,i} - \beta_k z_{k-2,i} \\ \nonumber \alpha_{k+1} &=& \frac{z_{k-1,k}}{k-1,k-1}-\frac{z_{k,k+1}}{z_{k,k}} \\ \nonumber \beta_{k+1} &=& \frac{z_{k,k}}{z_{k-1,k-1}}\\ \nonumber z_{-1,i} &=& 0 \\ \nonumber z_{0,i} &=& \mu_i = \int_0^\infty x^i\chi_{d,\lambda}^{'2}(x) dx \\ \nonumber \alpha_1 &=& -\frac{\mu_1}{\mu_0} \\ \nonumber \beta_1 &=& \mu_0 = 1\end{array}

The first n moments \mu_i(d, \lambda), i=0,..,n of the noncentral chi-squared distribution can be calculated using Mathematica and exported as plain C code in order to be integrated into QuantLib.

m[n_] := CForm[ Expectation[X^n, X \[Distributed]
NoncentralChiSquareDistribution[d, lambda]] // Simplify]

Solving the recurrence relation is often ill-conditioned when using double precision. Therefore the resulting equations will be solved using the Boost.Multiprecision package with the cpp_dec_float back-end (header only and dependency free). The results have been tested accordingly to the proposal in [3]. For any reasonable number of collocation points a precision of up to 100 digits seems to be sufficient. On the other side the computation even with 100 digits is really fast. The eigenvalue calculation is carried out using standard double precision.

The implementation of the calibration routine follows the same way as for the Ornstein-Uhlenbeck kernel process. Again the calibration is pretty fast and accurate compared with other structured models even though it involves multiple precision arithmetic.

Example: Forward volatility skew

  • Market prices are given by a Heston model with

S_0=100, r=0.1, q=0.05, \nu_0=0.09, \kappa=1.0, \theta=0.06, \sigma=0.4, \rho=-0.75.

  • SquareRoot-CLV kernel process parameters are given by

\kappa=0.2, \theta=0.09, \sigma=0.1, x_0=0.09

The diagram below shows the implied volatility of an forward starting European option with moneyness varying from 0.5 to 2 and maturity date six month after the reset date.

circlvforwardskew1

Same plot but with the following parameters of the square root kernel process

\kappa=0.1, \theta=0.09, \sigma=0.1, x_0=0.25

circlvforwardskew2

Source code is available here. The code is in an early state and needs some more testing/clean-up before a pull request can be made out of it. Next step will be to calibrate such a CIR-CLV model to the forward skew dynamics of an Heston Stochastic Local Volatility model.

[1] A. Grzelak, 2016, The CLV Framework – A Fresh Look at Efficient Pricing with Smile

[2] M. Morandi Cecchi and M. Redivo Zaglia, Computing the coefficients of a recurrence formula for numerical integration by moments and modified moments.

[3] Walter Gautschi, How and How not to check Gaussian Quadrature Formulae.

The Collocating Local Volatility Model

hpcquantlibLeave a comment

In a recent paper [1] Lech Grzelak has introduced his Collocating Local Volatility Model (CLV). This model utilises the so called collocation method [2] to map the cumulative distribution function of an arbitrary kernel process onto the true cumulative distribution function (CDF) extracted from option market quotes.

Starting point for the collocating local volatility model is the market implied CDF of an underlying S_t at time t_i:

F_{S(t_i)}(x) = 1 + e^{r t_i}\frac{\partial V_{call}(t_0, t_i, K)}{\partial K}\mid_{K=x} = e^{r t_i}\frac{\partial V_{put}(t_0, t_i, K)}{\partial K}\mid_{K=x}

The prices can also be given by another calibrated pricing model, e.g. the Heston model or the SABR model. To increase numerical stability it is best to use OTM calls and puts.

The dynamics of the spot process S_t should be given by some stochastic process X_t and a deterministic mapping function g(t, x) such that

S_t=g(t, X_t)

The mapping function g(t, x) ensures that the terminal distribution of S_t given by the CLV model matches the market implied CDF. The model then reads

\begin{array}{rcl} S_t &=& g(t,X_t)  \nonumber \\ dX_t &=& \mu(X_t)dt + \sigma(X_t)dW_t\nonumber \end{array}

The choice of the stochastic process X_t does not influence the model prices of vanilla European options – they are given by the market implied terminal CDF – but influences the dynamics of the forward volatility skew generated by the model and therefore the prices of more exotic options. It is also preferable to choose an analytical trackable process X_t to reduce the computational efforts.

The collocation methods outlined in [2] defines an efficient algorithm to approximate the mapping function g(t, x) based on a set of collocation points x_{i,j}=x_j(T_i) for a given set of maturities T_i, i=1,...,m and j=1,...,n interpolation points per maturity T_i. Let F_{S_{T_i}}(s) be the market implied CDF for a given expiry T_i. Then we get

\begin{array}{rcl}  F_{X_{T_i}}\left(x_{i,j}\right) &=& F_{S_{T_i}}\left(g(T_i, x_{i,j})\right) = F_{S_{T_i}}  \left(s_{i,j}\right) \nonumber \\ \Rightarrow s_{i,j}&=&F^{-1}_{S_{T_i}}\left(F_{X_{T_i}}(x_{i,j})\right) \nonumber \end{array}

for the collocation points with s_{i,j}=g(T_i, x_{i,j}).

The optimal collocation points are given by the abscissas of the Gaussian quadrature for the terminal distribution of X_{T_i}. The simplest choice is a normally distribute kernel process X_t like the Ornstein-Uhlenbeck process

dX_t = \kappa(\theta-X_t)dt + \sigma dW_t.

The corresponding collocation points of the Normal-CLV model are then given by

\begin{array}{rcl} x_j(t) &=& \mathbb{E}\left[X_t\right] + \sqrt{\mathbb{V}ar\left[X_t\right]} x_j^{\mathcal{N}(0,1)}  \nonumber \\ &=& \theta + \left(X_0 - \theta)e^{-\kappa t}\right) + \sqrt{\frac{\sigma^2}{2\kappa}\left(1-e^{-2\kappa t}\right)} x_j^{\mathcal{N}(0,1)}, \ j=1,...,n\end{array}

in which the collocation points x_j^{\mathcal{N}(0,1)} of the standard normal distribution can be calculated by QuantLib’s Gauss-Hermite quadrature implementation

Array abscissas = std::sqrt(2)*GaussHermiteIntegration(n).x()

Lagrange polynomials [3] are an efficient interpolation scheme to interpolate the mapping function g(t, x) between the collocation points

g(t, X_t) = \sum_{j=1}^N s_j (t)\prod_{k=1, j\neq k}^N \frac{X(t)-x_j(t)}{x_k(t)-x_j(t)}

Strictly speaking Lagrange polynomials do not preserve monotonicity and one could also use monotonic interpolation schemes supported by QuantLib’s spline interpolation routines. As outlined in [2] this method can also be used to approximate the inverse CDF of an “expensive” distributions.

Calibration of the Normal-CLV model to market prices is therefore pretty fast and straight forward as it takes the calibration of g(t, x_t).

Monte-Carlo pricing means simulating the trackable process X_t and evaluate the  Lagrange polynomial if the value of the spot process S_t is needed. Pricing via partial differential equation involves the one dimensinal PDE

\frac{\partial V}{\partial t} + \mu(x)\frac{\partial V}{\partial x} + \frac{1}{2}\sigma^2(x)\frac{\partial^2 V}{\partial x^2}-rV = 0

with the terminal condition at maturity time T

V(T, x_T) = \text{Payoff}\left(S_T=g(T,x_T)\right)

For plain vanilla options the upper and lower boundary condition is

\frac{\partial^2 V}{\partial x^2} = 0 \ \ \forall x\in\left\{x_{min},x_{max}\right\}

Example 1: Pricing error for plain vanilla options

  • Market prices are given by the Black-Scholes-Merton model with

S_0=100, r=0.1, q=0.04, \sigma=0.25.

  • Normal-CLV process parameters are given by

\kappa=1.0, \theta=0.1,\sigma=0.5,x_0=0.1

Ten collocation points are used to define the mapping function g(t, x) and the time to maturity is one year. The diagram below shows the deviation of the implied volatility of the Normal-CLV price based on the PDE solution from the true value of 25%

clvpriceerror

Even ten collocation points are already enough to obtain a very small pricing error. The advice in [2] to stretch the collocation grid has turned out to be very efficient if the number of collocation points gets larger.

Example 2: Forward volatility skew

  • Market prices are given by the Heston model with

S_0=100, r=0.1, q=0.05, \nu_o=0.09, \kappa=1.0, \theta=0.06, \sigma=0.4, \rho=-0.75.

  • Normal-CLV process parameters are given by

\kappa=-0.075, \theta=0.05,\sigma=0.25,x_0=0.05

The diagram below shows the implied volatility of an forward starting European option with moneyness varying from 0.5 to 2 and maturity date six month after the reset date.

hestonforward

The shape of the forward volatility surface of the Normal-CLV model shares important similarities with the surfaces of the more complex Heston or Heston Stochastic Local Volatility (Heston-SLV) model with large mixing angles \eta. But by the very nature of the Normal-CLV model, the forward volatility does not depend on the values of \theta, \sigma or x_0, which limits the variety of different forward skew dynamics this model can create. CLV models with non-normal kernel processes will support a greater variety.

Example 3: Pricing of Double-no-Touch options

  • Market prices are given by the Heston model with

S_0=100, r=0.02, q=0.01, \nu_o=0.09, \kappa=1.0, \theta=0.06, \sigma=0.8\rho=-0.8.

  • Normal-CLV process parameters are given by different \kappa values and

\theta=100,\sigma=0.15,x_0=100.0

Unsurprisingly the prices of 1Y Double-no-Touch options exhibit similar patterns with the Normal-CLV model and the Heston-SLV model as shown below in the “Moustache” graph. But the computational efforts using the Normal-CLV model are much smaller than the efforts for calibrating and solving the Heston-SLV model.

moustache.png

The QuantLib implementation of the Normal-CLV model is available as a pull request #117, the Rcpp based package Rclv contains the R interface to the QuantLib implementation and the demo code for all three examples.

[1] A. Grzelak, 2016, The CLV Framework – A Fresh Look at Efficient Pricing with Smile

[2] L.A. Grzelak, J.A.S. Witteveen, M.Suárez-Taboada, C.W. Oosterlee,
The Stochastic Collocation Monte Carlo Sampler: Highly efficient sampling from “expensive” distributions

[3] J-P. Berrut, L.N. Trefethen,  Barycentric Lagrange interpolation, SIAM Review, 46(3):501–517, 2004.