Swing Option: Linear vs. Dynamic Programming I

May 29, 2011February 18, 2014 hpcquantlib4 Comments

Finding the optimal exercise strategy for a swing option is a challenging task, even for simple payoff structures. In the absence of complicated time-integral constraints dynamic programming can be used to calculate the optimal exercise strategy. Linear programming applied ex post over the whole spot price path is capable to deal with more complicated time-integral constraints but the algorithm leads only to an upper bound of the real option price [1].

A basic payoff structure is e.g.: On the exercise dates $t_i, i\in\{1..N\}$ the decision variables $\beta_i\in\{0,1\}$ indicate whether to exercise $(\beta_i=1)$ or not to exercise $(\beta_i=0)$ the option at a given point $t_i$ in time. Payoff is either $\beta_i (P(t_i)-K)$ for a call or $\beta_i(K-P(t_i))$ for a put option. Here $K$ is the strike price and $P(t_i)$ is the spot price . The number of exercise opportunities is constraint by

$E_{min} \le \sum_{i=1}^N \beta_i \le E_{max}$ .

The price of a swing call option is then given by

$\text{npv} = \max_{\beta}\mathbb{E}^\mathbb{Q} \left[ \sum_{i=1}^N \beta_i \left( P(t_i) - K\right)e^{rt_i}\right]$

Let the dynamics of the spot price in the risk neutral measure $\mathbb{Q}$ be given by the Kluge model [2]:

$\begin{array}{rcl} S &=& exp(X_t + Y_t) \\ dX_t &=& \alpha(\mu(t)-X_t)dt + \sigma dW_t \\ dY_t &=& -\beta Y_{t-}dt + J_tdN_t \\ \omega(J)&=& \eta_u e^{-\eta_u J} \end{array}$

This model is composed of an Ornstein-Uhlenbeck (OU) process $X_t$ , a deterministic period function $\mu(t)$ characterizing the seasonality and a mean reverting process $Y_t$ with jumps to incorporate spikes. $N_t$ is a Poisson process with jump intensity $\lambda$ and the jump size itself is exponentially distributed. The random variables $W_t, N_t$ and $J_t$ are independent. The Monte-Carlo path simulation will be built on top of standard components for an OU process and a Poisson jump-diffusion process [3]. The dynamic programming approach can either be carried out using least squares Monte-Carlo ansatz (Longstaff Schwartz algorithm) or using finite difference methods. We prefer finite difference methods to avoid e.g. the problem of choosing an appropriate basis function set. The Feynman-Kac theorem leads to the corresponding partial integro differential equation (PIDE):

$\begin{array}{rcl} rV =&&\frac{\partial V}{\partial t}+\frac{\sigma^2}{2}\frac{\partial^2 V}{\partial x^2} + \alpha(\mu(t)-x)\frac{\partial V}{\partial x}-\beta y\frac{\partial V}{\partial y}\\&+&\lambda \int_\mathbb{R} \left( V(x,y+z,t)-V(x,y,t)\right )\omega(z) dz \end{array}$

A Gauss-Laguerre quadrature is appropriate to calculate the integral part of the PIDE. Beside this two-dimensional PIDE an additional dimension is needed to keep track of the already consumed exercise rights. The three-dimensional formulation together with Bellman’s principle of optimality transforms the global optimization problem into a local optimization problem.

Target of the linear programming approach is the upper bound

$\text{npv} \le \mathbb{E}^\mathbb{Q}\left[ \max_{\beta} \sum_{i=1}^N \beta_i\left(P(t_i) - K \right) e^{rt_i}\right]$

The linear programming algorithm will calculate the optimal exercise strategy $\left\{ \beta \right\}$ on each Monte-Carlo path separately (perfect foresight) with respect to the given constraints. For this basic type of swing option linear programming is sort of over-engineering because a simple sort algorithm will also reproduce the optimal exercise strategy. But introducing linear programming and mixed integer programming now will enable us later on to deal with more complicated time-integral constraints. Libraries for this task are freely available, e.g. the GNU Linear Programming Kit.

The test parameterization of the model is

