How to speed up the quantization tree algorithm with an application to swing options `s Anne Laure Bronstein ∗, Gilles Page † and Benedikt Wilbertz‡ November 23, 2009 hal-00273790, version 2 - 15 Dec 2009 Abstract In this paper, we suggest several improvements to the numerical implementation of the quantization method for stochastic control problems in order to get fast and accurate premium estimations. This technique is applied to derivative pricing in energy markets. Several ways of modeling energy derivatives are described and finally numerical examples including parallel execution on multi-processor devices are presented to illustrate the accuracy of these methods and their execution times. Keywords: Quantization, Gaussian process, L´evy process, Stochastic control, Swing options, Backward Dynamic programming. MSC: 60G15, 60E99. 1 Introduction In the last decades, optimal quantization arose as an efficient method for pricing complex financial contracts. Based on these quantizations, it is possible to construct a tree method, the quantization tree algorithm, which has successfully been applied to the pricing of derivatives like American multi-asset options (see [1]) or energy contracts like swing options (see [3] and [2]). In financial institutions, quickness of execution as well as high accuracy are important criteria in the choice of a pricing method. With this observation in mind, we suggest some improvements to the original quantization tree method. We mainly focus on the fast computations of the transition probabilities for the quantization method. In higher dimension, all procedures devoted to the computation of transition probabilities are based on Monte Carlo simulations. They all share the same core, namely repeated nearest neighbor searches. So we implemented a classical fast nearest neighbor procedure: the kd-tree algorithm. However, the first aim of this paper is to present and test several alternative approaches to transition probability estimation by simulation that can be massively parallelized. By “massively” we mean that the procedures can be parallelized beyond the “standard” parallelization of the plain-vanilla Monte Carlo method (consisting in generating scenarii of the underlying structure process “through” the quantization tree). These new approaches allow a bigger number of threads to run completely independent, which helps to avoid time consuming tasks like thread ∗ Laboratoire de Probabilit´es et Mod`eles al´eatoires, UMR 7599, Universit´e Paris 6, case 188, 4, pl. Jussieu, F-75252 Paris Cedex 5. E-mail: [email protected] † Laboratoire de Probabilit´es et Mod`eles al´eatoires, UMR 7599, Universit´e Paris 6, case 188, 4, pl. Jussieu, F-75252 Paris Cedex 5. E-mail: [email protected] ‡ Universit¨ at Trier, FB IV-Mathematik, D-54286 Trier, Germany. E-mail: [email protected] 1 hal-00273790, version 2 - 15 Dec 2009 synchronization and communication. Some extensive numerical experiments have been carried out to compare the performances of these different approaches in term of accuracy and execution time (on both CPU and GPU). As concerns the numerical aspects, we focused on swing options pricing, which are commonly dealt on energy markets. The underlying assets for this family of derivative products usually are forward contracts on index on gas, oil or electricity. We considered several more or less classical dynamics for the underlying asset: First some Gaussian 1-factor model, where the trinomial tree method from [14] served as reference procedure, and later on a Gaussian 2-factor model with up to 365 exercise dates. But we also considered time discretized jump dynamics, which are often used to model the spot electricity price, namely the exponential of a Normal Inverse Gaussian (NIG) L´evy process. The NIG distribution has been introduced in finance by Barndorff-Nielsen in the 90s and has recently been applied to energy markets by several authors, e.g. Benth, Frestad and Koekebakker and Benth and Saltyte-Benth (see [5] and [6]). They suggest that the NIG family fits empirical electricity return distribution very well and represents an attractive alternative to the family of normal distributions. It is the first time to our knowledge that such a jump model is implemented in a quantization framework. Doing so, we wish to strongly emphasize the fact, that quantization tree methods are efficient to produce spatial discretization of any (simulatable) Feller Markov chain. Indeed, any (usual) time discretization scheme of a diffusion driven by a L´evy process (or any L´evy process observed at some discrete times) is a Feller Markov chain. The third contribution of this paper is to introduce a Richardson-Romberg extrapolation based on two quantization trees of different sizes. This principle has already been successfully tested for “linear” European option pricing (see [17] or [18]), but it is the first time it is introduced in a “non linear” setting as swing option pricing and stochastic control. This idea, which is based on several heuristic considerations, seems quite promising from a numerical point of view (see section 6) and therefore strongly suggests that the theoretical expansion of the error as a function of the quantization tree size holds true. The paper is organized as follows. In Section 2, some short background on optimal quantization is provided. In Section 3, we recall the quantization tree algorithm associated to the swing option stochastic control problem (see [3] and [2]). Moreover, a special case of swing option, the Call Strip, is developed. In Section 4, some suggestions are addressed to improve the execution time as well as the accuracy of the algorithm. Section 5 is devoted to the description of the treated dynamics and in Section 6 we finally present numerical results on both a single core and multi-core CPU resp. GPU setting. 2 Optimal quantization Optimal quantization has been developed in the 1950s in the field of Signal Processing. Its main purpose consists in approximating a continuous signal by a discrete one in an optimal way. In the 1990s, its application has been used in the field of Numerical Integration to derive some cubature formulae (see [16]). Later in the early 2000s, this method has been extended to the computation of conditional expectations and applied to the field of Numerical Probability and Financial Mathematics. This extension has been motivated by the necessity of designing efficient methodologies for pricing and hedging more and more sophisticated financial products, especially multi-asset American options (see [1]). Let (Ω, A, P) be a given probability space and let X be a random vector defined on this probability space and taking valued in Rd . Let N be a positive integer, the first step of quantization ˆ taking at most N values in a grid consists in discretizing X by a σ(X)-measurable random vector X Γ = {x1 , . . . , xN } ⊂ Rd . This grid or codebook is called an N -quantizer of X. Let x = (x1 , . . . , xN ) 2 ˆ a Borel function qx : Rd → Γ called denotes the N -tuple induced by Γ. One can associate with X ˆ = qx (X). quantizer such that X Let p ∈ [1, ∞[ and let X ∈ Lp (Ω, A, P). Optimal quantization consists in studying the best Lp approximation of X by qx (X) when x runs over (Rd )N : it amounts to minimize in q the Lp -mean quantization error, i.e. solving the minimization problem inf{|| X − q(X)||p | q : Rd → Rd Borel and |q(X(Ω))| ≤ N }. (1) Noting that |X − qx (X)| ≥ dist(X, Γ) = min |X − xi |, 1≤i≤N one may restrict to a Voronoi quantization as defined below. Definition 1. Let x = (x1 , . . . , xN ) ∈ (Rd )N . A partition (Ci (x))1≤i≤N of Rd is a Voronoi tessellation of the N -quantizer x, if for every i ∈ {1, . . . , N }, Ci (x) is a Borel set satisfying Ci (x) ⊂ {u ∈ Rd | | u − xi | = min | u − xj |}, hal-00273790, version 2 - 15 Dec 2009 1≤j≤N (2) where | . | is the Euclidean norm on Rd . The nearest neighbor projection πx on x induced by a Voronoi partition (Ci (x))1≤i≤N is defined for every u ∈ Rd by πx (u) := N X xi 1Ci (x) (u). (3) i=1 It maps the random vector X into ˆx = X N X xi 1Ci (x) (X) = πx (X), (4) i=1 which is called a Voronoi quantization of X. Now, the minimization problem (1) turns into an optimization problem on x ∈ (Rd )N . This second step depends on p and the probability distribution PX of X. It always has at least one solution x∗ (see e.g. [12] or [16]). If moreover card(supp PX ) = ∞, then x∗ has pairwise distinct ˆ x ||p is decreasing to 0 as N goes to ∞. components and minx∈(Rd )N || X − X To be more precise, the Zador Theorem (see [12]) provides a rate for the optimal Lp -mean 0 quantization error, namely, if X ∈ Lp for some p0 > p, ˆ x ||p = O(N −1/d ) min || X − X x∈(Rd )N as N → ∞. (5) We refer to the following papers for a study of several algorithms designed to find optimal quantizers: (see [15, 11, 16, 17]). The problem dimension and the properties of the law of X might help to determine an efficient algorithm. Note that some optimized quantizers of the normal distribution N (0, Id ) have been computed and are available at the URL www.quantize.math-fi.com for dimensions up to 10 and sizes up to several thousands. ˆ x provides the best finite approximation to the distribution of Since an optimal quantization X ˆ x ) as an approximation for Ef (X), X in the least square sense, it becomes natural to use Ef (X where f : Rd → R is a Borel function. 3 ˆ x takes only finitely many values, we compute Ef (X ˆ x ) as the finite Note further, that since X sum N X P(X ∈ Ci (x))f (xi ). i=1 The weights pi := P(X ∈ Ci (x)) of this cubature formula are obtained as a by-product of optimization procedures to generate optimal quantizers like the ones used in [17] or [16]. Assume now that f exhibits some smoothness properties, i.e. f is differentiable with Lipschitz continuous differential Df . If x is a stationary quantizer, i.e. ˆ x = E(X|X ˆ x ), X (so is the case for every optimal L2 -quantizer), we can conclude from a second order Taylor expanˆ x ) satisfies sion of f , that the approximation error of Ef (X hal-00273790, version 2 - 15 Dec 2009 ˆ x )| ≤ [Df ]Lip kX − X ˆ x k22 . |Ef (X) − Ef (X Similarly, we can derive (see e.g. [18]) some error bounds for the approximation of E(f (X)|Y ) ˆ x )|Yˆ y ). by its optimal quantized counterpart E(f (X 3 The quantization tree algorithm 3.1 Abstract problem formulation We briefly recall the general stochastic control framework developed in [2] to model “abstract” swing options and the quantization tree algorithm devised to solve it. Let (tk )0≤k≤n−1 be a given, strictly increasing sequence of exercise dates, such that the first exercise date is set at the origin of the derivative and the last one happens strictly before the option maturity. Let (Ω, A, P) be a complete probability space and let p ∈ [1, ∞[. On this probability space, one defines an option with maturity T by a sequence of random variables (Vk )0≤k≤n−1 . For k ∈ {0, . . . , n−1}, Vk represents the option payoff obtained by the derivative holder at time tk . To model all the available market information, one introduces a filtration F := (Fk )0≤k≤n−1 on (Ω, A, P), such that, for k ∈ {0, . . . , n − 1}, the option payoff Vk is measurable with respect to the σ-field Fk . One common feature characterizing most of the options is that the product is built on an underlying asset. In our setting, the option owner has the possibility to choose, for each exercise date, the underlying asset amount included in the contract. These decision variables are represented by a sequence (qk )0≤k≤n−1 of F-adapted random variables and are subjected to constraints given by: (i) a local constraint: let qmin and qmax be two given positive constants. For k = 0, . . . , n − 1, qmin ≤ qk ≤ qmax P-a.s.. (ii) A global constraint: one defines a cumulative purchased process (¯ qk )0≤k≤n by q¯0 := 0 and q¯k := k−1 X l=0 4 ql , k = 1, . . . , n. (6) Then for any couple of F0 -measurable, non-negative random variables Q := (Qmin , Qmax ), the total cumulative volume purchased should satisfy Qmin ≤ q¯n ≤ Qmax P-a.s.. (7) A control process (qk )0≤k≤n−1 which satisfies these local and global conditions is called (F, Q)admissible. Now, one defines the residual global constraints Qk = (Qkmin , Qkmax ) at time tk by Qkmin = Qmin − q¯k and Qkmax = Qmax − q¯k . Then, the associated option value at time tk is given by ( ! ) n−1 X Pkn ((Qkmin , Qkmax )) := esssup E ql Vl |Fk , ql : (Ω, Fl ) → [qmin , qmax ], k ≤ l ≤ n − 1, q¯n − q¯k ∈ [Qkmin , Qkmax ] . hal-00273790, version 2 - 15 Dec 2009 l=k (8) Now, to simplify notations, we perform a decomposition of the problem into a simple “swap” part and a “normalized” one, where the control variables are [0, 1] valued, i.e. q = qmin + (qmax − qmin )q 0 for a [0, 1] valued r.v. q 0 , which leads to Pkn ((Qkmin , Qkmax )) = qmin n−1 X ˜k , Q ˜ k )) E (Vl |Fk ) + (qmax − qmin )Pk[0,1],n ((Q min max (9) l=k with k k ˜ k = Qmin − (n − k)qmin and Q ˜ k = Qmax − (n − k)qmin . Q min max qmax − qmin qmax − qmin One considers from now on a normalized framework, i.e. qmin = 0 and qmax = 1. In this setting, a couple of global constraints Qk at time tk is admissible if Qk is Fk -measurable and if 0 ≤ Qkmin ≤ Qkmax ≤ n − k P-a.s.. (10) For every couple of admissible global constraints Qk , the derivative price at time tk satisfies ( ! ) n−1 X ql Vl |Fk , ql : (Ω, Fl ) → [0, 1], k ≤ l ≤ n − 1, q¯n − q¯k ∈ [Qkmin , Qkmax ] . Pkn ((Qkmin , Qkmax )) := esssup E l=k (11) To design the quantization tree algorithm, we appeal to several results summarized below (see [2] for further details): (i) Without loss of generality, one may assume that at time tk the couple of admissible global constraints Qk is deterministic. (ii) The Backward Dynamic Programming Formula (BDP): Set Pnn ≡ 0. For every k ∈ {0, . . . , n−1} and every couple of admissible global constraints Qk = (Qkmin , Qkmax ) at time tk , we have n o n−k−1 n Pkn (Qk ) = sup xVk + E(Pk+1 (χn−k−1 (Qk , x))|Fk ), x ∈ IQ k with χM (Qk , x) := (Qkmin − x)+ , (Qkmax − x) ∧ M and M k + IQ ∧ 1, Qkmax ∧ 1]. k := [(Qmin − M ) 5 (iii) Bang-bang control: One defines the set of admissible deterministic global constraints at time 0 by T + (n) = {Q0 ∈ R2 | 0 ≤ Q0min ≤ Q0max ≤ n}. (12) Then, given any admissible integral global constraints at time 0 i.e. Q0 ∈ N2 ∩ T + (n), there exists an optimal bang-bang control process q ∗ = (qk∗ )0≤k≤n−1 with qk∗ ∈ {0, 1} P-a.s. for every k = 0, . . . , n − 1. hal-00273790, version 2 - 15 Dec 2009 (iv) The value function Q0 7→ P0n (Q0 ) is concave, continuous and piecewise affine on the tiling of T + (n). So in view of these results, it is sufficient to evaluate the option premium at integral values of T + (n). Then, it is possible to use a linear interpolation inside every tile to obtain the derivative price for all admissible constraints, i.e. for all Q0 ∈ T + (n). Hence we may reformulate the BDP Formula for an initial Q0 ∈ N2 ∩ T + (n) as Pnn ≡ 0 o n n−k−1 n , Pkn (Qk ) = max xVk + E(Pk+1 (χn−k−1 (Qk , x))|Fk ), x ∈ {0, 1} ∩ IQ k k = 0, . . . , n − 1. (13) In order to simulate the derivative prices obtained in (13), one appeals to the quantization method described in Section 2. At each exercise date tk , k = 0, . . . , n − 1, we perform a spatial discretization of the d-dimensional underlying asset and include the resulting discretized process in the BDP Formula. Finally, some convergence theorems are available to obtain some error bounds. However, these theorems require the following assumption: Assume that there exists a d-dimensional Markov structure process, say (Xk )0≤k≤n−1 , such that the payoffs Vk are function of Xk , i.e. Vk = vk (Xk ), k = 0, . . . , n − 1. (14) One obtains E(vk+1 (Xk+1 )|Fk ) = E(vk+1 (Xk+1 )|Xk ) = Θk (vk+1 )(Xk ), k = 0, . . . , n − 1, (15) where (Θk )0≤k≤n−1 is a sequence of Borel transition probabilities on (Rd , B(Rd )). Therefore (13) becomes Pnn ≡ 0 n o n−k−1 n Pkn (Qk ) = max xvk (Xk ) + E(Pk+1 (χn−k−1 (Qk , x))|Xk ), x ∈ {0, 1} ∩ IQ , k k = 0, . . . , n − 1. (16) ˆ k be a quantized version of Xk of size say Nk . An approximation of the price process Now, let X ˆ k in (16) and by forcing the Markov property on X ˆk . is designed by plugging X Pˆnn ≡ 0 n o ˆ k ) + E(Pˆ n (χn−k−1 (Qk , x))|X ˆ k ), x ∈ {0, 1} ∩ I n−k−1 Pˆkn (Qk ) = max xvk (X , k k+1 Q k = 0, . . . , n − 1. (17) 6 The resulting algorithm, called quantization tree, is then tractable on a computer. Indeed, it is possible to evaluate the above formula in a recursive backward way down to k = 0 where we arrive at a σ(X0 )-measurable solution. Note that if F0 = {∅, Ω}, the option premium is deterministic. From now on, we assume that F0 is trivial. The following results are concerned with the error resulting from the discretization of the underlying price process (Xk )0≤k≤n−1 (see [2]). Theorem 3.1. Assume that the Markov process (Xk )0≤k≤n−1 is Lipschitz Feller in the following sense: for every bounded Lipschitz continuous function g : Rd → R and every k ∈ {0, . . . , n − 1}, Θk (g) is a Lipschitz continuous function satisfying [Θk (g)]Lip ≤ [Θk ]Lip [g]Lip . Assume that every function vk : Rd → R is Lipschitz continuous with Lipschitz coefficient [vk ]Lip . Let p ∈ [1, ∞[ such that max0≤k≤n−1 |Xk | ∈ Lp (P). Then, there exists a real constant Cp > 0 such that k sup hal-00273790, version 2 - 15 Dec 2009 Q∈T + (n)∩N2 |P0n (Q) − Pˆ0n (Q)|kp ≤ Cp n−1 X ˆ k kp . kXk − X k=0 In view of Zador’s Theorem (5) and the F0 -measurability of P0n and Pˆ0n , this also reads sup Q∈T + (n)∩N2 3.2 |P0n (Q) − Pˆ0n (Q)| ≤ C n−1 X −1/d Nk . k=0 Complexity and implementationary notes For an analysis of the complexity of the quantization tree algorithm in its original form (17), we count the number of multiplications, which occur during the evaluation of P0n (Q0 ) for a given Q0 . We implement the quantization tree algorithm in a backward iterative manner, starting at layer k with the computation of Pkn (Qk ) for every possible residual global constraint Qk ∈ Qnk (Q0 ) with Qnk (Q0 ) = (Q0min − l)+ , (Q0max − l)+ ∧ (n − k) , l = 0, . . . , k for given initial global constraints Q = (Qmin , Qmax ) ∈ T + (n) ∩ N2 . This approach yields a complexity proportional to n−2 X #Qnk (Q0 )Nk Nk+1 + Nn−1 k=0 multiplications, which can be estimated from above by n−2 X (k + 1)Nk Nk+1 + Nn−1 . k=0 If we employ an equal grid size N in each layer for the quantization of the conditional expecta2 2 tions, we finally arrive at an upper bound proportional to n 2N multiplications. 3.3 The call strip: case Q = (0, n) We will study a special case of swing options. This class is characterized by particular values of global constraints. In this case the global constraint is always satisfied by the cumulative purchase process, i.e. Pn−1 l=0 ql ∈ [Qmin , Qmax ], as the control process (ql )0≤l≤n−1 is [0, 1]-valued. Hence, the optimal 7 control problem is reduced to a maximization problem without constraints: ( ! ) n−1 X Pkn (0, n) = esssup E ql Vl |Fk , ql : (Ω, Fl ) → [0, 1], k ≤ l ≤ n − 1 , k = 0, . . . , n − 1. l=k (18) An optimal control for (18) is clearly given by ( 1 on {Vl ≥ 0}, ql∗ := 0 elsewhere, and the associated price process by n−1 X Pkn (0, n) = E Vl + |Fk . hal-00273790, version 2 - 15 Dec 2009 l=k Calling upon the results of Section (3.1), one obtains the following quantized tree algorithm in a Markovian structure n−1 X E (vl (Xˆ l ))+ |Xˆ k . Pˆkn (0, n) = l=k Note that this quantity appears as a sum of European options. Since closed forms are available to evaluate European options, the call strip will be used as benchmark for numerical implementations. Remark 1. The case Q = (0, 1) leads in the same way to the payoff of a Bermudan option. 4 Two main improvements for the quantization tree algorithm 4.1 Computations of the transition probabilities: fast nearest neighbor search Solving the quantized BDP principle (17) amounts from a numerical point of view to compute for some Borel function f : Rd → R all the conditional expectations E(f (Xˆ k+1 )|Xˆ k ), which take values in the finite grid xk = (xk1 , . . . , xkNk ) of Rd . Moreover, we derive that Nk+1 E(f (Xˆ k+1 )|Xˆ k = xki ) = X f (xk+1 )πkij j j=1 where ˆ k+1 = xk+1 |X ˆ k = xk ), πkij = P(X i j k = 0, . . . , n − 2, denote the transition probabilities of the quantized process. By definition of the Voronoi quantizaˆ in (4), we arrive at tion X πkij = P(Xk+1 ∈ Cj (xk+1 )|Xk ∈ Ci (xk )) = P(Xk+1 ∈ Cj (xk+1 ), Xk ∈ Ci (xk )) . P(Xk ∈ Ci (xk )) (19) To estimate the transition probabilities by means of a naive Monte Carlo simulation, one gener˜ m )1≤m≤M , k = ates samples of the Markov chain (Xk )0≤k≤n−1 with size say M , which we denote by (X k ij 0, . . . , n − 1. So πk can be estimated by PM ˜m ˜m m=1 1Cj (xk+1 ) (Xk+1 )1Ci (xk ) (Xk ) ij π ˜k := , for k = 0, . . . , n − 2, (20) PM ˜ m) 1C (xk ) (X m=1 k i i = 1, . . . , Nk and j = 1, . . . , Nk+1 . 8 hal-00273790, version 2 - 15 Dec 2009 ˜ m (ω) and X ˜ m (ω) Given k, i and j, one has to check for each m = 1, . . . , M if the realizations X k k+1 k k+1 belong to the cells Ci (x ) or Cj (x ) of the Voronoi tessellation to determine the value of the indicator functions present in the estimations (20). Nearest neighbor search: This problem consists in finding the nearest neighbor for a query point q ∈ Rd in the N -tuple x = (x1 , . . . , xN ) with regard to the Euclidean distance. A commonly known algorithm to optimize this nearest neighbor search, which in high dimensions can become a time consuming procedure, is the kd-tree algorithm (see e.g. [8]). This procedure consists of two steps: The partitioning of the space, which has an average complexity cost of O(n log n) operations. This has to be done only once and is therefore independent of the size M of the Monte Carlo sample. The second part is the search procedure, which has an average complexity cost of O(log n) distance computations. The kd-tree data structure is based on a recursive subdivision of space Rd into disjoint hyperrectangular regions called boxes. Each node of the resulting tree is associated with such a region. To ensure a well-defined splitting rule on each node, we equip the root node with a bounding box, which is chosen to contain all data points. As long as the number of data points within such box is greater than a given threshold, the box is split-up into two boxes by an axis-orthogonal hyperplane that intersects the original box. These two boxes are then associated to two children nodes. When the number of points in a box is below the threshold, the associated node is declared a leaf node, and the points are stored within this node. To find the nearest neighbor for a query point q, we start at the root node and decide at each node due to the cutting hyperplane, to which child node we descend. When arriving at a leaf node, we compute the minimal distance from all data points of this node to the query point q and ascend back on the trial to the root node as long as it may be possible to find closer points than the previous found minimal distance. During this way back we obviously descend to the so far not visited children and update the minimal distance found so far. From now on, we assume for convenience that (Xk )0≤k≤n−1 can be written recursively, i.e. Xk+1 = Ak Xk + Tk Zk , k = 0, . . . , n − 2, X0 = x0 , (21) for some positive definite Matrices Ak and Tk , and a deterministic initial value x0 . Furthermore we suppose, that the distribution of Zk is either known through its density or can be simulated at some reasonable costs and that Zk is independent of the distribution of Xl , l = 0, . . . , k. These assumptions hold for example for the Black-Scholes model (where the driving process is a Brownian motion), the Ornstein-Uhlenbeck process, which we will later on represent as a Gaussian first order auto-regressive process or for certain L´evy processes, e.g. the NIG process. 4.2 4.2.1 Parallelization of the transition probability estimation Diffusion method Using the representation (21) of X as an auto-regressive process, the naive Monte Carlo approach from (20) now reads as follows. Generate a Monte Carlo sample of size M of the random variate Zk in each layer k, i.e m (Z˜0m )1≤m≤M , . . . , (Z˜n−2 )1≤m≤M , and define the Monte Carlo sample for (Xk )0≤k≤n−1 by ˜ 0m := x0 , X ˜ m := Ak X ˜ m + Tk Z˜ m , X k+1 k k 9 k = 0, . . . , n − 2 (22) for m = 1, . . . , M . Then we get an estimate for πkij , which writes PM π ˜kij := ˜m ˜m m=1 1Cj (xk+1 ) (Xk+1 )1Ci (xk ) (Xk ) . PM ˜m m=1 1Ci (xk ) (Xk ) hal-00273790, version 2 - 15 Dec 2009 It is possible to implement a parallel procedure of the diffusion method like for any “linear” Monte Carlo simulation. However the computational structure of the transition probabilities is not well adapted to this parallel implementation. Indeed, in this method, we have to proceed the ˜ m , that is, step by step from time 0 up to the computations along the discrete trajectories of X exercise date tn−1 . Therefore, the natural way to partition the tasks would be to divide the number of Monte Carlo simulations into several processes. Once the computations executed simultaneously by each process are completed and sent back to the main process, this process still would have to add up the partial sums for each transition in each layer k, which leads to some additional and significant work. So, a more advantageous approach from the parallelization viewpoint will be described in the section below. 4.2.2 Parallel Quantization Weight Estimation This method was introduced in [3] and heavily relies on the assumptions on X satisfying the auto-regressive structure (21). We generate in each layer k, independently, bivariate Monte Carlo samples by simulating directly from the distributions of Xk and Zk , i.e. ˜ m , Z˜ m )1≤m≤M , . . . , (X ˜ m , Z˜ m )1≤m≤M . (X 0 0 n−2 n−2 The estimated transition probabilities are then given by PM ˜m ˜m ˜ ˜m m=1 1Cj (xk+1 ) (Ak Xk + Tk Zk )1Ci (xk ) (Xk ) ij π ˜k := . PM ˜m m=1 1Ci (xk ) (Xk ) In this approach, the computations of the transitions for a layer k are completely independent from all the preceding and succeeding layers. So a parallel implementation with respect to the time layers tk becomes very straightforward. Additionally, all the ideas for a further parallel split-up with regard to the number of Monte Carlo simulations, as for the Diffusion method, remain valid. 4.2.3 Spray method The spray method is not an exact approach, since it aims at estimating only an approximation of ˆ k in the πkij . This means that we not only replace the random variable Xk by its quantization X k identity Xk+1 = Ak Xk + Tk Zk , but also in the conditional part Xk ∈ Ci (x ) of (19). ˆ k ∈ Ci (xk ) and X ˆ k = xk , this leads to the following approximation Using the equivalence of X i π ˜kij ˆ k + Tk Zk ∈ Cj (xk+1 )|X ˆ k ∈ Ci (xk )) := P(Ak X = = P(Ak xki + Tk Zk ∈ Cj (xk+1 )) P Zk ∈ Tk−1 (Cj (xk+1 ) − Ak xki ) . In low dimensions and when the density of Zk is known, this quantity can be computed by deterministic integration methods. Otherwise a Monte Carlo estimate would be given by π ˜kij := M 1 X 1T −1 (Cj (xk+1 )−Ak xk ) (Z˜km ), i k M m=1 10 where m (Z˜0m )1≤m≤M , . . . , (Z˜n−2 )1≤m≤M denote again the i.i.d. Monte Carlo samples of (Zk ). Note that such an approach can be extended to a non linear dynamic Xk+1 = F (Xk , Zk ), Zk i.i.d.. 4.3 Improving the convergence rate: a Romberg extrapolation approach Consider for instance the classical numerical integration by means of optimal quantization, i.e. Ef (X) is approximated by Ef (Xˆ N ) for sequence of stationary and rate optimal quantizers Xˆ N (see Section 2). Assume that f : Rd → R is a 3-times differentiable function with bounded derivatives. We also assume (which is still a conjecture) that hal-00273790, version 2 - 15 Dec 2009 E(D2 f (Xˆ N )(X − Xˆ N )⊗2 ) ∼ cN −2/d as N → ∞ and that E|X − Xˆ N |3 = O(N −3/d ). Then we can conclude from a Taylor expansion ˆ N )+Df (X ˆ N )(X−X ˆ N )+ 1 D2 f (X ˆ N )(X−X ˆ N )⊗2 + 1 D3 f (ξ)(X−X ˆ N )⊗3 , f (X) = f (X 2 6 ˆ N , X), ξ ∈ (X so that Ef (X) = Ef (Xˆ N ) + cN −2/d + o(N −3/d+ε ) for every ε > 0. Motivated by this identity and our numerical tests, we make the conjecture, that an analogous result holds also (at least approximately) true for the quantization tree algorithm P0n = Pˆ0n + cN −2/d + o(N −3/d+ε ) for every ε > 0. Hence we perform a Romberg extrapolation on the quantized swing option prices by computing Pˆ0n for two different values N1 , N2 (e.g. N1 ≈ 4N2 ) and denote them by Pˆ0n (N1 ) and Pˆ0n (N2 ). Thus we arrive at Pˆ n (N1 ) − Pˆ0n (N2 ) −2/d P0n = Pˆ0n (N1 ) + 0 −2/d N1 + o(N −3/d+ε ). −2/d N2 − N1 This suggests to consider as improved approximation for P0n the quantity P0n (N1 , N2 ) := Pˆ0n (N1 ) + 5 Pˆ0n (N1 ) − Pˆ0n (N2 ) −2/d N2 − −2/d N1 −2/d N1 . Application: swing options We will consider in the section three different ways of modeling for the underlying dynamics, which all fulfill the assumptions of the quantization tree algorithm. These are namely the 1- and 2-factor Gaussian model and an exponential L´evy model. From now on, we assume an equidistant time discretization 0 = t0 < t1 < . . . < tn = T , i.e. tk = k∆t, k = 0, . . . , n with ∆t := T /n. 11 5.1 Gaussian 1-factor model In this case, we consider a dynamic of the underlying forward curve (Ft,T )t∈[0,T ] given by the SDE dFt,T = σe−α(T −t) Ft,T dWt , which yields Z t St = F0,t exp(σ 0 for hal-00273790, version 2 - 15 Dec 2009 ∆2t = 1 e−α(t−s) dWs − ∆2t ) 2 σ2 1 − e−2αt . 2α R k∆t If we set Xk := 0 e−α(k∆t−s) dWs , k = 0, . . . , n − 1, (Xk )0≤k≤n−1 is a discrete time Markov process and we get Sk∆t −K = vk (Xk ) with Lipschitz-continuous vk (x) = (F0,k∆t exp(σx− 12 ∆2k∆t )− K). So the formal requirements for the application of the quantization tree algorithm are fulfilled. For the quantization and the computation of the transition probabilities of Xk we need the following results: R k∆t Proposition 5.1. The discrete time Ornstein-Uhlenbeck process Xk = 0 e−α(k∆t−s) dWs can be written as first order auto-regressive process r 1 − e−2α∆t Xk+1 = e−α∆t Xk + k , k = 0, . . . , n − 2, X0 := 0, 2α for an i.i.d sequence (k ), 1 ∼ N (0, 1). Especially we have 1 − e−2αtk . Xk ∼ N 0, 2α Since affine transformations of a one-dimensional random variate can be transformed one-toone on its optimal quantizers, the quantizers for Xk can be constructed as a dilatation of optimal q quantizers for N (0, 1) by the factor 5.2 1−e−2αtk . 2α Gaussian 2-factor model Furthermore we have also considered a Gaussian 2-factor model, where the dynamic of the forward curve (Ft,T )t∈[0,T ] is given by the SDE dFt,T = Ft,T σ1 e−α1 (T −t) dWt1 + σ2 e−α2 (T −t) dWt2 for two Brownian motions W 1 and W 2 with correlation coefficient ρ. This yields Z t Z t 1 2 −α1 (t−s) 1 −α2 (t−s) 2 St = F0,t exp σ1 e dWs + σ2 e dWs − ∆t 2 0 0 for σ12 σ2 σ1 σ2 1 − e−2α1 t + 2 1 − e−2α2 t + 2ρ 1 − e−(α1 +α2 )t . 2α1 2α2 α1 + α2 In this case we have to choose Z k∆t Z k∆t Xk := e−α1 (k∆t−s) dWs1 , e−α2 (k∆t−s) dWs2 , k = 0, . . . , n − 1, ∆2t = 0 0 12 as underlying Markov structure process with 1 vk (x1 , x2 ) = (F0,k∆t exp(σ1 x1 + σ2 x2 − ∆2k∆t ) − K). 2 Applying Proposition 5.1 on the two components of (Xk )0≤k≤n−1 allows us to write it as a first order auto-regressive process Xk+1 = Ak Xk + Tk k with A= and q T = q 1 2α1 1 2α2 hal-00273790, version 2 - 15 Dec 2009 (1 − e−2α1 ∆t ) r = ρq 1 α1 +α2 1 4α1 α2 (1 0 (1 − e−2α2 ∆t )r where 5.3 e−α1 ∆t 0 −α 0 e 2 ∆t q 1 2α2 √ (1 − e−2α2 ∆t ) 1 − r2 1 − e−(α1 +α2 )∆t . − e−2α1 ∆t )(1 − e−2α2 ∆t ) Normal Inverse Gaussian spot return In this case we assume, that the dynamics of the underlying is driven by the exponential of a special L´evy process, the so-called Normal Inverse Gaussian process (NIG). This means, that we assume, that the dynamic of the underlying is given by St = S0 exp(Lt ) for a NIG process (Lt )t∈[0,T ] . This special model for financial data has been first proposed by Barndorff-Nielsen in [4] and has been later applied to financial contracts on energy markets (see e.g. [5] and [6]) as well. The NIG process has beneath its L´evy property, which makes it a Markov process, the convenient property, that it is completely determined by its distribution at time t =1, the so-called NIG distribution, NIG(α, β, δ, µ), with parameters α > 0, |β| < α, δ > 0 and µ ∈ R, i.e. L1 ∼ NIG(α, β, δ, µ) and Lt ∼ NIG(α, β, tδ, tµ) Thus, as soon as we have knowledge on some properties of the NIG distribution like the density or the characteristic function, we know it already for (Lt )t∈[0,T ] at any timescale. So, the density of (Lt )t∈[0,T ] is for example given by NIG(α,β,δ,µ) ft (x) √ tδ = αδte p t2 δ 2 + (x − tµ)2 ) p , π t2 δ 2 + (x − tµ)2 α2 −β 2 +β(x−tµ) K1 (α where K1 denotes the Bessel function of third kind with index 1. These facts will help to the computation of optimal quantizers for Lt at given time-points {k∆t, k = 0, . . . , n − 1}. Another useful property of the NIG process is the fact that if we model our underlying as St = S0 exp(Lt ) 13 with respect to the real-world measure P, we can use the Esscher transform to construct an equivalent martingale measure Q (see [9] and [10]), for a first occurrence of this method, which preserves the NIG structure of the driving L´evy process. Hence the probability change from P to Q can be completely performed by adjusting the parameters α, β, δ and µ of the NIG process, so that from a numerical point of view it makes no difference, if we are modeling under the real-world measure or the risk-neutral one. 5.3.1 Computation of optimal quantizers for the NIG distribution Recall that we may rewrite the L2 -quantization problem of the NIG distribution as N Z X DN := min (ξ − xi )2 f NIG (ξ)dξ hal-00273790, version 2 - 15 Dec 2009 x∈RN i=1 (23) Ci (x) where {Ci (x), i = 1, . . . , N } denotes a Voronoi partition induced by x = (x1 , . . . , xN ). The fact, that in the one-dimensional setting the Voronoi cells Ci (x) are just intervals in R and the uniqueness of the optimal quantizer due to the unimodality of f NIG , makes the quantization problem in this case very straightforward to solve. Now assume x = (x1 , . . . , xN ) ∈ RN to be ordered increasingly and denote by xi + xi±1 for 2 ≤ i ≤ N − 1, xN +1/2 := +∞ 2 the midpoints between the quantizer elements respectively ±∞. A Voronoi partition of Γ is therefore given by x1/2 := −∞, xi±1/2 := C1 (x) = (−∞, x1+1/2 ], Ci (x) = (xi−1/2 , xi+1/2 ], 2 ≤ i ≤ N − 1, CN (x) = (xN −1/2 , +∞). so that (23) now reads min x∈RN N Z X i=1 xi+1/2 (ξ − xi )2 f NIG (ξ)dξ. (24) xi−1/2 This is an N -dimensional optimization problem, which can be easily solved by Newton’s method as soon as we have access to the first and second order derivatives of DN . In fact, we can calculate the first order derivative of DN (see e.g. [12], Lemma 4.10 or [19], Lemma C) as Z xi+1/2 ∂DN (x) = 2 (xi − ξ) f NIG (ξ)dξ, 1 ≤ i ≤ N. (25) ∂xi xi−1/2 Moreover the Hessian matrix turns out to be a symmetric tridiagonal matrix with diagonal entries Z x1+1/2 ∂ 2 DN x2 − x1 NIG (x) = 2 f NIG (ξ)dξ − f (x1+1/2 ) 2 2 ∂x1 x1/2 Z xi+1/2 ∂ 2 DN xi+1 − xi NIG xi − xi−1 NIG (x) = 2 f NIG (ξ)dξ − f (xi−1/2 ) − f (xi+1/2 ), 2 ≤ i ≤ N 2 2 2 ∂xi xi−1/2 Z xN +1/2 ∂ 2 DN xN − xN −1 NIG (x) = 2 f NIG (ξ)dξ − f (xN −1/2 ) 2 2 ∂xN xN −1/2 (26) 14 and the super- respectively sub-diagonals ∂ 2 DN ∂ 2 DN xi − xi−1 NIG (x) = (x) = − f (xi−1/2 ), ∂xi ∂xi−1 ∂xi−1 ∂xi 2 2 ≤ i ≤ N. As a consequence of (25), the remaining entries vanish ∂ 2 DN ∂ 2 DN (x) = (x) = 0, ∂xi ∂xj ∂xj ∂xi 1 ≤ i < j − 1 ≤ N − 1. hal-00273790, version 2 - 15 Dec 2009 For the evaluation of the integrals occurring in the expressions (25) and (26) we employed high-order numerical integration methods which gave satisfying results. But some attention should be paid to the initialization of the Newton’s method, since it converges only locally. We achieved a fast convergence using an equidistant placed N -tuple in the interval [ELtk − 2 Var Ltk , ELtk + 2 Var Ltk ] as starting vector, i.e. " ! !# δβ δα2 δβ δα2 tk µ + p −2 2 , tk µ + p +2 2 (α − β 2 )3/2 (α − β 2 )3/2 α2 − β 2 α2 − β 2 for the case of the NIG, so that we reached the stopping criterion of k∇DN (x∗ )k ≤ 10−8 within 20 − 30 iterations. 5.3.2 Simulation of ∆Lt For the simulation of a NIG(α, β, δ, µ) process it is useful to regard Lt as subordinated to a Brownian motion, i.e. Lt = βδ 2 It + δWIt + µt, for an Inverse Gaussian process It , that p is again a L´evy process with Inverse Gaussian distribution, IG(a, b), and parameters a = t and b = α2 − β 2 . A way to simulate the IG(a, b) distribution was proposed in [7], p.182. 6 Numerical results In the above explained models we performed some tests on the pricing of swing options. 6.1 The pricing procedure Given any initial global constraints Q ∈ N2 ∩T + (n), the quantized tree algorithm aims at backward solving the following program Pˆnn ≡ 0 n o ˆ k ) + E(Pˆ n (χn−k−1 (Qk , x))|X ˆ k ), x ∈ {0, 1} ∩ I n−k−1 Pˆkn (Qk ) = max xvk (X , k+1 Qk k = 0, . . . , n − 1, (27) for all admissible Qk from section 3.1. Here, vk is a local payoff with strike K, which is defined by 1 vk (x) = (F0,k∆t exp(σx − ∆2k∆t ) − K), 2 for the Gaussian 1-factor model, 1 vk (x1 , x2 ) = (F0,k∆t exp(σ1 x1 + σ2 x2 − ∆2k∆t ) − K), 2 15 for the Gaussian 2-factor model and vk (x) = S0 exp(x) − K for the NIG-model. Note that the only admissible Q0 is Q, and therefore P0n (Q) := P0n (Q0 ) gives then the premium of a swing option. In all test cases, we set the local constraints to qmin = 0 and qmax = 6 and give results for a number of exercise dates n = 30 resp. n = 365 and strikes K = 5, 10, 15, 20. hal-00273790, version 2 - 15 Dec 2009 As a benchmark served the call strip case Q = (0, qmax · n), ˜ = (0, n) and therefore has a reference price which reduces in view of (9) to the normalized case Q given by a sum of plain vanilla calls (see section 3.3). The tests were carried out in JAVA on a dual Intel Xeon QuadCore [email protected] and the CUDAresults on a NVIDIA GTX 280 GPU. Moreover, the Monte-Carlo sample size were in all cases chosen as M = 105 per time layer tk . 6.2 1-factor model For the Gaussian 1-factor model we have chosen the parameters σ = 0.7, α = 4, F0,t = 20 with respect to a 30-day period for the underlying dynamic. We give results for the diffusion method (diff), as representative for the MC-based approaches, and the deterministic spray method (dSpray), which uses a 63-point Gauss quadrature for the numerical integration. In addition, we implemented the trinomial tree method from [14], which is based on a tree dynamics for a mean reverting process from [13]. In fact, this approach fits exactly into our pricing framework, where only the grid points and the 3-way transitions have to be chosen accordingly to [14], section 4.1. Unfortunately, this design is only defined for a special number of grid points, which is 15 in our parameter setting. Hence, we could carry out comparisons to the trinomial tree method only for the case N = 15. Concerning the accuracy of these methods, we start with the call strip-case in Table 1. 16 N 15 50 100 200 diff 1799.02 1800.37 1800.49 1800.51 rel. Err. 0.073% 0.003% 0.009% 0.010% sec. 1.52 1.65 1.74 1.86 dSpray 1797.54 1799.97 1800.23 1800.30 N 15 50 100 200 diff 319.58 320.13 320.38 320.40 rel. Err. 0.210% 0.037% 0.041% 0.048% sec. 1.52 1.65 1.74 1.86 dSpray 316.61 319.57 320.10 320.21 K = 10 rel. Err. 0.155% 0.020% 0.005% 0.001% K = 20 rel. Err. 1.137% 0.211% 0.046% 0.013% sec. 0.08 0.45 1.65 6.48 trinom 1808.06 rel. Err. 0.430% sec. 0.001 sec. 0.08 0.45 1.65 6.48 trinom 327.36 rel. Err. 2.218% sec. 0.001 hal-00273790, version 2 - 15 Dec 2009 Table 1: Gaussian 1-factor model for n = 30 and Q = (0, 180) with computational times for the transition probabilities on a single CPU. For both strikes K = 10, 20 the quantization methods induce already for N = 15 a much smaller approximation error than the trinomial tree method, since they are based on a discretization which is fitted in an optimal way to the underlying distribution, in contrast to the trinomial tree method which is based on an equidistant discretization. Moreover, with a grid size of N = 50, both quantization methods provide relative errors in the region of a few ‰ or less, which is nearly the exact result. Since the diffusion approach contains, compared to the spray method, one approximation less, it outperforms the deterministic spray method for smaller N as shown in the table. In view of the computational time for the estimation of the transition probabilities, the trinomial tree approach, where each node can only take 3 different states, performs very fast. Nevertheless, the deterministic spray method executes for N = 15 in less than 1/10 second and for N = 50 in less than 1/2 second, which is still instantaneous and provides moreover a much higher accuracy than the trinomial tree method. Concerning MC-simulation based approaches, we here only give computation times for the diffusion approach, which is about 1 21 second and therefore cannot compete with the other approaches for small and moderate N . Once the computation of the transition probabilities is done, the pricing method has to iteratively traverse the quantization tree for the solution of the stochastic control problem in (27), whose execution time are given in Table 2. N sec. 15 0.02 50 0.10 100 0.17 200 0.29 Table 2: Computational times on a single CPU for the traversal of the quantization tree in the 1-factor model. These times, which are independent of the transition approach and which can be considered as the minor compute intensive problem, have to be added to those of Table 1 to arrive at the total time for the pricing of one contract. To illustrate the performance of these methods in case of a non-trivial control problem, we chose global constraints Q = (100, 150) in Table 3. 17 N 15 50 100 200 diff 1589.97 1588.95 1588.89 1588.86 N 15 50 100 200 diff 226.83 225.28 225.19 225.15 K = 10 dSpray 1587.66 1588.41 1588.51 1588.54 K = 20 dSpray 223.76 224.75 224.90 224.94 trinom 1602.48 trinom 241.29 hal-00273790, version 2 - 15 Dec 2009 Table 3: Gaussian 1-factor model for n = 30 and Q = (100, 150). Here again, the quantization methods are very close to each other and suggest that the true price for K = 20 is around 225, which means that these methods again outperform the trinomial tree method by several orders of magnitude. Moreover, the consistency with Table 1 justifies the use of the call strip case as a general benchmark for the performance. In all these examples the dSpray-method with N = 50 performs the pricing of a 30-day contract in about 1/2 second and with an error of only a few ‰, which seems to be a good compromise between speed and accuracy. 6.3 2-factor model For the numerical analysis in the Gaussian 2-factor model, we set parameters σ1 = 0.36, α1 = 0.21, σ2 = 1.11, α2 = 5.4, ρ = −0.11, F0,t = 20 with respect to a one year period. Since the 2-factor model is based on bivariate normal distribution in each time layer, deterministic approaches are not competitive anymore, and we therefore focus only on the MC-based methods, i.e. the diffusion method (diff), the parallel Quantization Weight Estimation-method (pQWE) and the spray method for MC-simulation (MCSpray). Starting with a 30-day contract in Table 4, N 100 250 500 diff 1796.57 1797.58 1797.94 rel. Err. 0.202% 0.146% 0.126% sec. 3.84 4.64 5.26 N 100 250 500 diff 263.72 265.50 266.16 rel. Err. 1.816% 1.153% 0.906% sec. 3.84 4.64 5.26 K = 10 rel. Err. 0.108% 0.016% 0.041% K = 20 pQWE rel. Err. 264.38 1.568% 266.64 0.728% 267.70 0.334% pQWE 1798.26 1799.91 1800.94 sec. 4.78 5.21 5.57 MCSpray 1793.11 1794.72 1795.98 rel. Err. 0.394% 0.305% 0.235% sec. 2.56 2.67 2.92 sec. 4.78 5.21 5.57 MCSpray 254.55 259.64 262.99 rel. Err. 5.229% 3.334% 2.086% sec. 2.56 2.67 2.92 Table 4: Gaussian 2-factor model for n = 30 and Q = (0, 180) with computational times for the transition probabilities on a single CPU. 18 hal-00273790, version 2 - 15 Dec 2009 the diffusion and the pQWE-method perform similarly with respect to execution time and accuracy, whereas the spray method induces about double the error for half the execution time. In this 2-factor setting we also present results for a 365-day contract in table 5. N 100 250 500 diff 21943.77 22006.16 22042.31 rel. Err. 0.662% 0.379% 0.216% sec. 88 118 145 pQWE 21960.16 22001.90 22052.94 N 100 250 500 diff 6398.35 6470.05 6508.88 rel. Err. 2.084% 0.987% 0.393% sec. 88 118 145 pQWE 6380.18 6437.34 6483.65 K = 10 rel. Err. 0.588% 0.399% 0.168% K = 20 rel. Err. 2.363% 1.488% 0.779% sec. 58 62 67 MCSpray 21758.67 21580.61 21756.86 rel. Err. 1.500% 2.306% 1.508% sec. 22 24 27 sec. 58 62 67 MCSpray 6090.92 5926.04 6100.28 rel. Err. 6.789% 9.312% 6.646% sec. 22 24 27 Table 5: Gaussian 2-factor model for n = 365 and Q = (0, 2190) with computational times for the transition probabilities on a single CPU. Here we observe execution time in the range of about 1/2 minute up to 2 minutes for the single core implementation, which becomes quite critical for time sensitive applications. 6.3.1 Parallel implementation To demonstrate the parallel performance of the transition estimation methods, we first implemented the pQWE-methods using 365 threads (for each time layer one) on 8 Xeon cores. This approach reduces the execution to about 7 seconds (cf. Table 6), which means that it scales linearly with respect to number of available processors (and even a bit better, since these 365 threads exploit the computing power of a single core more efficiently than a single thread). Moreover, we developed a CUDA-implementation of the pQWE-method, which contains beneath the layer-wise parallelization also a MC-wise one, i.e. 365 blocks with each containing 256 threads for the MC-simulation per layer. This finally enabled us to price the 365-day contract at N = 100 in less than 1 second with about 2% accuracy using the 240 ScalarProcessors of an NVIDIA GTX 280 GPU. N 8 × Xeon GTX 280 100 7.12 0.75 250 7.33 1.26 500 7.79 2.11 Table 6: Computational time in seconds for pQWE-method on the 2-factor model with n = 365. Note, that the dual Xeon has an overall computation power of about 192 GFLOPS, whereas the GTX 280 can perform up to 933 GFLOPS. 6.3.2 Romberg extrapolation A further way to improve the accuracy of the quantization method is the extrapolation method from section 4.3 here applied to the 30-day contract from Table 4. 19 hal-00273790, version 2 - 15 Dec 2009 N 100 − 250 250 − 500 diff 2698.04 2698.09 rel. Err. 0.073% 0.071% N 100 − 250 250 − 500 diff 1798.25 1798.30 rel. Err. 0.109% 0.106% N 100 − 250 250 − 500 diff 922.31 922.43 rel. Err. 0.234% 0.220% N 100 − 250 250 − 500 diff 266.68 266.82 rel. Err. 0.716% 0.663% K=5 pQWE rel. Err. 2700.79 0.029% 2701.77 0.065% K = 10 pQWE rel. Err. 1801.01 0.045% 1801.97 0.098% K = 15 pQWE rel. Err. 924.79 0.035% 925.92 0.157% K = 20 pQWE rel. Err. 268.14 0.168% 268.75 0.060% MCSpray 2695.64 2697.09 rel. Err. 0.162% 0.108% MCSpray 1795.79 1797.24 rel. Err. 0.246% 0.165% MCSpray 918.73 921.23 rel. Err. 0.624% 0.351% MCSpray 263.03 266.34 rel. Err. 2.114% 0.844% Table 7: Extrapolation for Gaussian 2-factor model with n = 30 and Q = (0, 180) Here, we observe for all the methods a boost up of the accuracy, which in particular holds for the pQWE-approach. 6.4 NIG-model To finally apply the method in a non-Gaussian setting, we tested the exponential L´evy-model from section 5.3 for daily parameters α = 50, β = −2.0, δ = 0.02, µ = 0.001, s0 = 20 and a 30 days contract. Here again as in the 1-factor model (cf. Table 1), the dSpray-method performs very well and even outperforms in this case the diffusion method. Moreover it reveals a very stable extrapolation behavior, which can be seen in table 9. Since the computational times agreed with those from the 1-factor setting, we did not reproduce them again. N 50 100 200 diff 1820.22 1820.23 1820.24 N 50 100 200 diff 111.92 111.96 111.97 K = 10 rel. Err. dSpray 0.040% 1820.59 0.039% 1820.90 0.039% 1820.94 K = 20 rel. Err. dSpray 0.356% 111.93 0.321% 112.22 0.313% 112.30 rel. Err. 0.020% 0.003% < 0.001% rel. Err. 0.347% 0.088% 0.022% Table 8: NIG-model with n = 30 and Q = (0, 180) 20 N 50 − 100 100 − 200 diff 1820.23 1820.24 N 50 − 100 100 − 200 diff 111.97 111.97 K = 10 rel. Err. dSpray 0.039% 1821.00 0.039% 1820.95 K = 20 rel. Err. dSpray 0.311% 112.32 0.312% 112.32 rel. Err. 0.003% < 0.001% rel. Err. 0.002% < 0.001% Table 9: Extrapolation for NIG-model with n = 30 and Q = (0, 180) 7 Conclusions hal-00273790, version 2 - 15 Dec 2009 We have systematically carried out numerical tests on the pricing of swing options in the quantization tree framework using different methods for the transition probability estimation. The results are very promising, since we could show that • already in the 1-factor model the quantization approach outperforms the trinomial tree methods due to its better fit to the underlying distribution, • the computational time for the transition probability estimation in the 2-factor model with 365 exercise dates can be reduced to less than 1 second using our parallel approaches on a nowadays GPU device, • the quantization framework also works very well for non-Gaussian, in our case a NIG L´evy process, Markov Feller underlyings. Moreover, we could speed up the convergence of all these approaches by means of a Romberg extrapolation. References [1] V. Bally, G. Pag`es, and J. Printems. A quantization tree method for pricing and hedging multidimensional American options. Math. Finance, 15(1):119–168, 2005. [2] O. Bardou, S. Bouthemy, and G Pag`es. When are Swing options bang-bang and how to use it. LPMA preprint, 2007. [3] O. Bardou, S. Bouthemy, and G Pag`es. Optimal Quantization for the Pricing of Swing Options. Applied Mathematical Finance, 16(2):183–217, 2009. [4] O. E. Barndorff-Nielsen. Processes of normal inverse Gaussian type. Finance Stoch., 2(1):41– 68, 1998. [5] F.E. Benth, D. Frestad, and S. Koekebakker. Modeling the term structure dynamics in the Nordic Electricity swap market. Preprint, 2007. [6] F.E. Benth and J. Saltyte-Benth. The normal inverse Gaussian distribution and spot price modeling in energy markets. Intern. J. Theor. Appl. FinanceJournal, 7(2):177 – 192, 2004. [7] R. Cont and P. Tankov. Financial modelling with jump processes. Financial Mathematics Series. Chapman & Hall, 2004. 21 [8] J. H. Freidman, J. L. Bentley, and R. A. Finkel. An algorithm for finding best matches in logarithmic expected time. ACM Trans. Math. Softw., 3(3):209–226, 1977. [9] H.U. Gerber and E.S.W. Shiu. Option pricing by Esscher transforms. Trans. Soc. Actuaries, 46:99–191, 1994. [10] H.U. Gerber and E.S.W. Shiu. Actuarial Approach to Option Pricing. In Proceedings of the 5th AFIR International Colloquium, volume 1, pages 43–96, 1995. [11] A. Gersho and R.M. Gray. Vector Quantization and Signal Compression. Kluwer, Boston, 1992. [12] S. Graf and H. Luschgy. Foundations of Quantization for Probability Distributions. Lecture Notes in Mathematics n0 1730. Springer, Berlin, 2000. hal-00273790, version 2 - 15 Dec 2009 [13] A. Hull, J. White. Numerical procedures for implementing term structure models I: SingleFactor Models. Journal of Derivatives, 2(1):7–16, 1994. [14] P. Jaillet, E. I. Ronn, and S. Tompaidis. Valuation of Commodity-Based Swing Options. Managment Science, 50(7):909–921, 2004. [15] J. C. Kieffer. Uniqueness of locally optimal quantizer for log-concave density and convex error weighting function. IEEE Trans. Inform. Theory, 29(1):42–47, 1983. [16] G. Pag`es. A space vector quantization method for numerical integration. Journal of Applied and Computational Mathematics, 89:1–38, 1997. [17] G. Pag`es and J. Printems. Optimal quadratic quantization for numerics: the gaussian case. Monte Carlo Methods and Applications, 9(2):135–166, 2003. [18] G. Pag`es and J. Printems. Optimal Quantization for Finance: From Random Vectors to Stochastic Processes. In P.G. Ciarlet, editor, Handbook of Numerical Analysis, volume XV of Mathematical Modeling and Numerical Methods in Finance, pages 595–648. North-Holland, 2009. [19] D. Pollard. A central limit theorem for k-means clustering. Ann. Probab., 10(4):919–926, 1982. 22

© Copyright 2019