Accelerate VPP Pricing on the GPU using Preconditioner

January 22, 2012January 30, 2012 hpcquantlibLeave a comment

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

VPP Pricing IV: Variance Reduction for Perfect Foresight

August 12, 2011May 4, 2013 hpcquantlibLeave a comment

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

August 7, 2011February 18, 2014 hpcquantlib3 Comments

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

June 23, 2011June 12, 2019 hpcquantlibLeave a comment

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

$\begin{array}{rcl} \beta_i &\ge& \sum_{t=i-t_{up}+1}^{t=i} s_t \\[7pt] \beta_i &\le& 1-\sum_{t=i+1}^{t=i+t_{down}} s_t \end{array}$

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

$P\& L = \sum_{t=1}^N\left[\left(\gamma_iP_{max} + P_{min}(\beta_i-\gamma_i)\right) \left(P_i - \frac{G_i+CO_2^i}{\zeta}\right) - s_i\left(\eta + \theta (G_i+CO_2^i)\right)\right]$

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

June 13, 2011October 17, 2019 hpcquantlib5 Comments

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]:

$p_t = \ln F_0^t -X_0 e^{-\alpha t}-Y_o e^{-\beta t} -\frac{\sigma_x^2}{4\alpha}\left(1-e^{-2\alpha t} \right ) - \frac{\lambda}{\beta}\ln\left( \frac{\eta-e^{-\beta t}}{\eta-1}\right), \eta\ge 1$

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

$g_t = \ln H_0^t -U_0 e^{-\kappa t}-\frac{\sigma_u^2}{4\kappa}\left(1-e^{-2\kappa t} \right )$

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

June 4, 2011May 4, 2013 hpcquantlib4 Comments

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

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.

[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

HPC-QuantLib

Category: Commodity