Swing Option: Linear vs. Dynamic Programming I

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

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

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

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:

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

Non-linear (trusted region) Minimization:

  1. nlm:Newton style minizer
  2. bfgs: quasi-Newton method published by Broyden, Fletcher, Goldfarb and Shanno
  3. l-bfgs-b: limited memory BFGS algorithm incl.box constraints
  4. cg: conjugate gradient algorithm
  5. nm: Nelder-Mead method
  6. bobyqa:trust region method that forms quadratic models by interpolation.
  7. newuoa: trust region method that forms quadratic models by interpolation.
  8. 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

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