The two-dimensional tree–grid method

MA CM Bergische Universität Wuppertal Fachbereich Mathematik und Naturwissenschaften Institute of Mathematical Modelling, Analysis and Computational Mathematics (IMACM) Preprint BUW-IMACM 18/02 Igor Kossaczký, Matthias Ehrhardt and Michael Günther The Two-dimensional Tree-Grid Method February 2018 http://www.math.uni-wuppertal.de The Two-dimensional Tree-Grid Method Igor Kossaczkýa,∗, Matthias Ehrhardta , Michael Günthera a Lehrstuhl für Angewandte Mathematik und Numerische Analysis, Fakultät für Mathematik und Naturwissenschaften, Bergische Universität Wuppertal, Gaußstrasse 20, 42119 Wuppertal, Germany Abstract In this paper we introduce a novel two-dimensional Tree-Grid method for solving stochastic control problems with two space dimensions and one time dimension or equivalently, the corresponding Hamilton-Jacobi-Bellman equation. This new method can be regarded as a generalization of the Tree-Grid method for the stochastic control problems with one space dimension that was recently developed by the authors. This method overcomes two typical shortcomings: it is explicit and unconditionally stable. We prove the convergence of the method and exemplify it in the application on the two-factor uncertain volatility model. Keywords: 2D Tree-Grid Method, Hamilton-Jacobi-Bellman equation, Stochastic control problem, two-factor uncertain volatility model 1. Introduction Stochastic control problems (SCP) arise in many application where the stochastic (Itô) process is controlled with respect to time and state in order to optimize the expectation of some utility or cost function. The problem of 5 searching for an optimal control strategy can be often treated by solving the ∗ Corresponding author Email addresses: [email protected] (Igor Kossaczký), [email protected] (Matthias Ehrhardt), [email protected] (Michael Günther) URL: http://www-num.math.uni-wuppertal.de/en/amna/people/igor-kossaczky-msc.html (Igor Kossaczký), http://www-num.math.uni-wuppertal.de/~ehrhardt/ (Matthias Ehrhardt), http://www-num.math.uni-wuppertal.de/~guenther/ (Michael Günther) Preprint submitted to - February 23, 2018 underlying partial differential equation (PDE), the so-called Hamilton-JacobiBellman (HJB) equation, cf. [1]. 1.1. Overview of the numerical methods Either way, solving the SCP directly or solving the HJB equation, numerical 10 methods are often needed as closed form solution are rarely known. For the onedimensional SCP (with one controlled stochastic process for one state variable) various numerical methods were developed. Based on the approximation of the time-continuous stochastic process with a discrete chain, Markov chain approximation methods are presented in [1]. Implicit finite difference methods (FDM) 15 for solving HJB equations were presented for example in [2]. The advantage of these methods in contrast to the explicit FDMs is their unconditional stability and monotonicity. In [3], a method based on a transformation of the HJB equation was developed. The advantage of using this approach is that we don’t need to solve the optimization problem in each time-layer. 20 Recently a new explicit but still unconditionally stable and monotone scheme, the so-called Tree-Grid method, was developed by the authors [4]. This method combines the tree structure of the trinomial tree methods for option pricing, with the grid used typically in FDMs and can be regarded as some special explicit FDM or a Markov chain approximation. A useful modification of this method 25 was presented in [5]. As the explicitness and great flexibility of the method (unconditional stability with any grid spacing) exhibit clear advantages of this method, we are interested here in generalizing this method to two state-space dimensions (two controlled stochastic processes). In two space dimensions, a generalization of the implicit unconditionally 30 stable method of Wang and Forsyth [6] was presented in [7] and later used in [8]. The main idea of this method is combining the wide and the fixed stencil depending on the correlation in the particular time-space node. Alternatively, explicit methods based on the ideas of Kushner and Dupuis [1] presented in the papers [9, 10, 11] can be used. These are wide stencil schemes stable under some 35 CFL condition. Moreover, in these methods a linear interpolation of the grid 2 values is needed. Our generalization of the Tree-Grid method presented in this paper also falls into this class of wide stencil schemes, however it overcomes two shortcomings: it is unconditionally stable for any grid and no interpolation of the grid values is needed. Furthermore, as this two-dimensional (2D) Tree-Grid 40 method is explicit, it is also suitable for parallelization. 1.2. Problem formulation We are concerned with searching for the value function V (x, y, t) of the following 2-dimensional general stochastic control problem (SCP): ! Z Z T V (x, y, t) = max θ(x,y,t)∈Θ̄ k exp E r(∗l )dl f (∗k )dk t + exp Z t ! T ! r(∗k )dk VT (XT , YT ) Xt = x, Yt = y , t (1) dXt =µx (∗t )dt + σx (∗t )dWtx , (2) dYt =µy (∗t )dt + σy (∗t )dWty , (3) dWtx dWty =σxy (∗t )/(σx (∗t )σy (∗t )) (4) with ∗t =(Xt , Yt , t, θ(Xt , Yt , t)), 0 < t < T, x, y ∈ R, where x, y are the state variables and t is time. Here, Θ̄ denotes the space of all suitable control functions (see e.g. [1, 12]) from R2 × [0, T ] to a set Θ. For our purpose, we will assume Θ to be discrete. If this is not the case, we can easily achieve this property by its discretization. Now following the Bellman’s principle, the dynamic programming equation holds: V (x, y, tj ) = max θ(x,y,t)∈Θ̄tj + exp Z tj+1 E Z tj+1 exp tj ! Z k ! r(∗l )dl f (∗k )dk tj ! r(∗k )dk V (Xtj+1 , Ytj+1 , tj+1 ) Xtj = x, Ytj = y , tj (5) where 0 ≤ tj < tj+1 ≤ T are some time-points and Θ̄tj is a set of control functions from Θ̄ restricted to the R2 × [tj , tj+1 ) domain. Using the dynamic programming equation (5), it can be shown [12], that solving the SCP (1)–(4) 3 is equivalent to solving the two-dimensional Hamilton-Jacobi-Bellman (HJB) equation:
∂V + max LV + r(·)V + f (·) = 0, θ∈Θ ∂t ∂2V ∂V σy (·)2 ∂ 2 V ∂V σx (·)2 ∂ 2 V + σxy (·) + µx (·) + + µy (·) LV = 2 2 ∂x ∂x∂y 2 ∂y 2 ∂x ∂y V (x, y, T ) = VT (x, y), 0 < t < T, (6) (7) (8) x, y ∈ R, where σ(·), µ(·), r(·), f (·) are functions of x, y, t, θ. We should note that the maximum operator in (1) and (6) can be replaced by a minimum operator and the whole following analysis will hold analogously. A possible generalization 45 of the SCP with corresponding HJB equation is a stochastic differential game with the corresponding Hamilton-Jacobi-Bellman-Isaac equation [13]. Further generalizations can be found in [14]. 2. Construction of the 2D Tree-Grid algorithm In this section we will derive the two-dimensional Tree-Grid algorithm for 50 solving the problem (1)–(4). We will use the ideas that were widely explained in [4], [5] and we refer interested readers to these papers. We will work on a three-dimensional rectangular domain with two space dimensions and one time dimension. For a fixed control θ, the candidate for a value in each node V (xi , yj , tk ) will be computed from seven values from the next time layer tk+1 . 55 Figure 1 illustrates this approach. We will denote these seven values in this context simply as stencil located at (xi , yj , tk ). The weights of these seven values can be interpreted as probabilities and therefore we demand, that the moments of such discrete random variable are matching with the moments of the increment of the two-dimensional stochastic process (2)–(4) with the fixed 60 control θ at least asymptotically. Then, the actual value V (xi , yj , tk ) will be computed as a maximum of the candidates. For handling nodes close to the boundary we suppose that we know, how the solution behaves in the outer 4 neighbourhood of the boundaries, and that we can describe this behaviour with a boundary function BC(x, y, t). The terminal condition in the time-layer tM = T 65 reads V (x, y, tM ) = VT (x, y). Figure 1: Illustration of the two-dimensional Tree-Grid structure. Only the red nodes in the later time layer impact the blue node. On the other hand, the red nodes can be interpreted as possible future states if we are in the blue node state. The figure illustrates a situation with positive correlation and a variance that is larger in the x-direction than in the y-direction. The stencil size in each direction is roughly proportional to the square root of the variance coefficient in that direction multiplied by the discretization parameter h. 2.1. Notation At first, before discussing how to choose the stencil nodes around node (xi , yj , tk ), and the proper weights, we present here the notation used in the sequel: 70 • x1 , x2 , . . . , xNx -space discretization in the first space dimension. y1 , y2 , . . . , yNy -space discretization in the second space dimension. t1 , t2 , . . . , tM -time discretization . • ∆k t = tk+1 − tk -(current) time-step. We will use equidistant time stepping, ∆k t = ∆t for all k = 1, 2, . . . , M − 1. A generalization to non- 75 equidistant time stepping is straightforward. • ∆i x = xi+1 − xi , ∆j y = yj+1 − yj space-steps in the first and second dimension (possibly non-equidistant). • ∆x = maxi ∆i x, ∆y = maxi ∆i y. 5 • h = max(K max(∆x, ∆y), ∆t), where K > 0 is a parameter used for 80 regulating the stencil size. In the non-equidistant case, ∆t should be replaced by ∆k t. • b = h/∆t. In the non-equidistant case, this should be replaced by b = h/∆k t in the following algorithm. • (xi , yj , tk ) -the node for which the stencil is constructed. 85 • E∗ = µ∗ (xi , yj , tk , θ)∆t for ∗ = x, y.
• V ar∗ = σ∗2 (xi , yj , tk , θ) + O(h) ∆t for ∗ = x, y, will be determined later. • ρ(xi , yj , tk , θ) = σxy (xi , yj , tk , θ)/(σx (xi , yj , tk , θ)σy (xi , yj , tk , θ)) √  h , will be determined later. • σ̃xy (xi , yj , tk , θ) = σxy (xi , yj , tk , θ) + O  √  • Cov = σxy (xi , yj , tk , θ) + O h ∆t, will be determined later. 90 • W∗ = V ar∗ + E∗2 for ∗ = x, y. • Wxy = Cov + Ex Ey . k -numerical approximation of V (xi , yj , tk ) • vi,j 2.2. Choosing the stencil nodes Next, we describe how to choose the stencil nodes around an arbitrary node 95 (xi , yj , tk ) for a fixed control θ. The values in these nodes will impact the value in the node (xi , yj , tk ). If Ex > 0, j k p x− = xi − 2Wx b , x else if Ex < 0, l m p x+ = xi + 2Wx b , x else (Ex = 0), l  m p x+ = max xi + 2Wx b, xi + (xi − x− ) , (9) x j  k p x− = min xi − 2Wx b, xi − (x+ − xi ) , (10) x j k p x− = xi − 2Wx b , m l p x+ = xi + 2Wx b , x x 6 (11) where ⌈⌉x resp. ⌊⌋x denotes rounding to the nearest greater resp. smaller element from x1 , x2 , . . . , xNx . If such element does not exist, ⌈z⌉x resp. ⌊z⌋x will return just z. If Ey > 0, j k p y− = yj − 2Wy b , l  m p y+ = max yj + 2Wy b, yj + (yj − y− ) , (12) l m p y+ = yj + 2Wy b , j  k p y− = min yj − 2Wy b, yj − (y+ − yj ) , (13) y else if Ey < 0, y else (Ey = 0), y y l m p y+ = yj + 2Wy b , j k p y− = yj − 2Wy b , y 100 y where ⌈⌉y , ⌊⌋y are defined analogously. (14) The following nodes with the respective weights (probabilities) will be used in the stencil located at (xi , yj , tk ): • node (xi , yj , tk+1 ) with the probability po , • nodes (x+ , yj , tk+1 ) and (x− , yj , tk+1 ) with the probabilities px+ and px− , 105 • nodes (xi , y+ , tk+1 ) and (xi , y− , tk+1 ) with the probabilities py+ and py− , • nodes (x+ , y+ , tk+1 ) and (x− , y− , tk+1 ), both with the probability pxy if Wxy is non-negative, and nodes (x+ , y− , tk+1 ) and (x− , y+ , tk+1 ) both with the probability pxy if Wxy is negative. In the following algorithm we will focus on the case of non-negative Wxy , the case of a negative Wxy is 110 treated analogously. Moreover, we define the difference operators ∆+ x = x+ − xi , ∆− x = x i − x − , (15) ∆ + y = y+ − yi , ∆− y = yi − y− . (16) 7 2.3. Choosing the stencil weights (probabilities) To match the first two moments of the approximative increment of the stochastic process and of the increment of the discrete process defined by the “stencil nodes” and their probabilities and to ensure that the probabilities are positive and sum up to 1, we demand: po + px+ + px− + py+ + py− + 2pxy = 1, (17) (px+ + pxy)∆+ x − (px− + pxy)∆− x = Ex , (18) (py+ + pxy)∆+ y − (py− + pxy)∆− y = Ey , (19) (px+ + pxy)(∆+ x)2 + (px− + pxy)(∆− x)2 = Wx , (20) (py+ + pxy)(∆+ y)2 + (py− + pxy)(∆− y)2 = Wy , (21) pxy (∆+ x∆+ y + ∆− x∆− y) = Wxy . (22) po , px+ , px− , py+ , py− , pxy ≥ 0 (23) For negative Wxy , only the condition (22) changes to pxy (∆+ x∆− y + ∆− x∆+ y) = |Wxy |. (24) The solution of the six equations (17)–(22) reads ∆+ x∆− x∆+ y∆− y − ∆+ x∆− xWy − ∆+ y∆− yWx , ∆+ x∆− x∆+ y∆− y Ey (∆+ y − ∆− y) Ex (∆+ x − ∆− x) 2|Wxy | , + + + ∆+ y∆− y ∆+ x∆− x ∆c |Wxy | W x + Ex ∆ − x , − = 2 (∆+ x) + ∆+ x∆− x ∆c W x − Ex ∆ + x |Wxy | = , − 2 (∆− x) + ∆+ x∆− x ∆c |Wxy | W y + Ey ∆ − y , − = 2 (∆+ y) + ∆+ y∆− y ∆c |Wxy | W y − Ey ∆ + y , − = 2 (∆− y) + ∆+ y∆− y ∆c |Wxy | = , ∆c po = px+ px− py+ py− pxy (25) (26) (27) (28) (29) (30) where ∆c = ∆+ x∆+ y + ∆− x∆− y for non-negative (σ̃xy (xi , yj , tk , θ)∆t + Ex Ey ) and ∆c = ∆+ x∆− y + ∆− x∆+ y for negative (σ̃xy (xi , yj , tk , θ)∆t + Ex Ey ). 8 Following the construction of x+ , x− , y+ , y− , it is easy to check that po is 115 always non-negative. The same holds for pxy . To ensure also the non-negativity of px+ , px− , py+ and py− we have to properly choose the variables V arx , V ary , and Cov to get non-negative weights, while still remaining consistent with the original problem as h → 0. This is done in the next Section 2.4. 2.4. Artificial diffusion and covariance adjustment Let us assume Ex ≥ 0. Now, the first fraction of px+ is positive, however the first fraction of px− may be negative. We will set V arx (and hence Wx ) in such way, that it will be also positive. It holds p 2b(V arx + Ex2 ) + 2∆x). (31)   p Ax = 1/2 |Ex | 4b2 Ex2 + 16b∆x|Ex | − (2b − 2)Ex2 + 4∆x|Ex | . (32) Wx − Ex ∆+ x > V arx + Ex2 − Ex ( The right-hand side of the inequality (31) is 0 for V arx = Ax and greater than 0 for V arx > Ax with We replaced here Ex with |Ex | to cover also the analogous case Ex < 0 and possibly negative px+ . Now, if we set V arx = max(σx2 (xi , yj , tk , θ)∆t, Ax , Ex2 ), (33) V ary = max(σy2 (xi , yj , tk , θ)∆t, Ay , Ey2 ), (34) where Ay is defined analogously, the first fraction in px+ , px− , py+ and py− will be non-negative. The possible difference between V arx resp. V ary and the variances from the original problem σx2 ∆t resp. σy2 ∆t is called artificial diffusion. Now, taking into account the correlation coefficient in the current node, following these definitions of variances V arx , V ary we introduce also the new covariance σ̃xy p V arx V ary =ρ ∆t 9 (35) and Cxy := min  W x − Ex ∆ + x W x + Ex ∆ − x ∆c , ∆c , (∆+ x)2 + ∆+ x∆− x (∆− x)2 + ∆+ x∆− x W y − Ey ∆ + y W y + Ey ∆ − y ∆c , ∆c , (∆+ y)2 + ∆+ y∆− y (∆− y)2 + ∆+ y∆− y  |σ̃xy (xi , yj , tk , θ)∆t + Ex Ey | . Now, we define Cov as   Cxy − Ex Ey Cov =  −Cxy − Ex Ey if (σ̃xy (xi , yj , tk , θ)∆t + Ex Ey ) ≥ 0, (36) (37) if (σ̃xy (xi , yj , tk , θ)∆t + Ex Ey ) < 0. This covariance Cov is consistent with the variances V arx , V ary defined by (33), (34) as for σ̃xy (xi , yj , tk , θ)∆t + Ex Ey ≥ 0 it holds − p p V arx V ary ≤ −Ex Ey ≤ Cov ≤ σ̃xy (xi , yj , tk , θ)∆t ≤ V arx V ary , (38) p p V arx V ary ≤ σ̃xy (xi , yj , tk , θ)∆t ≤ Cov ≤ −Ex Ey ≤ V arx V ary . (39) and for σ̃xy (xi , yj , tk , θ)∆t + Ex Ey < 0 it holds − Now it also holds |Wxy | = Cxy , which implies that px+ , px− , py+ and py− are all positive. It is easy to check that it holds V arx = σx2 (xi , yj , tk , θ) + O(h), ∆t V ary = σy2 (xi , yj , tk , θ) + O(h), ∆t (40) and V arx = σx2 (xi , yj , tk , θ), ∆t resp. V ary = σy2 (xi , yj , tk , θ), ∆t (41) for σx2 (xi , yj , tk , θ) 6= 0 resp. σy2 (xi , yj , tk , θ) 6= 0 and h small enough. It holds  √  Wx + Ex ∆− x = σx2 + O h ∆t, (42) (∆+ x)2 + ∆+ x∆− x = 4hσx2 + O(h3/2 ), ∆c = 4hσx σy + O(h3/2 ). (43) (44) It follows that  √  W x + Ex ∆− x h ∆t ∆c = σx σ y + O 2 (∆+ x) + ∆+ x∆− x 10 (45) and the same holds also for the second, third and fourth maximum candidate in (36). Moreover,  √  |σ̃xy ∆t + Ex Ey | = |σ̃xy |∆t + O((∆t)2 ) = |σxy | + O h ∆t. (46) Therefore,  √  h ∆t, Cxy = |σxy | + O ⇒ √  Cov h = σxy + O ∆t (47) and it is easy to check that Cov = σxy ∆t 120 (48) for σx2 6= 0, σy2 6= 0, |σxy | 6= σx σy and h small enough. Following (40),(47), the modified variances and the modified covariance are consistent with the variances and the covariance from the original problem. Moreover, for σx , σy 6= 0 and |σxy | < σx σy , the modified variances and the covariance will be equal to the original ones for h small enough. 125 2.5. Setting the parameter K and the stencil size reduction In the formula for h = max(K max(∆x, ∆y), ∆t), the part K max(∆x, ∆y) is responsible for the consistency of the correlation in the numerical model with the correlation of the original problem. Here, the parameter K > 0 can be chosen arbitrarily, however for a large (relative to the problem parameters) K, 130 the stencil is large, what typically increases the error. On the other hand, for too small K, the correlation may start being exact (or exact enough) only for very fine grids, while not being sufficiently exact on the coarse grids, resulting in larger errors on these coarse grids. For K = 0 we can’t guarantee the convergence of the correlation. However, the correlation in the numerical model 135 may match with the correlation from the original problem even for small K or K = 0. This motivates the following multiple K modification of the Tree-Grid method; For each control in each node: 1. set l = 1 and use K = K0 , K0 ≥ 0 to compute x∗ , y∗ , p∗ 11 140 2. compute the correlation of the numerical model: ρ̃ = Cov/ p V arx V ary 3. if ρ̃ 6= ρ: recompute x∗ , y∗ , p∗ with K = Kl , Kl > Kl−1 and set l := l + 1 else break. 4. if l < lmax : go to step 2, else break. 145 Using this modification, we will use a smaller stencil size as much as possible. This approach can be seen as some analogy to the approach of Ma and Forsyth [7], where a fixed (and thus small) stencil is used as much as possible. However, here we will not increase the convergence rate, but possibly reduce the error. Another approach, the non-constant K modification, is to use a non-constant 150 K = K(x, y, θ). This K can be smaller in nodes with large volatilities σx , σy and larger in nodes with smaller σx , σy , regulating the stencil size to not explode in case of large volatilities. Both modifications can be also combined and it is easy to check, that they do not harm the consistency. In our numerical simulation these modifications 155 didn’t lead to a significant improvement and therefore we used simply a constant K = 1/400. However, a non-constant K modification may be useful for other stochastic control problems. 3. The 2D Tree-Grid Algorithm Finally, we can use the stencil nodes and weights to construct the 2D Tree160 Grid algorithm. We define the function v k+1 : If (x, y) ∈ {x1 , x2 , . . . , xNx } × {y1 , y2 , . . . , yNy } : k+1 v k+1 (x, y) = vi,j so that (x, y) = (xi , yj ), else: v k+1 (x, y) = BC(x, y, tk+1 ). 12 (49) For a given space-node (xi , yj ) and a given control θ we define vok+1 = v k+1 (xi , yj ), k+1 vx+ = v k+1 (x+ , yj ), k+1 vy+ = v k+1 (xi , y+ ), k+1 vy− = v k+1 (xi , y− ), k+1 vx− = v k+1 (x− , yj ), (50) (51) If Wxy ≥ 0 : k+1 vxy+ = v k+1 (x+ , y+ ), k+1 vxy− = v k+1 (x− , y− ). (52) If Wxy < 0 : k+1 vxy+ = v k+1 (x+ , y− ), k+1 vxy− = v k+1 (x− , y+ ). (53) k fi,j (θ) = f (xi , yj , tk , θ), k ri,j (θ) = r(xi , yj , tk , θ). (54) Now, the discretized version of the dynamic programming equation for i = 2, 3, . . . , Nx − 1, j = 2, 3, . . . , Ny − 1, k = 1, 2, . . . , M − 1 reads k k = max wi,j (θ) vi,j (55) θ∈Θ k k k (θ) =fi,j (θ)∆t + (1 + ri,j (θ)∆t) wi,j  k+1 k+1 k+1 k+1 + py+ vy+ + py− vy− + px− vx− · px+ vx+  k+1 k+1 + po vok+1 + pxy (vxy+ + vxy− ) . (56) For the boundary nodes (xi , yj ) (i ∈ {1, Nx } or j ∈ {1, Ny }) we employ the boundary condition: k vi,j = BC(xi , yj , tk ). (57) The terminal condition is defined as M vi,j = VT (xi , yj ). (58) Finally we can summarize the whole algorithm of the 2D Tree-Grid Method: 13 Algorithm 1 The 2D Tree-Grid method M = VT (xi , yj ) for i = 1, 2, . . . , Nx , j = 1, 2, . . . , Ny ; 1: Set vi,j 2: 3: 4: 5: for k = M − 1, M − 2, . . . , 1 do for i = 1, 2, 3, . . . , Nx do for j = 1, 2, 3, . . . , Ny do if i ∈ {1, Nx } or j ∈ {1, Ny } then k Compute vi,j according to (57); 6: 7: else for θ ∈ Θ do 8: Determine Ex , Ey according to Section 2.1; 9: 10: Compute Ax , Ay according to (32); 11: Compute V arx , V ary according to (33), (34); 12: Determine Wx , Wy according to Section 2.1; 13: Determine x+ , x− , y+ , y− according to (9)–(14); 14: Determine ∆+ x, ∆− x, ∆+ y, ∆− y according to (15)–(16); 15: Compute σ̃xy according to (35); 16: Compute Cxy according to (36); 17: Compute Cov according to (37); 18: Determine Wxy according to Section 2.1; 19: Compute po , px+ , px− , py+ , py− , pxy according to (25)–(30); 20: Using (49), compute all variables defined by (50)–(54); 21: k Compute wi,j (θ) according to (56); 22: end for 23: k Compute vi,j according to (55); 24: 25: 26: 27: end if end for end for end for 14 4. Convergence In this section we will prove the convergence of the 2D Tree-Grid method. First we will quickly summarize the classical convergence theory of Barles and 165 Souganidis [15] and in the second part of this section we will use this theory to prove the convergence of our scheme. Let us note that the algorithm was derived by discretizing the dynamics in the original SCP (1)–(4), but we will prove that the scheme is consistent with the HJB equation (6)–(8). 4.1. The convergence theory Let U denote some suitable function space. Let us define some nonlinear differential operator F F : U → R, V (s) → F V (s). We suppose there exists a viscosity solution (see [16]) of the equation F V (s) = 0, and denote this solution simply by V (s). To find some approximation of the viscosity solution we define a discrete approximation scheme
Gv(s) = G v(s), v(s + b1 h), v(s + b2 h), . . . , v(s + bn h) , 170 (59) where v(s), s ∈ RK is defined as a (possibly) multidimensional function, bi ∈ RK , i = 1, 2, . . . , n and h ∈ R+ . Let us consider the system of sets called discretized domains Sh = {si ∈ RK |i = 1, 2, . . . , Nh }, defined for different values of h, which is often referred as step size. Definition 1 (Numerical scheme). The system of equations Gv(s) = 0 with s ∈ Sh depending on a parameter h ∈ R+ is called a numerical scheme. 175 The numerical scheme is well-defined, if it possess an unique solution. We will assume that this condition is met for any feasible h. By v(s), we will denote an approximation of the solution of F V (s) = 0, computed by solving the system of equations Gv(s) = 0, s ∈ Sh . In order to distinguish between approximations with different h, we will sometimes denote v(s) as vh (s). 15 180 Definition 2 (Monotonicity). A discrete approximation scheme
Gv(s) = G v(s), v(s + b1 h), v(s + b2 h), . . . , v(s + bn h) is monotone, if the function G is non-increasing in v(s + bi h) for bi 6= 0, i = 1, . . . , n. Definition 3 (Consistency). The scheme Gv(s) = G(v(s), v(s + b1 h), v(s + b2 h), . . . , v(s+bn h)) is a consistent approximation of F V (s), if limh→0 |F φ(s)− 185 Gφ(s)| = 0, for any C ∞ -smooth test function φ(s). A scheme is consistent on a numerical domain, if it is consistent in all points of this numerical domain. In such case we will call the scheme consistent. In literature, often C 2 -smooth test functions are used. However, as shown for example in [17], this leads to an equivalent definition. 190 Definition 4 (Stability). The numerical scheme defined by the system of equations Gvh (s) = 0, s ∈ Sh with solution vh (s) is stable, if there exists some constant C so that kvh (s)k∞ < C, ∀h > 0. The following Theorem of Barles and Souganidis [15] is the key for proving convergence of a numerical scheme approximating a nonlinear PDE: 195 Theorem 1 (Barles-Souganidis). If the equation F V (s) = 0 satisfies the strong uniqueness property (see [15]) and if the numerical scheme Gvh (s) = 0, s ∈ Sh approximating the equation F V (s) = 0 is monotone, consistent and stable, its solution vh (s) converges locally uniformly to the solution V (s) of F V (s) = 0 with h → 0. 200 The above mentioned strong uniqueness property [15] is a property of the problem and not of the numerical scheme. Therefore, we will simply assume that our problem possess this property without actually proving it. 4.2. Convergence of the 2D Tree-Grid method In this section, we will prove the convergence of the 2D Tree-Grid method. For the purpose of this convergence analysis we will rewrite equations (55), (56) 16 as k+1 k+1 k+1 k+1 k+1 k+1 Gv(xi , yj , tk ) =G(vok+1 , vx+ , vx− , vy+ , vy− , vxy+ , vxy− )  
1 k k k vi,j − max fi,j (θ)∆t + 1 + ri,j (θ)∆t = θ∈Θ ∆t  k+1 k+1 k+1 k+1 · px+ vx+ + px− vx− + py+ vy+ + py− vy−  k+1 k+1 + po vok+1 + pxy (vxy+ + vxy− ) = 0. (60) Using the theory of the previous Section 4.1, our goal is to show that equation (60) is a monotone, consistent, and stable approximation of the nonlinear differential operator F V defined by the PDE (6): F V (x, y, t) = −   ∂V − max LV + r(·)V + f (·) . θ∈Θ ∂t (61) Let us define the difference operators ∆o+ z = ∆+ z, ∆o− z = −∆− z, (62) for z ∈ {x, y}. At first, we will show the consistency of the scheme in an 205 arbitrary point (xi , yj , tk ): Lemma 1 (Consistency). The discrete scheme (60) is consistent with the PDE operator (61). Proof. Let φ : R2 × [0, T ] → R be a C ∞ -smooth function. Let us define φk (x, y) = φ(x, y, tk ) and use the short notation defined by (50)–(54) for φk+1 instead of v k+1 . At first, let us sketch the main idea of the proof of consistency in an arbitrary point (xi , yj , tk ): We will write all values φk+1 for α α ∈ {x+, x−, y+, y−, xy+, xy−} as Taylor expansions around φki,j . We will substitute these Taylor expansions into the discrete scheme (60), group terms with the same derivatives together and estimate the sum of the coefficients in front of them. We will end up with the PDE operator (61) and some terms of order O(hλ ), where λ = 1/2 if |ρ| = 1 or σx = 0 or σy = 0 and λ = 1 else. Let us suppose that Wxy ≥ 0 in (xi , yj , tk ). The case of negative Wxy can be 17 treated analogously. For ∗ ∈ {+, −}, let us define the operators: ∂φ ∂φ 1 ∂2φ ∂2φ 2 ∆o∗ x + ∆o∗ y + (∆ x) + ∆o∗ x∆o∗ y o∗ ∂x ∂y 2 ∂x2 ∂x∂y 1 ∂2φ (∆o∗ y)2 , (63) + 2 ∂y 2 1 ∂3φ 1 ∂3φ 1 ∂3φ 3 2 Bxy∗ φ = (∆ x) + (∆ x) ∆ y + ∆o∗ x(∆o∗ y)2 o∗ o∗ o∗ 6 ∂x3 2 ∂x2 ∂y 2 ∂x∂y 2 1 ∂3φ + (∆o∗ y)3 , (64) 6 ∂y 3 1 ∂4φ 1 ∂4φ 1 ∂4φ 4 3 (∆ x) + (∆ x) ∆ y + (∆o∗ x)2 , (∆o∗ y)2 Cxy∗ φ = o∗ o∗ o∗ 24 ∂x4 6 ∂x3 ∂y 4 ∂x2 ∂y 2 1 ∂4φ 1 ∂4φ + ∆o∗ x(∆o∗ y)3 + (∆o∗ y)4 . (65) 3 6 ∂x∂y 24 ∂y 4 Axy∗ φ = k Now we can expand φk+1 xy∗ around φi,j as follows: k φk+1 xy∗ =φi,j +
∂ ∂ k Axy∗ φki,j ∆t + Bxy∗ φki,j φi,j ∆t + Axy∗ φki,j + ∂t ∂t ∂ (Bxy∗ Rxy∗ ) ∆t + Cxy∗ Rxy∗ + O((∆t)2 ), ∂t (66) Rxy∗ = φ(xi + ǫxy∗ ∆o∗ x, yj + δxy∗ ∆o∗ y, tk + γxy∗ ∆t), (67) + where k+1 for some ǫxy∗ , δxy∗ , γxy∗ ∈ [0, 1]. We expand φk+1 x∗ , φy∗ , ∗ ∈ {+, −} in the k+1 same manner: for φk+1 x∗ we only need to substitute ∆o∗ y with 0 and for φy∗ we only need to substitute ∆o∗ x with 0 in the Taylor expansion (66) and change the index xy∗ to x∗ resp. y∗ in all expressions (63)–(67). Now, by substituting the Taylor expansions into the scheme Gφki,j defined by (60) we get: Gφki,j = 
1  k k k φi,j − max fi,j (θ)∆t + 1 + ri,j (θ)∆t θ∈Θ ∆t k   ∂φ i,j ∆t + (Lφki,j + O(hλ ))∆t · φki,j + ∂t = F φki,j + O(hλ ), (68) where λ = 1/2 if |ρ| = 1 or σx = 0 or σy = 0 and λ = 1 else. In the first equation of (68) we used the estimates of the linear combinations of higher order terms 210 from the Taylor expansions. The coefficients of these linear combinations are 18 the probabilities po , px+ , px− , py+ , py− , pxy+ , pxy− (according to the definition of G (60)). We describe here these estimates: 1. For the linear combinations of the terms included in Aα φki,j , α ∈ {x+, x−, y+, y−, xy+, xy−} we used the properties of the scheme 215 (17)–(22). After summing all expressions obtained by applying (17)–
(22) we get Lφki,j + O(h1/2 ) ∆t. Moreover, following Section 3, if |ρ| = 6 1, σx 6= 0, σy 6= 0, for h small enough we end up just with Lφki,j ∆t. 2. In the same way, for the linear combinations of the terms included in


∂ λ 2 k k ∂t Aα φi,j ∆t we get L ∂φi,j /∂t + O(h ) (∆t) = O(h)∆t. 3. For the linear combinations of the terms included in Bα φki,j we used the following estimates: (px+ + pxy+ )(∆o+ x)3 + (px− + pxy+ )(∆o− x)3 = Wx (∆o+ x + ∆o− x) − Ex ∆o+ x∆o− x = O(h)∆t, (69) pxy+ (∆o+ x)2 ∆o+ y + pxy− (∆o− x)2 ∆o− y =
Wxy (∆o+ x)2 ∆o+ y + (∆o− x)2 ∆o− y = O(h)∆t. ∆o+ x∆o+ y + ∆o− x∆o− y 220 (70) Here we used that ∆o+ x+∆o− x = O(h), (∆o+ x)2 ∆o+ y+(∆o− x)2 ∆o− y = O(h2 ). We used analogous estimates also for the terms where x and y are switched symmetrically. 4. As (∆o∗ x)2 = O(h), (∆o∗ y)2 = O(h) it is clear that all terms in ∂ ∂t (Bα Rα ) ∆t are of order O(h)∆t. 5. For the linear combinations of the terms included in Cα Rα we constructed the estimate in the following way: Let us define ψ = (px+ ax+ + pxy+ axy+ )(∆o+ x)2 + (px− ax− + pxy− axy− )(∆o− x)2 . Then, it holds:

m σx2 + O(h) ∆t ≤ ψ ≤ M σx2 + O(h) ∆t, m = min (ax+ , ax− , axy+ , axy− ) , 19 M = max(ax+ , ax− , axy+ , axy− ). 225 If m = O(h), M = O(h) then ψ = O(h)∆t. Now for (ax+ , ax− , axy+ , axy− ) equal to • 1 ∂4 24 ∂x4 • 1 ∂4 6 ∂x3 ∂y • 230 4 1 ∂ 4 ∂x3 ∂y
(∆o+ x)2 Rx+ , (∆o− x)2 Rx− , (∆o+ x)2 Rxy+ , (∆o− x)2 Rxy− , (0, 0, ∆o+ x∆o+ yRxy+ , ∆o− x∆o− yRxy− ),
0, 0, (∆o+ y)2 Rxy+ , (∆o− y)2 Rxy− , which are all O(h), we get that the linear combinations of the first three terms in Cα Rα are of order O(h)∆t. Using analogous estimates, also the linear combinations of the fourth and fifth term are of order O(h)∆t. Using the estimates above, we proved (68) which is the first order consistency of our scheme if |ρ| < 1, σx > 0, σy > 0 and consistency of order 1/2 else. 235 Now the lemma establishing monotonicity of the scheme follows: k Lemma 2 (Monotonicity). If 1 + ri,j (θ)∆t ≥ 0 for all possible i, j, k, θ, then the discrete scheme (60) is monotone. Proof. In this case, monotonicity is a direct implication of the non-negativity of po , px+ , px− , py+ , py− , pxy+ , pxy− . 240 k Remark 1. As already mentioned by the authors in [4], even if 1+ri,j (θ)∆t < 0 for some i, j, k, θ, we can obtain a monotone scheme if we substitute 1 + k ri,j (θ)∆t by k 1/(1 − ri,j (θ)∆t) in (60) for these parameters i, j, k, θ. Note that this change does not harm the consistency, nor the stability of the scheme. 245 The stability analysis of the 2D Tree-Grid method is basically identical to the stability analysis of the one-dimensional Tree-Grid method done in [4]. Therefore we state here just the stability condition and the Lemma about stability of the scheme. Property 1 (Stability condition of the problem). We suppose that: 250 1. There exist constants Cf , Cr such that for all x, y, t, θ ∈ [x1 , xNx ]×[y1 , yNy ]× [t1 , tM ] × Θ holds: |f (x, y, t, θ)| < Cf , |r(x, y, t, θ)| < Cr . 20 2. There exist a constant CBC such that |BC(x, y, t)| < CBC holds for all t ∈ [t1 , tM ] and for all possible outer-domain values (x, y) of the variables (x+ , yj ), (x− , yj ), (xi , y+ ), (xi , y− ), (x+ , y+ ), (x+ , y− ), (x− , y+ ), (x− , y− ) for 255 any grid. The lemma about stability of the scheme follows: Lemma 3 (Stability). If the problem satisfies the stability conditions defined in Property 1, then the scheme (60) is stable. Proof. The proof of this lemma can be found in [4]. 260 Finally, the convergence theorem follows: Theorem 2 (Convergence of the 2D Tree-Grid Method). The approximation computed with the 2D Tree-Grid Method (defined by Algorithm 1) for solving the SCP (1)–(4) and the corresponding HJB equation (6)–(8) satisfying the strong uniqueness property (see [15]) and the stability conditions defined in Property 1 265 converges to the viscosity solution of this SCP (and the HJB equation). Proof. The proof follows from the Theorem 1 and Lemmas 1, 2, 3. 5. Application to option pricing models 5.1. Two-factor uncertain volatility model Here we will use the 2D Tree-Grid method for pricing options on two different risky assets under uncertain volatility. In this setting, the volatilities and the correlation of the assets are only known to lie in certain bounds. In this case the maximal option price can be computed as the solution of the HJB equation ∂V + max (LV − rV ) = 0, θ∈Θ ∂t σy2 2 ∂ 2 V σ2 ∂ 2 V ∂2V ∂V ∂V LV = x x2 2 + ρσx σy xy + y + rx + ry , 2 ∂x ∂x∂y 2 ∂y 2 ∂x ∂y V (x, y, T ) = VT (x, y), 0 < t < T, (71) (72) (73) x, y ∈ R. 21 Here, x, y, t denote the first and the second asset price and the time, r denotes the risk-free interest rate, T is the maturity time, the terminal condition VT (x, y) is the payoff function and for the control variable holds θ = (σx , σy , ρ), where σx , σy , ρ are uncertain volatilities of the first and the second asset, and their correlation. For the set of possible controls holds: Θ = [σx,min , σx,max ] × [σy,min , σy,max ] × [ρmin , ρmax ]. (74) To obtain the minimal option price, we have to replace ’max’ by ’min’ in the HJB equation (71). This model was discussed for example in [18] and later solved with a hybrid (combining fixed and wide stencils) implicit method in [7]. As explained in [7], the optimal control lies in a subset of Θ, Θ̃ =  ({σx,min , σx,max } × [σy,min , σy,max ])  ∪ ([σx,min , σx,max ] × {σy,min , σy,max }) × {ρmin , ρmax }. (75) Therefore, we will search for the control in Θ̄, a discretized version of Θ̃. To 270 verify the implementation, we will also solve a problem with an one-element set Θ̄ = {(σx , σy , ρ)}. In this case, the equation (71) is simply the 2-dimensional Black-Scholes equation [19] for which the analytical solution is known. Terminal condition: As a terminal condition, we use the payoff function in the form of a butterfly spread on the maximum of two assets: + VT (x, y) = (M − K1 )+ − 2 (M − (K1 + K2 )/2) + (M − K2 )+ M = max(x, y), K1 = 34, K2 = 46. Boundary conditions: We will use Dirichlet boundary conditions. On the upper and the right boundary, we will set the value to 0: BC(x > xmax , y, t) = 0, BC(x, y > ymax , t) = 0, xmax = ymax = 144. To verify that our computational domain is large enough, we solved the HJB equation and the Black-Scholes equation also for xmax = ymax = 400 and obtained in node (x, y, t) = (40, 40, 0) the same values as for xmax = ymax = 144. 22 On the lower (y = 0) and left (x = 0) boundary, the equation (71) degenerates to a HJB equation for the one-dimensional uncertain volatility model from [2]. We solve it with the one-dimensional Tree-Grid method, as in paper [4]. For the case that the values from outside of the computational domain [0, xmax ] × [0, ymax ] are needed, we artificially define the solution for x < 0, y < 0 to have the same values as the solution on the boundary: BC(x < 0, y) = BC(0, y), BC(x, y < 0) = BC(x, 0), BC(x < 0, y < 0) = 0. 5.2. Numerical results Here we present the experimental convergence results of the 2D Tree-Grid 275 method applied to the Black-Scholes model and the uncertain volatility model. We implemented our method in Python 1 and tested on an Intel Core i7-4770 CPU 3.40GHz computer with 16 GB RAM. We performed the simulations on four different sets of parameters: • Black-Scholes model with moderate volatility and correlation: σx = 0.3, 280 σy = 0.5, ρ = 0.4. Results of the simulation are in Table 1. • Black-Scholes model with extreme parameters: σx = 0.05, σy = 0.05, ρ = −0.95. Results of the simulation are in Table 2. • Uncertain volatility model, maximal option price (worst case scenario) with parameters σx,min = σy,min = ρmin = 0.3, σx,max = σy,max = 285 ρmax = 0.5. Results of the simulation are in Table 3. • Uncertain volatility model, minimal option price (best case scenario). The max operator is replaced by min in equation (71), all other parameters are the same as in the worst case scenario. Results of the simulation are in Table 4. 290 In all models we used the parameters T = 0.25, r = 0.05. For each model we computed the approximations of the solution on different refinement levels. Let 1 The code can be downloaded from https://github.com/igor-vyr/Tree-Grid-method 23 us denote the final time-layer (t = 0) of the approximation of the solution on the k-th refinement level as Ak , and final time-layer of the reference solution as Aref . We measured the approximation error on each refinement level in two 295 different ways: 1. the L1 error - the error was computed using the formula: Error Ak = kAk − Aref k1 , (76) 2. the error in (x = 40, y = 40), was computed using the formula Error Ak = |Ak (40, 40) − Aref (40, 40)|. (77) We use this error and the value in (x = 40, y = 40) to compare our results with results from the paper [7], where the same parameters are used for the uncertain volatility model. To compute the experimental order of convergence (EOC) we used the following formula: EOC Ak = log(ErrAk−1 ) − log(ErrAk ) . log(hk−1 ) − log(hk ) (78) In all refinement levels we used a (rectangular) grid with equidistant time300 stepping and non-equidistant space-stepping with more nodes near to the nonsmooth regions of the terminal conditions. The refinements were done uniformly. From Table 2 we can deduce, that the numerical solution converges also in the case of very small volatility and large correlation, although the convergence is not as smooth as in the case of moderate parameters (Table 1). As we can 305 see in Tables 1, 3, the point-wise convergence may be quite non-smooth, even if the approximation is converging relatively smoothly in L1 . In the uncertain volatility model, the values in (x = 40, y = 40) are similar to the values of Ma and Forsyth [7], what verifies our method. The experimental order of convergence in the Tables 1, 3, 4 seems to be slightly better than the theoretical order 310 1, the experimental order of convergence in Table 2 is quite non-smooth. The increase of the experimental order of convergence in the last refinement level in the uncertain volatility model results from using a solution computed on 24 Table 1: Black-Scholes model, σx = 0.3, σy = 0.5, ρ = 0.4. M -number of time steps, N -number of space nodes. Value, error and experimental order of convergence (EOC) in the final time layer in the point (x, y) = (40, 40) and error, EOC in L1 norm in the final time layer. As reference solution the exact solution (computed with the R-library fExoticOptions) was used. M 25 50 100 200 400 N 35 2 69 2 Convergence in (x, y) = (40, 40) Convergence in L1 Value Error EOC Error EOC 1.9910 1.77E-01 - 1.21E-02 - 1.8211 7.02E-03 4.66 1.84E-03 2.72 2 1.8229 8.88E-03 -0.34 7.16E-04 1.36 2 1.8177 3.67E-03 1.27 3.02E-04 1.25 137 273 2 545 1.8141 5.31E-05 6.11 1.24E-04 1.28 800 1089 2 1.8138 1.96E-04 -1.89 5.11E-05 1.28 1600 21772 1.8139 1.38E-04 0.51 2.54E-05 1.01 Table 2: Black-Scholes model, σx = 0.05, σy = 0.05, ρ = −0.95. M -number of time steps, N -number of space nodes. Value, error and experimental order of convergence (EOC) in the final time layer in the point (x, y) = (40, 40) and error, EOC in L1 norm in the final time layer. As reference solution the exact solution (computed with the R-library fExoticOptions) was used. M 25 50 N 35 2 69 2 Convergence in (x, y) = (40, 40) Convergence in L1 Value Error EOC Error EOC 3.3619 1.28E+00 - 3.60E-02 - 3.8702 7.71E-01 0.73 1.96E-02 0.88 2 4.3095 3.31E-01 1.22 7.62E-03 1.37 200 2 273 4.6339 6.91E-03 5.58 5.53E-04 3.78 400 5452 4.6465 5.73E-03 0.27 3.86E-05 3.84 800 10892 4.6468 5.97E-03 -0.06 4.61E-05 -0.26 1600 21772 4.6416 8.35E-04 2.84 9.85E-06 2.22 100 137 a fine grid as reference solution (that is disproportionally closer to the solution on the last refinement level in contrast to the solutions on previous refinement 25 Table 3: Uncertain volatility model, worst case scenario (maximization). M -number of time steps, N -number of space nodes, Q -number of controls. Value, error and experimental order of convergence (EOC) in the final time layer in the point (x, y) = (40, 40) and error, EOC in L1 norm in the final time layer. As reference solution, a solution computed on a fine grid with 400 time steps, 5452 space nodes and 256 controls was used. M 25 50 100 200 N 2 35 2 69 Convergence in (x, y) = (40, 40) Convergence in L1 Q Value Error EOC Error EOC 16 2.8364 1.59E-01 - 1.17E-02 - 32 2.6619 1.53E-02 3.38 3.64E-03 1.68 137 2 64 2.6784 1.19E-03 3.68 1.17E-03 1.64 273 2 128 2.6784 1.20E-03 -0.01 3.07E-04 1.93 Table 4: Uncertain volatility model, best case scenario (minimization). M -number of time steps, N -number of space nodes, Q -number of controls. Value, error and experimental order of convergence (EOC) in the final time layer in the point (x, y) = (40, 40) and error, EOC in L1 norm in the final time layer. As reference solution, a solution computed on a fine grid with 400 time steps, 5452 space nodes and 256 controls was used. M 25 50 100 200 315 N 2 35 2 69 Convergence in (x, y) = (40, 40) Convergence in L1 Q Value Error EOC Error EOC 16 0.9847 6.95E-02 - 4.29E-03 - 32 0.9475 3.23E-02 1.11 2.01E-03 1.09 137 2 64 0.9270 1.18E-02 1.46 7.59E-04 1.41 273 2 128 0.9173 2.07E-03 2.50 1.31E-04 2.53 levels). In Figure 2, we see the final time layer (t = 0) of the approximation of option prices under the uncertain volatility model for both best and worst case scenario. 26 Worst case scenario Best case scenario 0.8 0.6 0.6 0.4 0.4 0.5 0.2 0.2 0.0 0.0 2.0 1.5 1.5 1.0 1.0 0.5 0.0 60 y 40 20 0 0 20 60 40 x V(x,y) V(x,y) 0.8 2.0 100 80 1.0 1.0 2.5 2.5 100 80 100 80 0.0 60 y 40 20 0 0 20 60 40 x 100 80 Figure 2: Final time layers (t=0) of the approximations of the worst case option price (maximization) and of the best case option price (minimization) from the uncertain volatility model computed with the 2D Tree-Grid method with 50 time steps, 692 space nodes and 32 controls. 6. Conclusion 320 In this paper we generalized the one-dimensional Tree-Grid method to 2 space dimensions. Although the ideas behind this two-dimensional method are similar to the one-dimensional case, the algorithm construction is more challenging, mainly because of the correlation between the two spatial variables. Our 2D Tree-Grid method is explicit and hence suitable for parallelization, but 325 still unconditionally stable and monotone in contrast to other explicit methods [1, 9, 10]. Moreover, unlike other wide-stencil methods [7, 9, 10], it doesn’t use any interpolation. These properties make the method very flexible and, by simply following of the Algorithm 1, also easy to implement. We proved the convergence of the 2D Tree-Grid method using the theory of 330 Barles and Souganidis [15]. We tested the method on the two-factor uncertain volatility model and on the Black-Scholes model for option pricing and verified the convergence also experimentally. Acknowledgements This research was supported within the bilateral German-Slovakian Project 335 ENANEFA - Efficient Numerical Approximation of Nonlinear Equations in Financial Applications financed by DAAD and the Slovakian Ministry of Education (01/2018–12/2019). 27 References [1] H. Kushner, P. G. Dupuis, Numerical methods for stochastic control prob340 lems in continuous time, Vol. 24 of Applications of Mathematics, Springer Science & Business Media, 2013. [2] P. A. Forsyth, K. R. Vetzal, Numerical methods for nonlinear PDEs in finance, in: Handbook of Computational Finance, Springer, 2012, pp. 503– 528. 345 [3] S. Kilianová, D. Ševčovič, A transformation method for solving the Hamilton-Jacobi-Bellman equation for a constrained dynamic stochastic optimal allocation problem, The ANZIAM Journal 55 (01) (2013) 14–38. [4] I. Kossaczký, M. Ehrhardt, M. Günther, A new convergent explicit TreeGrid method for HJB equations in one space dimension, Applied Mathe- 350 matics Letters 11 (2018) 1–29. [5] I. Kossaczký, M. Ehrhardt, M. Günther, The Tree-Grid Method with control-independent stencil, in: Proceedings of Equadiff 2017 Conference, Bratislava, 2017, pp. 79–88. [6] J. Wang, P. A. Forsyth, Maximal use of central differencing for Hamilton- 355 Jacobi-Bellman PDEs in finance, SIAM Journal on Numerical Analysis 46 (3) (2008) 1580–1601. [7] K. Ma, P. Forsyth, An unconditionally monotone numerical scheme for the two-factor uncertain volatility model, IMA Journal of Numerical Analysis 37 (2) (2016) 905–944. 360 [8] Y. Chen, J. W. Wan, Monotone mixed narrow/wide stencil finite difference scheme for Monge-Ampère equation, arXiv preprint arXiv:1608.00644. [9] J. F. Bonnans, H. Zidani, Consistency of generalized finite difference schemes for the stochastic HJB equation, SIAM Journal on Numerical Analysis 41 (3) (2003) 1008–1021. 28 365 [10] K. Debrabant, E. Jakobsen, Semi-Lagrangian schemes for linear and fully non-linear diffusion equations, Mathematics of Computation 82 (283) (2013) 1433–1462. [11] K. Debrabant, E. R. Jakobsen, Semi-Lagrangian schemes for linear and fully non-linear Hamilton-Jacobi-Bellman equations, arXiv preprint 370 arXiv:1403.1217. [12] J. Yong, X. Y. Zhou, Stochastic controls: Hamiltonian systems and HJB equations, Vol. 43 of Stochastic Modelling and Applied Probability, Springer Science & Business Media, 1999. [13] S. Mataramvura, B. Øksendal, Risk minimizing portfolios and HJBI equa375 tions for stochastic differential games, Stochastics An International Journal of Probability and Stochastic Processes 80 (4) (2008) 317–337. [14] Q. Song, Convergence of Markov chain approximation on generalized HJB equation and its applications, Automatica 44 (3) (2008) 761–766. [15] G. Barles, P. E. Souganidis, Convergence of approximation schemes for fully 380 nonlinear second order equations, in: Asymptotic Analysis, no. 4, 1991, pp. 2347–2349. [16] M. G. Crandall, H. Ishii, P.-L. Lions, Users guide to viscosity solutions of second order Partial Differential Equations, Bulletin of the American Mathematical Society 27 (1) (1992) 1–67. 385 [17] P.-L. Lions, Optimal control of diffusion processes and Hamilton–Jacobi– Bellman equations part 2: viscosity solutions and uniqueness, Communications in Partial Differential Equations 8 (11) (1983) 1229–1276. [18] D. M. Pooley, P. A. Forsyth, K. R. Vetzal, Numerical convergence properties of option pricing PDEs with uncertain volatility, IMA Journal of 390 Numerical Analysis 23 (2) (2003) 241–267. [19] F. Black, M. Scholes, The pricing of options and corporate liabilities, Journal of Political Economy 81 (3) (1973) 637–654. 29