Either way – solving the SCPs directly or solving the HJB equation – numerical methods are often needed, as closed-form solutions are rarely known. For the one-dimensional (1D) 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 Kushner and Dupuis (2013). Implicit finite-difference methods (FDMs) for solving HJB equations were presented, for example, in Forsyth and Vetzal (2012). The advantage of these methods in contrast with the explicit FDMs is their unconditional stability and monotonicity. In Kilianová and Ševčovič (2013), a method based on a transformation of the HJB equation was developed. The advantage of using this approach is that we do not need to solve the optimization problem in each time layer. Recently, a new explicit but still unconditionally stable and monotone tree–grid method was developed in Kossaczký et al (2018). This method combines the tree structure from the trinomial tree methods for option pricing with the grid typically used in FDMs, and it can be seen as some special explicit FDM or a Markov chain approximation. A useful modification of the method was presented in Kossaczký et al (2017). As the explicitness and great flexibility of this method (unconditional stability with any grid spacing) are clear advantages, we are interested in generalizing it to two state-space dimensions (two controlled stochastic processes).

A generalization of the implicit unconditionally stable method from Wang and Forsyth (2008) in two space dimensions was presented in Ma and Forsyth (2016) and later used in Chen and Wan (2016). The main idea behind this method is to combine the wide and the fixed stencil, depending on the correlation in a particular time–space node. Alternatively, explicit methods based on the ideas of Kushner and Dupuis (2013) presented in Bonnans and Zidani (2003) and Debrabant and Jakobsen (2013, 2014) can be used. These are wide-stencil schemes that are stable under some Courant–Friedrichs–Lewy (CFL) condition. Moreover, bilinear interpolation of the grid values is needed in these methods. Our generalization of the tree–grid method presented in this paper also falls into the class of wide-stencil schemes; however, it is unconditionally stable and no two-dimensional (2D) interpolation of the grid values is needed. The unconditional stability is achieved by making the stencil size dependent not only on the space steps but also on the time step. In order to get a monotone seven-point-wide stencil, we possibly need to artificially increase the variance and decrease the covariance of the system. However, the approximation is still consistent, as this artificial increase/decrease vanishes with grid refinement and, in standard cases, even reaches zero at some point. Note that even in standard wide-stencil schemes (Debrabant and Jakobsen 2013) the variance is increased and the covariance is altered due to the use of interpolation. Moreover, in that case, this behavior is present for any grid refinement, which contrasts with the case of the tree–grid method presented in this work. As this 2D tree–grid method is explicit, it is also suitable for parallelization.

We are concerned with searching for the value function $V(x,y,t)$ of the following 2D general SCP:

where $x$ and $y$ are state variables and $t$ is time. Here, $\overline{\mathrm{\Theta }}$ is the space of all suitable control functions (see, for example, Kushner and Dupuis 2013; Yong and Zhou 1999) from ${ℝ}^2\times [0,T]$ to a set $\mathrm{\Theta }$. For our purpose, we will assume $\mathrm{\Theta }$ to be discrete. If this is not the case, we can easily achieve this property by its discretization. Now, following Bellman’s principle, the dynamic programming equation holds:

where are some time points and ${\overline{\mathrm{\Theta }}}_{t_j}$ is a set of control functions from $\overline{\mathrm{\Theta }}$ restricted to the ${ℝ}^2\times [t_j,t_{j+1})$ domain. Using (1.5), it can be shown (Yong and Zhou 1999) that solving the SCP (1.1)–(1.4) is equivalent to solving the 2D HJB equation:

where ${\sigma }_x(\cdot )$, ${\sigma }_y(\cdot )$, ${\sigma }_{xy}(\cdot )$, ${\mu }_x(\cdot )$, ${\mu }_y(\cdot )$, $r(\cdot )$ and $f(\cdot )$ are functions of $x$, $y$, $t$ and $\theta$. Note that the maximum operator in (1.1) and (1.6) can be replaced by a minimum operator and the entire following analysis will hold analogously. A possible generalization of the SCP with a corresponding HJB equation is the stochastic differential game with a corresponding Hamilton–Jacobi–Bellman–Isaac equation (Mataramvura and Øksendal 2008). Further generalizations can be found in Song (2008).

First, before discussing how to choose the stencil nodes around node $(x_i,y_j,t_k)$ and the proper weights, we present the notation used in this paper.