$\alpha=1.0, \sigma=200\%, \mu(t)=3.0,\beta=5.0, \eta=2.0, \lambda=1.0, X_0=3.0, Y_0=0.0$ ,

the example swing call option has maturity of one year, strike price is equal 40, one exercise opportunity per month and the minimum number of exercise rights is equal to the maximum number of exercise rights, $E_{min} = E_{max}$ . Interest rates are at $r=10\%$ .

As the diagram above shows for this set-up the upper bound for the swing option price calculated using linear programming differs significantly from the correct price calculated based on dynamic programming on a lattice. For twelve exercise rights both algorithms have to provide the same results because the constraint forces to always exercise on each exercise date. Next to come is to rerun the simulation with a more realistic forward curve and process parameters.

The code is available here. It depends on the GNU Linear Programming Kit, the Boost Thread library for parallelization and at the time of writing on the latest QuantLib version from the SVN trunk. If you want to generate the plot directly out of the C++ program you also need R, RCPP and RInside.

[1] M. Burger, B. Graeber, G. Schindlmayr, Managing Energy Risk, ISDN 978-0-470-ß2962-6

[2] T. Kluge, Pricing Swing Options and other Electricity Derivatives

[3] P. Glasserman, Monte Carlo Methods in Financial Engineering. ISBN-0387004513

Calibrate the Heston Model on a GPU: CUDA 3.x and 4.0 Update

May 19, 2011May 29, 2011 hpcquantlibLeave a comment

Please find here a CUDA based program (upgraded for CUDA 3.x and 4.0) to price and calibrate a piecewise constant time dependent Heston model based on the QuantLib 1.0. The example code is in test-suite/hestonmodel.cpp. Some implementation details are outlined in the following paper, A. Bernemann, R. Schreyer, K. Spanderen: Accelerating Exotic Option Pricing and Model Calibration Using GPUs.

Calibration: Comparison of deterministic Optimizers

May 17, 2011May 29, 2011 hpcquantlibLeave a comment

The gaol is to calibrate a time dependent Heston model defined by the following SDE

$\begin{array}{rcl} dS(t, S) &=& \mu_t S dt + \sqrt{v} S dW_1 \\ dv(t, S) &=& \kappa_t (\theta_t - v) dt + \sigma_t \sqrt{v} dW_2 \\ dW_1 dW_2 &=& \rho_t dt \end{array}$

The parameter set $\{\kappa_t, \theta_t, ,\sigma_t, \rho_t\}$ is supposed to be piecewise constant in time. This model has a semi-closed solution for plain vanilla European put/call options based on the characteristic function method [1].

The DAX implied volatility surface based on July 5, 2002 data taken from [2] defines the “benchmark” calibration problem. The benchmark model parameters for the optimization problem are given by

$\Theta = \{\nu_0, \kappa_{0\leq t < 0.25}, \kappa_{0.25 \leq t},\theta, \sigma, \rho\}$ .

The non-linear least square optimization problem is defined by the goodness of fit measure

$\zeta=\min_{\Theta} \sum_{i=1}^N \left( \sigma_i^{market}(K,T) - \sigma_i^{model}(K,T, \Theta) \right)^2$

where $\sigma_i^{market}(K,T)$ is the market implied volatility for strike K and maturity T and $\sigma_i^{model}(K,T, \Theta)$ is the corresponding Black-Scholes volatility implied from the model price. The optimal solution for this problem is

$\Theta_{min} = \{ 0.2231, 39.651, 7.546, 0.0954, 5.1865, -0.5004 \}$

leading to a goodness of fit measure of $\zeta_{min}=74.4731$ (Please keep in mind that this result is the outcome of a naive calibration procedure. Due to the large $\kappa$ and $\sigma$ values I’d not use these parameters to price a derivative.).

The diagram above shows the “goodness of fit”-surface for the parameter sets in

$\Theta = \{ 0.2231, 39.651, \kappa_{0.25 \leq t} \in [0, 16], 0.0954, \sigma \in [0, 16], -0.5004 \}$

To be able to compare a larger number of deterministic optimizers the model calibration will be carried out using R and with help of the additional packages minpack.lm and minqa.

