Accelerate VPP Pricing on the GPU using Preconditioner

As outlined in Solving Sparse Linear Systems using CUSP and CUDA a VPP pricing problem can be solved on a GPU using finite difference methods based on the implicit Euler scheme. The matrix inversion is carried out on the GPU using the BiCGStab solver of the numerical library cusp [1]. The  optimal exercise problem for the VPP can be solved using dynamic programming for every time slice on the GPU (see [2]).

The lower convergence order in \Delta t of the implicit Euler scheme is not a disadvantage because the hourly optimization step forces \Delta t = 1h anyhow. This step size is small enough for the implicit Euler scheme to obtain a very good accuracy. The explicit Euler scheme is  numerically unstable, even if \Delta t = 1h. It converged only for large \Delta x, \Delta y and \Delta u.

Clearly the matrix inversion dominates the overall runtime of the pricing routine. Preconditioning is a standard technique to improve the performance of the involved BiCGStab algorithm. The cusp library supports a set of suited preconditioner for this problem.

  • Diagonal: preconditioning is very inexpensive to use, but has often limited effectiveness.
  • Approximate Inverse: Non-symmetric Bridson with different dropping strategies.
  • Smoothed Aggregation: Algebraic Multi-grid preconditioner based on smoothed aggregation.

On a CPU QuantLib supports preconditioning based on the Thomas algorithm. This algorithm solves a tridiagonal matrix with O(8n) operations. The corresponding preconditioner can be seen as a natural extension of the diagonal preconditioner. Experiments on a CPU indicate that  the tridiagonal preconditioner leads to a significant speed-up for our VPP pricing problem. The upcoming CUSPARSE library of the CUDA 4.1 release contains a tridiagonal solver for the GPU. With the help of a small interface class the CUDA tridiagonal solver can be used as a cusp preconditioner. The source code is available here. The code depends on the latest QuantLib version from the SVN trunk.

The testbed is a VPP contract with t_{minUp} = t_{minDown} = 2h and a maturity of 4 weeks (corresponds to 4*7*24=672 exercise steps). Therefore in addition to the three dimensions for the stochastic model the problem has a fourth dimension of size 2t_{minUp} + t_{minDown} = 6  to solve the optimization problem via dynamic programming.

The following diagram shows the different run-times of the algorithms on a GPU and a CPU.

The fastest algorithm on the GPU (implicit Euler scheme with BiCGStab and tridiagonal preconditioner) outperforms the best CPU implementation (operator splitting based on the Douglas scheme) roughly by a factor of ten. The non-symmetric Bridson preconditioner with the dropping strategy from Lin & More on the GPU is the second best algorithm but it ran out of memory on larger grids. The diagonal preconditioner slows down the overall performance on the GPU and the BiCGStab algorithm with the smoothed aggregation preconditioner did not converge.

[1] cusp, Generic Parallel Algorithms for Sparse Matrix and Graph Computations

[2] Pricing of a Virtual Power Plant on a GPU

Advertisements

VPP Pricing IV: Variance Reduction for Perfect Foresight

Even though perfect foresight provides only an upper bound to the real VPP value the differences are often neglectable and the implementation efforts  are small compared with “exact” pricing based on finite difference methods or least square Monte-Carlo. Perfect foresight is the method of choice in conjunction with a linear programming optimizer if the problem contains time-integral constraints. Therefore it is worth to test the efficiency of two standard variance reduction techniques, namely antithetic sampling and quasi Monte-Carlo (QMC) together with a Brownian Bridge. Both methods are explained in [1], antithetic sampling in chapter 4.2 and quasi Monte-Carlo in section 5. Randomized QMC is used to calculate the error estimates for QMC as it is outlined in chapter 5.4.

Using the parameterization of the previous section VPP Pricing III, QMC in conjunction with a Brownian Bridge clearly out-performance the other algorithms for a 6 month contract as can be seen in the diagrams below. The code is available here. It depends on the latest QuantLib version from the SVN trunk or the upcoming QuantLib 1.2 release. If you want to generate the plots you’ll also need R.

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

VPP Pricing III: Exact Pricing based on Finite Difference Methods

The total value of a virtual power pant (VPP) can be decomposed in an intrinsic part plus an extrinsic part. The intrinsic value is given by the cash-flows that the VPP would generate based on the current power and gas forward curve. Therefore the intrinsic value can be calculated without defining a stochastic model for the power and gas prices using either linear or dynamic optimization methods. Calculating the extrinsic value implies pricing the VPP exactly to calculate the total value for a given stochastic process. The model in use here is outlined in the article VPP Pricing I: Stochastic Processes & Partial Integro Differential Equation. Exact pricing can  be done using least square Monte-Carlo or finite difference methods using dynamic programming for the local optimization [1].