• •

  $x_1,x_2,\mathrm{\dots },x_{N_x}$: space discretization in one space dimension. $y_1,y_2,\mathrm{\dots },y_{N_y}$: space discretization in two space dimensions. $t_1,t_2,\mathrm{\dots },t_M$: time discretization.

• •

  ${\mathrm{\Delta }}_{k}t=t_{k+1}-t_k$: (current) time step. We will use equidistant time stepping ${\mathrm{\Delta }}_{k}t=\mathrm{\Delta }t$ for all $k=1,2,\mathrm{\dots },M-1$. Generalization to nonequidistant time stepping is straightforward.

• •

  ${\mathrm{\Delta }}_{i}x=x_{i+1}-x_i$, ${\mathrm{\Delta }}_{j}y=y_{j+1}-y_j$: space steps in one dimension and two dimensions (possibly nonequidistant).

• •

  $\mathrm{\Delta }x={\max }_i{\mathrm{\Delta }}_{i}x$, $\mathrm{\Delta }y={\max }_i{\mathrm{\Delta }}_{i}y$.

• •

  $h=\max (K\max (\mathrm{\Delta }x,\mathrm{\Delta }y),\mathrm{\Delta }t)$, where $K>0$ is a parameter used for regulating the stencil size. In the nonequidistant case, $\mathrm{\Delta }t$ should be replaced by ${\mathrm{\Delta }}_{k}t$.

• •

  $b=h/\mathrm{\Delta }t$: in the nonequidistant case, this should be replaced by $b=h/{\mathrm{\Delta }}_{k}t$ in the following algorithm.

• •

  $(x_i,y_j,t_k)$: the node for which the stencil is constructed.

• •

  $E_{*}={\mu }_{*}(x_i,y_j,t_k,\theta )\mathrm{\Delta }t$ for $*=x,y$.

• •

  ${\mathrm{Var}}_{*}=({\sigma }_{*}^2(x_i,y_j,t_k,\theta )+?(h))\mathrm{\Delta }t$ for $*=x,y$ (to be determined later).

• •

  $\rho (x_i,y_j,t_k,\theta )={\sigma }_{xy}(x_i,y_j,t_k,\theta )/({\sigma }_x(x_i,y_j,t_k,\theta ){\sigma }_y(x_i,y_j,t_k,\theta ))$.

• •

  ${\tilde{\sigma }}_{xy}(x_i,y_j,t_k,\theta )={\sigma }_{xy}(x_i,y_j,t_k,\theta )+?(\sqrt{h})$ (to be determined later).

• •

  $\mathrm{Cov}=({\sigma }_{xy}(x_i,y_j,t_k,\theta )+?(\sqrt{h}))\mathrm{\Delta }t$ (to be determined later).

• •

  $W_{*}={\mathrm{Var}}_{*}+E_{*}^2$ for $*=x,y$.

• •

  $W_{xy}=\mathrm{Cov}+E_x{E}_y$.

• •

  $v_{i,j}^k$: numerical approximation of $V(x_i,y_j,t_k)$.

Next, we describe how to choose the stencil nodes around an arbitrary node $(x_i,y_j,t_k)$ for a fixed control $\theta$. The values in these nodes will affect the value in the node $(x_i,y_j,t_k)$.

If $E_x>0$,

else if ,

else ($E_x=0$)

where ${\lceil \cdot \rceil }_x$ (respectively, ${\lfloor \cdot \rfloor }_x$) denotes rounding to the nearest greater (respectively, smaller) element from $x_1,x_2,\mathrm{\dots },x_{N_x}$. If such an element does not exist, ${\lceil z\rceil }_x$ (respectively, ${\lfloor z\rfloor }_x$) will return just $z$.

If $E_y>0$,

else if ,

else ($E_y=0$)

where ${\lceil \cdot \rceil }_y$ and ${\lfloor \cdot \rfloor }_y$ are defined analogously. Note that the operators ${\lfloor \cdot \rfloor }_x$, ${\lceil \cdot \rceil }_x$, ${\lfloor \cdot \rfloor }_y$ and ${\lceil \cdot \rceil }_y$ formally correspond to some special (left or right neighbor) constant interpolation in one dimension. However, in the standard wide-stencil schemes for 2D problems, constant, (bi)linear or higher-order interpolations in two dimensions are used. In such cases, the use of higher-order interpolations is motivated by a need to reduce the interpolation error. In our scheme, however, this 1D constant interpolation does not introduce any additional interpolation error, as it is only needed to identify the stencil nodes. As we will see later, the moments of our approximation match exactly with the moments of the original stochastic process for a nondegenerate problem and small enough $h$.