Non-linear Least Square Optimization:

nls.lm: Levenberg-Marquardt algorithm(based on MINPACK, also available in QuantLib)
nls: Gauss-Newton algorithm
nl2sol:based on the PORT library.

Non-linear (trusted region) Minimization:

nlm:Newton style minizer
bfgs: quasi-Newton method published by Broyden, Fletcher, Goldfarb and Shanno
l-bfgs-b: limited memory BFGS algorithm incl.box constraints
cg: conjugate gradient algorithm
nm: Nelder-Mead method
bobyqa:trust region method that forms quadratic models by interpolation.
newuoa: trust region method that forms quadratic models by interpolation.
uobyqa:trust region method that forms quadratic models by interpolation.

The particular result dependents on the starting vector but the following diagram shows a common outcome.

The best methods are the non-linear least square algorithms nls.lm and nls followed by nl2sol. Algorithms not included in the diagram have performed even worse than “nm” for this problem.

The goodness of fit measure is calculated in C++ based on the QuantLib and exposed to R using RCPP. The C++ code and the R scripts to perform the optimizations and to create the plots can be found here.

1] A. Elices, Models with time-dependent parameters using transform methods: application to Heston’s model,

[2] A. Sepp, Pricing European-Style Options under Jump Diffusion Processes with Stochastic Volatility: Applications of Fourier Transform.

Scrambled Sobol Sequences

May 11, 2011June 4, 2011 hpcquantlib2 Comments

The authors of [1][2] presented an efficient way to create scrambled Sobol numbers and have shown that their scrambling algorithm can lead to a significantly higher convergence speed. The upcoming CUDA 4.0 Release provides an implementation of the scrambling algorithm outlined in [2].

The aim is now to compare the order of convergence of a PRNG (Merssenne-Twister), a QRNG (Sobol) and this scrambled-QRNG for the pricing of a portfolio of geometric Asian options. To number of observation dates of the Asian options in the benchmark portfolio varies between 45 and 120 observation dates. The dimensionality of the pricing problem is considerably large. We’ll use a Brownian bridge to construct the paths. This ensures that the first coordinates of a Sobol sequence are used to sample the most important modes.

The error analysis is not an easy task. Simply applying another randomized shift is not sufficient as this transformation changes the discrepancy of a point-set (A shifted (t,m,d)-net does not need to be a (t,m,d)-net.). In [3] this problem was overcome by looking on the RMS relative error

$\sqrt{\frac{1}{m}\sum_{i=1}^m\left(\frac{\hat{C}_i -C_i}{C_i}\right)^2}$

of a set of benchmark options with true prices given by $C_1,...,C_m$ and the Monte-Carlo approximations $\hat{C}_1,...,\hat{C}_m$ .

We are in particular interest in ratios between the PRNG, QRNG and scrambled QRNG errors and in the order $\alpha$ of the convergence. For PRNG we are expecting a convergence of order $err \simeq N^{\alpha}$ with $\alpha=-0.5$ and $N$ being the number of paths.

As can be seen in the diagram the scrambled Sobol sequence provides the smallest overall error bounds but the convergence order between scrambled and non-scrambled Sobol does not differ significantly. As the implementation comes more or less for free it’s definitely worth to check scrambled Sobol sequences out.

C++ and CUDA code is available here. It depends on QuantLib 1.1 and CUDA 4.0. If you want to generate the plot directly out of the C++ program you also need R, RCPP and RInside.

[1] J. Dick, Higher order scrambled digital nets achieve the optimal rate of the root mean square error for smooth integrands. ArXiv, 8 July 2010.

[2] A.B. Owen, Local antithetic sampling with scrambled nets. ArXiv, 28 May 2008

[3] P. Glasserman, Monte Carlo Methods in Financial Engineering. ISBN-0387004513

HPC-QuantLib

Month: May 2011

Swing Option: Linear vs. Dynamic Programming I

Calibrate the Heston Model on a GPU: CUDA 3.x and 4.0 Update

Calibration: Comparison of deterministic Optimizers

Scrambled Sobol Sequences