The focus is on the finite difference method, which involves solving the three-dimensional partial integro differential equation

\begin{array}{rcl} rV&=&\frac{\partial V}{\partial t}+\frac{\sigma_x^2}{2}\frac{\partial^2 V}{\partial x^2}-\alpha x\frac{\partial V}{\partial x}-\beta y \frac{\partial V}{\partial y} \\[6pt]&+&\frac{\sigma_u^2}{2}\frac{\partial^2 V}{\partial u^2}- \kappa u\frac{\partial V}{\partial u} +\rho\sigma_x\sigma_u\frac{\partial^2 V}{\partial x\partial u}\\[6pt] &+&\lambda\int_\mathbb{R}\left(V(x,y+z,u,t)-V(x,y,u,t) \right )\omega(z)dz \\ \end{array}

An additional fourth dimension is needed to keep track of the different states of the VPP. Based on the characteristics of the VPP type outlined in the article VPP Pricing II: Mixed Integer Linear Programming a VPP with three possible load levels P \in \{0, P_{min}, P_{max} \} has 2 t_{up} + t_{down} different states.

Pricing via Monte-Carlo and perfect foresight involves simulating the stochastic process

\begin{array}{rcl} P_t &=& \exp(p_t + X_t + Y_t) \\ dX_t &=& -\alpha X_tdt + \sigma_x dW_t^x \\ dY_t &=& -\beta Y_{t-}dt+J_tdN_t \\ \omega(J) &=& \eta e^{-\eta J} \\ G_t &=& \exp(g_t + U_t) \\ dU_t &=& -\kappa U_tdt+\sigma_udW_t^u \\ \rho &=& \mathrm{corr} (dW_t^x, dW_t^u) \end{array}

and optimize the power plant load schedule for each path separately in the same way the intrinsic value is calculated. This procedure will result in an upper bound for the exact price of the VPP. The prices for a 4 weeks VPP contract based on this two methods and with the parameters outlined in VPP Pricing II show almost no differences between the “exact” finite difference method and the perfect foresight  upper bound value beside the Monte-Carlo error (s. diagram below, compare with [1]).

The size of the lattice for the finite difference method was

(t, nPower, nJump, nGas, nStates) = (672, 101, 51, 20, 14).

The code is available here. It depends on the latest QuantLib version from the SVN trunk or the upcoming QuantLib 1.2 release. If you want to generate the plot you’ll also need R

[1] H. van Dijken,  D. van Abbena, H.S. Los, C. de Jong, The value of starting up the power plant.

VPP Pricing II: Mixed Integer Linear Programming

The next two steps are defining a simple VPP contract (or a simplified gas-run power plant) and setting up a mixed integer linear programming optimization (MIP) to calculated the intrinsic value and an upper bound for the extrinsic value based on a Monte-Carlo simulation and assuming perfect foresight. The third step outlined in the next part will then be the “exact” pricing of the extrinsic value using dynamic programming and finite difference methods.

The set-up of the simplified gas-run power plant is similar to the one explained in chapter 4.2.3 of the text-book [1]. In general the power plant has three power output level:

  • Plant is off, P=0
  • Generation at minimum load P=P_{min}
  • Generation at maximum load P=P_{max}

The power plant has a fixed efficiency rate

\zeta=\frac{MWh Power Output}{MWh Heat}.

Ramp rates will be neglected, but the power plant has a minimum uptime t_{up} and a minimum downtime t_{down}. The start-up costs are given by a fixed start-up cost \eta (in €) and the price of the gas needed to produce the start-up heat \theta (in MWh).

The mixed integer linear optimization is running in one hour blocks and is using three decision variables per hour i. The binary decision variable \beta_i is true if the power plant is running at minimum load P_{min} or at maximum load P_{max} and \beta_i is false if the plant is off. The real decision variable 0\le s_i \le 1 is equal to one if the plant is started in hour i, which is implied by the following constraint

\beta_i - \beta_{i-1} \le s_i.

The minimum up-time t_{up} and the minimum down-time t_{down} is a consequence of the constraints

The real decision variable 0\le \gamma \le 1 is equal to one if the power plant is running at maximum load P_{max} and zero if the power plant is either running at minimum load P_{min} or if the plant is off, that means

\beta_i \ge \gamma_i

Let P_i be the power price, G_i be the gas price and CO_2^i be the carbon dioxide price at hour i. The objective function is then given by