The following nodes with their respective weights (probabilities) will be used in the stencil located at $(x_i,y_j,t_k)$:

• •

  node $(x_i,y_j,t_{k+1})$ with probability $p_o$;

• •

  nodes $(x_{+},y_j,t_{k+1})$ and $(x_{-},y_j,t_{k+1})$ with probabilities $p_{x+}$ and $p_{x-}$;

• •

  nodes $(x_i,y_{+},t_{k+1})$ and $(x_i,y_{-},t_{k+1})$ with probabilities $p_{y+}$ and $p_{y-}$; and

• •

  nodes $(x_{+},y_{+},t_{k+1})$ and $(x_{-},y_{-},t_{k+1})$, both with probability $p_{xy}$ if $W_{xy}$ is nonnegative, and nodes $(x_{+},y_{-},t_{k+1})$ and $(x_{-},y_{+},t_{k+1})$, both with probability $p_{xy}$ if $W_{xy}$ is negative (in the following algorithm, we will focus on the case of nonnegative $W_{xy}$; the case of negative $W_{xy}$ is treated analogously).

In the next section, we will impose six conditions on these probabilities in equation form in order to match with the moments of the discrete and the original process as well as ensure that the probabilities sum up to $1$. Therefore, to get a unique solution to these equations, we need, in the seven nodes, only six different probabilities, which explains why we chose the same probability for nodes $(x_{+},y_{+},t_{k+1})$ and $(x_{-},y_{-},t_{k+1})$ (respectively, $(x_{+},y_{-},t_{k+1})$ and $(x_{-},y_{+},t_{k+1})$). Our motivation for choosing these particular nodes to have the same probability is that we try to handle the drift and diffusion in the $x$ (respectively, $y$) direction in a similar manner to the 1D tree–grid method (Kossaczký et al 2018, 2017), using the nodes on the $x$ (respectively, $y$) axis with origin at $(x_i,y_j)$. Then, only one probability remains for the two remaining nodes.

First, let us define the following difference operators:

To match with 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 the following:

For negative $W_{xy}$, only the condition (2.14) changes to

The solutions of (2.9)–(2.14) are

where ${\mathrm{\Delta }}_c={\mathrm{\Delta }}_{+}x{\mathrm{\Delta }}_{+}y+{\mathrm{\Delta }}_{-}x{\mathrm{\Delta }}_{-}y$ for nonnegative $({\tilde{\sigma }}_{xy}(x_i,y_j,t_k,\theta )\mathrm{\Delta }t+E_x{E}_y)$ and ${\mathrm{\Delta }}_c={\mathrm{\Delta }}_{+}x{\mathrm{\Delta }}_{-}y+{\mathrm{\Delta }}_{-}x{\mathrm{\Delta }}_{+}y$ for negative $({\tilde{\sigma }}_{xy}(x_i,y_j,t_k,\theta )\mathrm{\Delta }t+E_x{E}_y)$.

Following the construction of $x_{+}$, $x_{-}$, $y_{+}$ and $y_{-}$, it is easy to check that $p_o$ is always nonnegative. The same holds for $p_{xy}$. To ensure the nonnegativity of $p_{x+}$, $p_{x-}$, $p_{y+}$ and $p_{y-}$, we have to properly choose the variables ${\mathrm{Var}}_x$, ${\mathrm{Var}}_y$ and $\mathrm{Cov}$ to get nonnegative weights while remaining consistent with the original problem as $h\to 0$. This is done in the next section.

Let us assume $E_x\ge 0$. Now, the first fraction of $p_{x+}$ is positive, but the first fraction of $p_{x-}$ may be negative. We will set ${\mathrm{Var}}_x$ (and, hence, $W_x$) in such way that the latter will be also positive. It holds that

The right-hand side of the inequality (2.23) is $0$ for ${\mathrm{Var}}_x=A_x$ and greater than $0$ for ${\mathrm{Var}}_x>A_x$ with

Here, we replaced $E_x$ with $|E_x|$ to cover the analogous case and possibly negative $p_{x+}$. Now, if we set

where $A_y$ is defined analogously, the first fraction in $p_{x+}$, $p_{x-}$, $p_{y+}$ and $p_{y-}$ will be nonnegative. The possible difference between ${\mathrm{Var}}_x$ (respectively, ${\mathrm{Var}}_y$) and the variances from the original problem ${\sigma }_x^2\mathrm{\Delta }t$ (respectively, ${\sigma }_y^2\mathrm{\Delta }t$) is called artificial diffusion. Now, taking into account the correlation coefficient in the current node and following these definitions of variances ${\mathrm{Var}}_x$ and ${\mathrm{Var}}_y$, we define a new covariance:

and

Now, we define $\mathrm{Cov}$ as

This covariance $\mathrm{Cov}$ is consistent with variances ${\mathrm{Var}}_x$ and ${\mathrm{Var}}_y$ defined by (2.25) and (2.26). As for ${\tilde{\sigma }}_{xy}(x_i,y_j,t_k,\theta )\mathrm{\Delta }t+E_x{E}_y\ge 0$, it holds that

and for it holds that

Now, it also holds that $|W_{xy}|=C_{xy}$, which implies that $p_{x+}$, $p_{x-}$, $p_{y+}$ and $p_{y-}$ are all positive. It is easy to check that the following holds:

In addition, for ${\sigma }_x^2(x_i,y_j,t_k,\theta )\ne 0$ (respectively, ${\sigma }_y^2(x_i,y_j,t_k,\theta )\ne 0$) and small enough $h$,

holds. It also holds that

It follows that

and the same holds for the second, third and fourth maximum candidates in (2.28). Moreover,

Therefore,

and it is easy to check for ${\sigma }_x^2\ne 0$, ${\sigma }_y^2\ne 0$, $|{\sigma }_{xy}|\ne {\sigma }_x{\sigma }_y$ and small enough $h$ that

holds. Following (2.32) and (2.39), the modified variances and modified covariance are consistent with the variances and covariance from the original problem. Moreover, for ${\sigma }_x,{\sigma }_y\ne 0$ and , the modified variances and covariance will be equal to their originals for small enough $h$.

In the formula for $h=\max (K\max (\mathrm{\Delta }x,\mathrm{\Delta }y),\mathrm{\Delta }t)$, the $K\max (\mathrm{\Delta }x,\mathrm{\Delta }y)$ part is responsible for keeping the correlation in the numerical model consistent with the correlation of the original problem. Here, the parameter $K>0$ can be chosen arbitrarily; however, for a large (relative to problem parameters) $K$, the stencil is also large, which typically increases the error. On the other hand, for too small $K$, the correlation may start being exact (or exact enough) for very fine grids only and not sufficiently exact for coarse grids, resulting in larger errors on these coarse grids. For $K=0$, we cannot guarantee the convergence of the correlation. However, the correlation in the numerical model 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. (1)

   set $l=1$ and use $K=K_0$, $K_0\ge 0$ to compute $x_{*}$, $y_{*}$, $p_{*}$

2. (2)

   compute the correlation of the numerical model, $\tilde{\rho }=\mathrm{Cov}/\sqrt{{\mathrm{Var}}_x{\mathrm{Var}}_y}$

3. (3)

   if $\tilde{\rho }\ne \rho$, recompute $x_{*}$, $y_{*}$, $p_{*}$ with $K=K_l$, $K_l>K_{l-1}$ and set $l:=l+1$

   else break.

4. (4)

   if , go to step (2)

   else break.

With this modification, we will use a smaller stencil size as much as possible. This approach can be seen as somewhat analogous to the approach in Ma and Forsyth (2016), where a fixed (and thus small) stencil is used as much as possible. However, we will not increase the convergence rate, but we might reduce the error.

Another approach, the nonconstant $K$ modification, involves using nonconstant $K=K(x,y,\theta )$. This can be smaller in nodes with large volatilities ${\sigma }_x$ and ${\sigma }_y$, and larger in nodes with smaller ${\sigma }_x$ and ${\sigma }_y$, thus regulating the stencil size to not explode in case of large volatilities.

Both modifications can also be combined, and it is easy to check that they do not harm the consistency. In our numerical simulation, these modifications did not lead to significant improvements; therefore, we decided to just use a constant $K=1/400$. However, they may be useful for other SCPs.