For a one year span the problem consists of 3\cdot 365 \cdot 24 = 26280 decision variables \{\beta_i, \gamma_i, s_i\} and 4 \cdot 365\cdot 24 = 35040 constraints. This comparable small problem can be solved using e.g. the Gnu Linear Programming Kit (GLPK). For an overview on open source linear/mixed integer programming solver see [2].

The model parameters and the example forward curves are outlined in the previous entry VPP Pricing I. The diagram below shows the intrinsic value and the upper bound for the total value (intrinsic plus extrinsic value) based on Monte-Carlo, perfect foresight and MIP for different  power plant efficiencies \zeta. The parameters of the VPP contract are given by

t_{up}=t_{down}=2h, P_{min}=8MW, P_{max}=40MW, \eta=300 EUR, \theta=20MWh,

the (fixed) carbon dioxide price is 3.0€ per MWh heat.

The source code is available here. It depends on GLPK and the latest QuantLib version from the SVN trunk or the next QuantLib 1.2 release.

It is now quite easy to add and price time-integral constraints, e.g. the following constraint restricts the number of starts within a year to be less than or equal to a given number

\sum_{t=1}^N s_i \le \#Starts.

The following diagram shows the results for \#Starts \le 25 and a minimum load P_{min}=25MW.

The source code is available here. It depends on QuantLib 1.1 and if  you want to generate the plot directly from the C++ program you’ll also need R, RCPP and RInside.

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

[2] S. R. Thorncraft, Evaluation of Open-Source LP Optimization Codes in Solving Electricity Spot Market Optimization Problems.

VPP Pricing I: Stochastic Processes & Partial Integro Differential Equation

The basis for the virtual power plant (VPP) pricing algorithm is a joint stochastic process for the power price and the gas price. Instead of defining a  stochastic differential equation for the spark spread directly the spark spread is then given by the difference between the power price and the gas price times the VPP heating rate.

The Kluge model will be used to define the power price process [1], an exponential Ornstein-Uhlenbeck will describe the gas price [2]

where N_t is a Poisson process with jump intensity \lambda. To match the power forward curve F_0^t the seasonal function p_t is given by [1]:

To be consistent with the gas forward curve H_0^t the seasonal function g_t is defined by

The Feynman-Kac theorem can be applied to derive the corresponding three-dimensional  partial integro differential equation

In general at least one further dimension is needed to keep track of the state of the virtual power plant. Therefore solving this model using finite difference methods will lead to a four-dimensional PIDE problem.

The following diagram shows an one year example path for the power and the gas price based on the freely available Kyos example forward curves. The sample parameterization is affected by [1] for the power process and [2] for the gas process.

\beta=200, \eta=2.5, \lambda=4, \alpha=7, \sigma_x=1.4, \kappa=4.45, \sigma_u=\sqrt{1.3}, \rho=70\%

The code is available here. It depends 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] T. Kluge, Pricing Swing Options and other Electricity Derivatives

[2] G. Fusai, A. Roncoroni, Implementing Models in Quantitative Finance: Models and Cases, Chapter 19, ISBN: 978-3-540-22348-1

Swing Option: Linear vs. Dynamic Programming II

In order to get an improved impression on the differences between Monte-Carlo/linear programming and  finite difference methods/dynamic programming  a more realistic swing option with hourly payoff profile is evaluated based on a German hourly forward curve. The forward curve is taken from the Kyos example download page. The parameterization of the Kluge model is outlined in [1].

With maturity of twelve weeks the example swing option provides 12*7*24=2016 exercise opportunities. The number of overall exercised swing opportunities is constraint by

100 \le \sum_{i=1}^{2016} \beta_i \le 500.

The size of two dimensions of the finite difference method (dynamic programming) is therefore already been given. 2016 steps are needed in time direction and 500 steps are needed in the “consumed exercises” direction. Together with the two other directions – one for the power price and one dimension for the jump process – and without further simplidsfications this forms a pretty large finite difference problem. The Monte-Carlo based linear programming approach reduces the computational burden but will lead to an upper bound of the swing option price. The diagram below shows the corresponding results.



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. To utilize all CPU cores please use the -DNTHREADS=(number of CPU cores)  compiler switch. In addition to run the program you have to download the forward curve German power from the Kyos web page and convert it into a text file EEX_2010-2013.txt of the folliwing format

4-Oct-2010 36.73 32.09 27.32 23.22 25.71 35.60 47.32 …

5-Oct-2010 38.12 31.12 22.76 25.65 27.87 34.60 50.01 …

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

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