A fitted finite volume method for stochastic optimal control problems in finance

1.   Introduction

We consider the numerical approximation of the following controlled Stochastic Differential Equation (SDE) defined in $ \mathbb{R}^{ n} $ ($ n\geq 1 $) by

where

is the drift term and

the $ d $-dimensional diffusion coefficients. Note that $ W_t $ are $ d $-dimensional independent Brownian motion on $ \left(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \geq 0}, \mathbb{P} \right) $, $ \alpha = (\alpha_t)_{t \geq 0} $ is an $ \mathbb{F} $-adapted process, valued in $ \mathcal{A} $ compact convex subset of $ \mathbb{R}^m \, (m \geq 1) $ is the set of admissible controls satisfying some integrability conditions and/or state constraints. Precise assumptions on $ b $ and $ \sigma $ to ensure the existence of the uniqueness solution $ x_t $ of (1.1) can be found in ^[9].

Given a function $ g :\mathbb{R}^n \rightarrow \mathbb{R} $ and $ f :[0, T]\times \mathbb{R}^n \times \mathcal{A} \rightarrow \mathbb{R} $, the value function is defined by

Remember that $ x $ is the initial condition in (1.1). Using the dynamic programming approach, the problem (1.2) can be transformed into the following Hamilton-Jacobi-Bellman (HJB) equation (see ^[5,12] for more explanations)

where

where $ a(t, x, \alpha) = (\dfrac{1}{2}(\sigma (t, x, \alpha))(\sigma(t, x, \alpha))^T)_{i, j} $ the symmetric positive semi-definite diffusion coefficient matrix. Eq (1.3) is a initial value problem. Although this initial value problem (1.3) is defined on the unbounded region $ \mathbb{R}^n $, often it is restricted to a bounded region

with $ p_k $ and $ q_k $ are constants for $ k = 1, 2, \cdots, n $ for computational reasons. This initial value problem with $ \mathbb{R}^n $ in (1.3) being replaced by $ \mathcal{U} $ has been discussed in the literature (see for example ^[11,21]). There are two unknown functions in this equation, the value function $ v $ and the optimal control $ \alpha $. However, in most practical situations, it is not analytically solvable therefore numerical approximations are the only tools appropriate to provide reasonable approximations. Numerical approximation of HJB-equation of type (1.3) is therefore an active research area and has attracted a lot of attentions ^{[11,13,14,18,20]}. When solving numerically HJB equation, the keys challenge are the low regularity of the solution of HJB equation and the lack of appropriate numerical methods to tackle the degeneracy of the differential operator in HJB equation. Indeed adding to the standard issue that we usually have when solving degenerate PDE, we need to couple with an optimization problem at every grid point and every time step. In practice, a feedback optimal control solution is defined in a region containing the optimal trajectory is much preferred. This feedback solution gives a global optimal control defined over a time-space region.

In terms of existing numerical methods, there are three basic threads of literature concerning HJB equations. A standard approach is based on Markov chain approximation. In financial terms, this approach is equivalent to an explicit finite difference method. However, these methods are well-known to suffer from time step limitations due to stability issues ^[7]. A more recent approach is called semi-Lagrangian schemes, these schemes work for general diffusions with coefficient matrices that may be nondiagonal dominant and arbitrarily degenerate but also use finite difference approximations ^[22].

For many stochastic optimal control problems such as Merton's control problem, the linear operator is degenerated when the spatial variables approach the region near to zero. This degeneracy has an adverse impact on the accuracy when the finite difference method is used to solve the PDE (see ^[8], chapter 26) as the monotonicity of the scheme is usually lost. However, when solving HJB equation, the monotonicity also plays a key role to ensure the convergence of the numerical scheme toward the viscosity solution. Indeed in the two dimensional Merton's control problem, the matrix in the diffusion part is of rank 1 near the origin and it has been found in ^[2,3] that the standard finite difference schemes become non monotone and may not converge to the viscosity solution of the HJB.

The current work aims to propose an alternative novel monotone scheme based on fitted technique in dimensions $ 1 $ and $ 2 $. This fitting technique is based on the idea proposed by ^[1] for convection-diffusion equations and was upgraded in ^[4] to solve simple degenerated Black Scholes equations. The fitted technique has for feature to tackle the degeneracy in the HJB equation. Our method is coupled with implicit time-stepping method for temporal discretization method and the iterative method presented in ^[11] for optimization problem at every time step. Note that to the best of our knowledge, such method has not been used to solve the stochastic optimal control problem (1.3). The merit of the method is that it is absolutely stable in time because of the implicit nature of the time discretization, easy to understand and implement, and the corresponding matrices after spatial and temporal discretization are positive-definite $ M $-matrices. Therefore our alternative scheme is monotone. Numerical simulations prove that our proposed method can be more accurate than the standard method based on finite difference spatial discretization, thanks to the monotonicity properties of our scheme. The novel contribution of our paper over the existing literature can be summarized as

● We have upgraded the fitted finite volume technique to spatial discretization of a more generalized HJB equation coupled with the implicit time-stepping method for temporal discretization method and the iterative method for the optimization problem at every time step. To best of our knowledge such combination has not yet proposed so far to solve stochastic optimal control problems in the literature.

● We have proved that the corresponding matrices after spatial and temporal discretization are positive-definite $ M $-matrices, then our scheme is monotone. Indeed this is the key feature of our novel scheme comparing to standard finite difference scheme where the monotonicity is lost because of degeneracy and the convergence toward the viscosity solution is not longer guarantee ^[2,3].

● We have demonstrated by numerical experiments that the proposed scheme can be more accurate than the standard finite difference scheme, and that accuracy increase in two dimensional domain where the matrix in the diffusion part is of rank $ 1 $ near the origin. We have also applied our method to approximate realistic optimal problem in finance, more precisely the optimal cash management problem.

The rest of this article is organized as follows. In section 2, we will recall some results about the well posedness problem. In section 3, we present the finite volume method with the fitting technique for dimensions $ 1 $ and $ 2 $. We will also show that the system matrix of the resulting discrete equations is an $ M $-matrix. In section 4, we will present the temporal discretization and optimization problem in dimensions $ 1 $ and $ 2 $. Numerical experiments using Matlab software will be performed in section 5 to demonstrate the accuracy of the proposed numerical method. We conclude the work at section 6 by summarizing our finding.

2.   Well posedness problem

In this part, we will present results about existence and uniqueness of the value function. Firstly, assume that $ b $ and $ \sigma $ are continuous and for every $ \alpha \, \in\, \mathcal{A} $, $ b(\cdot, \cdot, \alpha) $ and $ \sigma(\cdot, \cdot, \alpha) $ are in $ C^1([0, T] \times \mathbb{R}^n) $. To ensure the existence and uniqueness of the value function solution, we make the following assumptions.

Assumption 1. We assume that there exists $ C \geq 0 $ such that for all $ \alpha \, \in\, \mathcal{A} $, $ x, y \, \in\, \mathbb{R}^n $ and $ t, s \, \in\, [0, T], $

Assumption 2. We assume that there exists $ C \geq 0 $ such that for every $ \alpha \, \in\, \mathcal{A} $, $ x \, \in\, \mathbb{R}^n $ and $ t, s \, \in\, [0, T], $

Let $ T > 0 $, $ \mathcal{U}\subset \mathbb{R}^n $ be an open bounded set and $ S(n, \mathbb{R}) $ the set of symmetric matrix $ n\times n $. We set $ \mathcal{O} = (0, T) \times\mathcal{U} $ and $ \partial\mathcal{O} $ the parabolic boundary of $ \mathcal{O} $ defined by : $ \partial\mathcal{O} = \partial \mathcal{U}\times (0, T)\cup (\overline{\mathcal{U}})\times \left\lbrace T \right\rbrace $.

Assumption 3. Let $ f : \overline{\mathcal{O}} \times \mathcal{A} \longrightarrow \mathbb{R} $ and $ g :\overline{\mathcal{O}} \longrightarrow \mathbb{R} $ be two continuous functions such that there exists $ C \geq 0 $ such that, for all $ (t, x), \, (s, y)\, \in\, \overline{\mathcal{O}} $, $ \alpha\, \in\, \mathcal{A} $

The value function $ v $ from $ \overline{\mathcal{O}} \rightarrow \mathbb{R} $ is now given by

where the diffusion process $ x_s $ is on the form (1.1). Using the dynamic programming approach (see ^[12]), the value function $ v $ leads to the following HJB equation

where

$ a(t, x, \alpha) = \bigg(\dfrac{1}{2}(\sigma (t, x, \alpha))(\sigma(t, x, \alpha))^T \bigg)_{i, j} $ the symmetric positive semi definite diffusion coefficient matrix.

Let us consider the following equation without final condition

We denote

Definition 1. ^[17] A function $ u\, \in\, USC \left(\overline{\mathcal{O}}\right) $ is a viscosity subsolution of Eq (2.7) if and only if $ u $ is such that for every test function $ \rho \in C^{\infty} \left(\overline{\mathcal{O}}\right) $, $ u-\rho $ has a strict local maximum at $ (t, x) \in \mathcal{O} $ with $ u(t, x) = \rho(t, x) $, implies

A function $ u\, \in\, LSC \left(\overline{\mathcal{O}}\right) $ is a viscosity supersolution of Eq (2.7) if and only if $ u $ is such that for every test function $ \rho \in C^{\infty} \left(\overline{\mathcal{O}}\right) $, $ u-\rho $ has a strict local minimum at $ (t, x) \in \mathcal{O} $ with $ u(t, x) = \rho(t, x) $, implies

A function $ u\, \in\, C\left(\overline{\mathcal{O}}\right) $ is a viscosity solution of Eq (2.7) if it is both, a viscosity subsolution and a viscosity supersolution.

Having stated the notion of a viscosity solution to a parabolic PDE, we now turn towards the notion of a viscosity solution to a parabolic final value problem with Dirichlet boundary data on the parabolic boundary.

Definition 2. A function $ v\, \in\, USC \left(\overline{\mathcal{O}}\right) $ is a viscosity subsolution of (2.11) if $ v $ is a viscosity subsolution in the sense of the definition (1) and $ v \leq g $.

Likewise, a function $ v\, \in\, LSC \left(\overline{\mathcal{O}}\right) $ is a viscosity supersolution of (2.11) if $ v $ is a viscosity subsolution in the sense of the definition (1) and $ v \geq g $ on $ \partial \mathcal{O} $.

A function $ v $ on $ C\left(\overline{\mathcal{O}}\right) $ is a viscosity solution if $ v $ is a viscosity supersolution and a viscosity subsolution.

Theorem 1. Let assumptions (1), (2) and (3) hold, there is at most one viscosity solution in the sense of definition (2) to the HJB final value problem.

Proof. The proof of existence and uniqueness can be found in ^[15,17,19].

3.   Spatial discretization

As we already know, the resolution of the HJB Eq (2.13) involves a spatial discretization, a temporal discretization and an optimisation problem at every grid point and each time step. The goal of this section is to provide the spatial discretization of the HJB Eq (2.13) solving our stochastic optimal control problem. Details in this section can be found in ^[6], where such methods have been used to solve the degenerated Black Scholes equation for option pricing with constant coefficients.

3.1.   Spatial discretization based on fitted finite volume method in dimension 1

Consider the more generalized HJB Eq (2.13) in dimensions 1 ($ n = 1 $) which can be written in the form

where $ a(t, x, \alpha) > 0 $ is bounded. Indeed this divergence form is not a restriction as the differentiation is respect to $ x $ and not respect to the control $ \alpha $, which may be discontinuous in some applications.

As usual, we truncate the problem in the finite interval $ I = [0, x_{\text{max}}] $. Let the interval $ I = [0, x_{\text{max}}] $ be divided into $ N_1 $ sub-intervals $ I_i : = (x_i, x_{i+1}), \, \, \, \, \, i = 0\cdots N_1-1 $ with $ 0 = x_0 < x_1 < \cdots \cdots < x_{N_1} = x_\text{max} $. We also set $ x_{i+1/2} = \dfrac{x_{i} +x_{i+1} }{2} $ and $ x_{i-1/2} = \dfrac{x_{i-1} +x_{i} }{2} $ for each $ i = 1\cdots N_1-1 $. These mid-points form a second partition of $ [0, x_{\text{max}}] $ if we define $ x_{-1/2} = x_0 $ and $ x_{N_1+1/2} = x_\text{max} $. Integrating both size of (3.1) over $ J_i = \left(x_{i-1/2}, \, x_{i+1/2}\;\right) $ and taking $ \alpha_i = \alpha(x_i, t) $, we have

Applying the mid-points quadrature rule to the first and the last point terms, we obtain the above

for $ i = 1, 2, \cdots N_1-1 $, where $ l_i = x_{i+1/2} - x_{i-1/2} $ is the length of $ J_i $. Note that $ v_i $ denotes the nodal approximation to $ v(t, x_i) $ and $ \rho(v) $ is the flux associated with $ v $ defined by

Clearly, we now need to derive approximation of the flux defined above at the mid-point $ x_{i+1/2} $, of the interval $ I_i $ for $ i = 0, 1, \cdots N_1-1. $ This discussion is divided into two cases for $ i \geq 1 $ and $ i = 0 $ on $ I_0 = [0, x_1] $.

$\underline{\bf Case\; I}:$ Approximation of $ \rho $ at $ x_{i+1/2} $ for $ i\geq 1 $.

The term $ \left(a (t, x, \alpha_i)\, x \dfrac{\partial v}{\partial x}+ b(t, x, \alpha_i)\, v \right) $ is approximated by solving the two point boundary value problem

Integrating (3.5) yields the first-order linear equations

where $ C_1 $ denotes an additive constant. As in ^[6], the solution is given by

Note that in this deduction we have assumed that $ b(t, x_{i+1/2}\;, \alpha_i) \neq 0 $. By setting $ \beta_i(t) = \dfrac{b(t, x_{i+1/2}\;, \alpha_i)}{a(t, x_{i+1/2}\;, \alpha_i)} $, using the boundary conditions in (3.5) yields

Solving the following linear system with respect to $ C_1 $ and $ C_2 $ yields

yields

$ \rho_i (v)(t) $ provides an approximation to the $ \rho(v)(t) $ at $ x_{i+1/2} $. Similarly the approximation of $ \rho(v)(t) $ at $ x_{i-1/2} $ is given by

$\underline{\bf Case\; II}:$ This is the degenerated zone. The aims here is to approximate $ \rho $ at $ x_{1/2} $ in the sub-interval $ I_0 = [0, x_1] $ for $ i = 0 $. In this case, the following problem is considered

where $ C_2 $ is an unknown constant to be determined. Following ^[6], integrating (3.12) yields

Since $ x_{1/2} = \dfrac{x_1+ x_0}{2} $ with $ x_0 = 0 $, we have $ C_2\, x_1 = (a(x_{1/2}, t, \alpha_1)+b(t, x_{1/2}, \alpha_1))(v_1(t)-v_0(t)) $. Therefore we have

By replacing $ \rho $ by its approximated value 3.10, 3.11 and 3.14, (3.3) becomes for $ i = 0, 1, \cdots, N_1-1 $

By setting $ \tau = T-t $ and including the boundary conditions, we have the following system of Ordinary Differential Equation (ODE) coupled with optimisation problem.

which can be rewritten as

where $ \mathcal{A}^{N_{1}-1} = \underset{N_{1}-1}{\underbrace{\mathcal{A} \times \mathcal{A}\times \cdots \times \mathcal{A}}} $, $ \textbf{v}(\tau) = (v_1(\tau), \cdots, v_{N_1-1}(\tau)) $ and $ F(\alpha, \tau) = (F_1(\alpha_1, \tau), \cdots, F_{N_1-1}(\alpha_{N_1-1}, \tau)) $ includes all Dirichlet boundary and final conditions, $ A (\alpha, \tau) = -E(\alpha, \tau) $ and $ G(\alpha, \tau) = - F(\alpha, \tau) $ are defined as for $ i = 1, \cdots, N_1-1 $

$ G(\alpha, \tau) = \left[\begin{array}{c} -\dfrac{1}{4\, l_1}\, x_{1}(a_{1/2}(\tau, \alpha_1)-b_{1/2}(\tau, \alpha_1))\, v_0 \\ 0\\ \vdots\\ 0\\ -x_{N_1-1/2}\dfrac{b_{N_1-1/2}(\tau, \alpha_{N_1-1})\, x_{N_1}^{\beta_{N_1-1}\; (\tau)} }{l_{N_1-1}\, (x_{N_1}^{\beta_{N_1-1}\; (\tau)}-x_{N_1-1}^{\beta_{N_1-1}\; (\tau)})}v_{N_1} \end{array} \right]. $

Theorem 2. Assume that $ c $ given in (3.1) is negative and let $ h = \underset{1\leq i\leq N_1} {\max} l_i $. If $ h $ is relatively small then the matrix $ E (\alpha, \tau) $ in the system (3.17) is an M-matrix for any $ \alpha \, \in\, \mathcal{A}^{N_{1}-1} $.

Proof. Let us show that $ E(\alpha, \tau) $ has positive diagonal, non-positive off diagonal, and is diagonally dominant. We first note that

for $ i = 1, \cdots, N_1-1 $, and all $ b_{i+1/2}(\tau, \alpha) \neq 0, \; b_{i-1/2}(\tau, \alpha) \neq 0, $ with $ a_{i+1/2}(\tau, \alpha) > 0 $ and $ a_{i-1/2}(\tau, \alpha) > 0 $.

This also holds when $ b_{i+1/2}(\tau, \alpha) \rightarrow 0 $ and $ b_{i-1/2}(\tau, \alpha) \rightarrow 0 $, that is

Using the definition of $ E (\alpha, \tau) $ given above, we see that

because $ x_{i+1}^{\beta_i(\tau)} \approx x_{i}^{\beta_i(\tau)} + x_{i}^{\beta_i(\tau)-1}\, {\beta_i(\tau)}\, h $, $ x_{i-1}^{\beta_{i-1}\;(\tau)} \approx x_{i}^{\beta_{i-1}\;(\tau)} - x_{i}^{\beta_{i-1}\;(\tau)-1}\, {\beta_{i-1}(\tau)}\, h $ and

Note that for $ i = 1 $, we have $ E_{1, 1}\geq 0 $ if $ {a_{1/2}(\tau, \alpha_1)} + b_{1/2}(\tau, \alpha_1) $, are nonnegative and $ c_1(\tau, \alpha_1) < 0 $. So $ E(\alpha, \tau) $ is diagonally dominant and is therefore an $ \textbf{M} $-matrix.

Indeed the assumption $ c < 0 $ can be relaxed, but if positive it should not be more that a certain threshold $ c_{0} > 0 $ depending of our numerical scheme. Indeed this is not restrictive as from numerical experiments the threshold $ c_{0} > 0 $ is always the best comparing to finite difference scheme where its not monotonicity does not depend on coefficients.

3.2.   Spatial discretization based on fitted finite volume method in dimension $ 2 $

Here we consider the generalized HJB Eq (2.13) in dimension 2 which can be written in the form

where $ k (t, x, y, \alpha) = A(t, x, y, \alpha)\cdot\nabla v(t, x, y)+ b\, v(t, x, y) $ is the flux,

We will assume that $ a_{21} = a_{12} $. We also define the following coefficients, which will help us to build our scheme smoothly $ a_{11}(t, x, y, \alpha) = a(t, x, y, \alpha)\, x^2, a_{22}(t, x, y, \alpha) = \overline{a}(t, x, y, \alpha) y^2 $ and $ a_{1 2} = a_{2 1} = d_1(t, x, y, \alpha)\, x\, y $. Although this initial value problem (3.25) is defined on the unbounded region $ \mathbb{R}^2 $, we often restrict our consideration to a bounded region. As usual the two dimensional domain is truncated to $ I_{x} \times I_{y} $, where $ I_{x} = [0, x_{\text{max}}] $ and $ I_{y} = [0, y_{\text{max}}] $ be divided into $ N_1 $ and $ N_2 $ sub-intervals:

with $ 0 = x_{0} < x_{1} < \cdots \cdots < x_{N_1} = x_{\text{max}} $ and $ 0 = y_{0} < y_1 < \cdots \cdots < y_{N_2} = y_{\text{max}} $. This defines a mesh on $ I_{x} \times I_{y} $ with all the mesh lines perpendicular to one of the axes.

We also set

for each $ i = 1, \cdots, N_1-1 $ and each $ j = 1, \cdots, N_2-1 $. We denote $ N = (N_1-1)\times (N_2-1) $. These mid-points form a second partition of $ I_{x} \times I_{y} $ if we define $ x_{-1/2} = x_{0} $, $ x_{N_1+1/2} = x_{\text{max}} $, $ y_{-1/2} = y_{0} $ and $ y_{N_2+1/2} = y_{\text{max}} $. For each $ i = 0, 1, \cdots, N_1 $ and $ j = 0, 1, \cdots, N_2 $, we set $ h_{x_i} = x_{i+1/2} - x_{i-1/2} $ and $ h_{y_j} = y_{j+1/2} - y_{j-1/2} $.

We now discuss the finite volume method for (3.25). Integrating both size of (3.25) over $ \mathcal{R}_{i, j} = \left[ x_{i-1/2}, x_{i+1/2} \right] \times \left[ y_{j-1/2}, y_{j+1/2}\right] $, we have

for $ i = 1, 2, \cdots, N_1-1 $, $ j = 1, 2, \cdots, N_2-1 $. Applying the mid-points quadrature rule to the first and the last point terms, we obtain the above

for $ i = 1, 2, \cdots N_1-1 $, $ j = 1, 2, \cdots N_2-1 $ where $ l_{i, j} = \left(x_{i+1/2}-x_{i-1/2}\; \right) \times \left(y_{j+1/2}-y_{j-1/2}\right) $ is the length of $ \mathcal{R}_{i, j} $, and $ v_{i, j}(t) $ denotes the nodal approximation to $ v(t, x_{i}, y_{j}) $. We now consider the approximation of the middle term in (3.28). Let $ \bf n $ denote the unit vector outward-normal to $ \partial \mathcal{R}_{i, j} $. By Ostrogradski Theorem, integrating by parts and using the definition of flux $ k(v) $, we have

We shall look at (3.29) term by term. For the first term we want to approximate the integral by a constant as

To achieve this, it is clear that we now need to derive approximations of the $ k(v(x, y, t, \alpha_{i, j})) \cdot \bf n $ defined above at the mid-point $ \left(x_{i+1/2}\;, y_{j} \right) $, of the interval $ I_{x_{i}} $ for $ i = 0, 1, \cdots N_1-1 $. This discussion is divided into two cases for $ i \geq 1 $ and $ i = 0 $ on $ I_0 = [0, x_{1}] $. This is really an extension of the one dimensional fitted finite volume presented in the previous section.

$\underline{\bf Case\; I}:$ For $ i\geq 1 $.

Remember that $ a_{11}(t, x, y, \alpha) = a(t, x, y, \alpha)\, x^2 $, we approximate the term $ \left(a_{11} \dfrac{\partial v}{\partial x}+ x\, b_1\, v \right) $ by solving the following two points boundary value problem

Integrating (3.31) yields the first-order linear equations

where $ C_1 $ denotes an additive constant. Following the one dimensional fitted finite volume presented in the previous section, we have

Therefore,

where $ {b_1}(t, x_{i+1/2}\;, y_j, \alpha_{i, j}) = {b_1}_{i+1/2, j}(t, \alpha_{i, j}) $, $ {a}(t, x_{i+1/2}\;, y_j, \alpha_{i, j}) = {a}_{i+1/2, j}(t, \alpha_{i, j}) $, $ \beta_{i, j}(t) = \dfrac{{b_1}_{i+1/2, j}(t, \alpha_{i, j})}{{a}_{i+1/2, j}(t, \alpha_{i, j})} $ and $ a_{12} = a_{21} = d_{1}(t, x, y, \alpha)\, x\, y $. Finally, we use the forward difference,

Finally,

Similarly, the second term in (3.29) can be approximated by

$\underline{\bf Case\; II}:$ For $ j\geq 1 $.

For the third term we want to approximate the integral by a constant, that is

Following the first case of this section, we have

Similarly, the fourth term in (3.29) can be approximated by

for $ j = 1, \cdots, N_2-1 $, where $ \bar{\beta}_{i, j}(t) = \dfrac{{b_2}_{i, j+1/2}(t, \alpha_{i, j})}{\bar{a}_{i, j+1/2}(t, \alpha_{i, j})} $ with $ a_{22}(t, x, y, \alpha) = \bar{a}(t, x, y, \alpha)\, y^2 $.

$\underline{\bf Case\; III}:$ Approximation of the flux at $ I_0 $. Note that the analysis in case I does not apply to the approximation of the flux on $ [0, x_1] $ because (3.31) is degenerated. Therefore, we reconsider the following form

where $ C_2 $ is an unknown constant to be determined. Integrating (3.39), we can deduce that

$\underline{\bf Case\; IV}:$ Approximation of the flux at $ J_0 $.

Using the same procedure for the approximation of the flux at $ I_0 $, we deduce that

By replacing the flux by his value for $ i = 1, \cdots, N_1-1 $ and $ j = 1, \cdots, N_2-1 $, Eq (3.28) becomes

By setting $ \tau = T-t $ and including the boundary conditions, we have the following system

where for $ i = 1, \cdots, N_1-1 $, $ j = 1, \cdots, N_2-1 $ and $ N = (N_1-1)\times(N_2-1) $, $ \mathcal{A}^{N} = \underset{(N_{1}-1)\times(N_2-1)}{\underbrace{\mathcal{A} \times \mathcal{A}\times \cdots \times \mathcal{A}}} $. We have

As for one dimension case, (3.43) can be rewritten as the ODE coupled with optimization problem

or

where $ A(\tau, \alpha) = - E (\tau, \alpha) $, $ E (\tau, \alpha) (I, J) = \left(e_{i', j'}^{i, j}\right) $, $ i', i = 1, \cdots, N_1-1 $, $ j', j = 1, \cdots, N_2-1 $, $ n_1 = N_1-1, \; n_2 = N_2-1; \; \; I: = I (i, j) = i + (j-1)n_1 $ and $ J: = J (i', j') = i' + (j'-1)n_1 $, $ \textbf{v} = \left({v}_{1, 1}, \cdots, {v}_{1, N_2-1}, \cdots, {v}_{N_1-1, 1}, \cdots, {v}_{N_1-1, N_2-1}\right) $ and $ G (\tau, \alpha) = -F(\tau, \alpha) $ includes boundary condition.

Theorem 3. Assume that the coefficients of $ A $ given by (3.26) are positive, $ c < 0 $ and let $ h = \underset{\underset{1\leq j\leq N_2}{1\leq i\leq N_1}} {\max} \{ h_{x_i}, \, h_{y_j} \} $. If $ h $ is relatively small then the matrix $ E (\tau, \alpha) = \left(e_{i, j}^{i, j}\right)_{\underset{j = 1, \cdots, N_2-1, }{i = 1, \cdots, N_1-1}} $ in (3.44) is an M-matrix for any $ \alpha\, \in\, \mathcal{A}^N $.

Proof. The Proof follows the same lines of that of Theorem 2.

4.   Temporal discretization and optimization

This section is devoted to the numerical time discretization method for the spatially discretized optimization problem using the fitted finite volume method. Let us re-consider the differential equation coupled with optimization problem given in (3.16) or (3.45) by

For temporal discretization, we use a constant time step $ \Delta t > 0 $, of course variable time steps can be used. The temporal grid points given by $ \Delta t = \tau_{n+1}-\tau_n $ for $ n = 1, 2, \ldots m-1 $. We denote $ \textbf{v}(\tau_n) \approx \textbf{v}^n $, $ A^n (\alpha) = A(\tau_n, \alpha) $ and $ G^n (\alpha) = G(\tau_n, \alpha). $

For $ \theta \, \in \left[\frac{1}{2}, 1\right] $, following ^[11], the $ \theta $-Method approximation of (4.1) in time is given by

this also can be written as

As we can notice, to find the unknown $ \textbf{v}^{n+1} $ we need also to solve an optimization problem. Let

Then, the unknown $ \textbf{v}^{n+1} $ is solution of the following equation

Note that when $ \theta = \frac{1}{2} $ the time-stepping scheme becomes the Crank-Nicolson scheme and when $ \theta = 1 $ it is the Implicit scheme. Both of these schemes are unconditionally stable, and they are second and first-order accuracy respectively. Unfortunately (4.2) and (4.3) are nonlinear and coupled and we need to iterate at every time step. The following iterative scheme close to the one in ^[11] is used.

1. Let $ \left(\textbf{v}^{n+1}\right)^0 = \textbf{v}^{n} $,

2. Let $ \hat{\textbf{v}}^{k} = \left(\textbf{v}^{n+1}\right)^k $,

3. For $ k = 0, 1, 2 \cdots $ until convergence ($ \Vert \hat{\textbf{v}}^{k+1}-\hat{\textbf{v}}^{k}\Vert \leq \epsilon $, given tolerance) solve

4. Let $ k_l $ being the last iteration in step 3, set $ \textbf{v}^{n+1}: = \hat{\textbf{v}}^{k_l} $, $ \alpha^{n+1}: = \alpha^{k_l} $.

Indeed to find the solution at time $ \tau_{n+1} $, we use the solution at time $ \tau_n $ and the iterations in our algorithm. We have used the Matlab optimisation toolbox and more precisely the Matlab function " fminbnd" in (4.5). Note that the optimization problem (4.5) is solved in each grid point at every iteration.

The monotonicity of system matrix of (4.2), more precisely $ [I + \Delta t \, \theta E^{n+1}] $ is given in the following theorem.

Theorem 4. Under the same assumptions as in Theorem 2 and 3, for any given $ n = 1, 2, \cdots, m - 1 $, the system matrix $ [I + \Delta t \, \theta E^{n+1}] $ in (4.2) is an M–matrix for each $ \alpha \in \mathcal{A}^N. $

Proof. The proof is obvious. Indeed as in Theorem 2 and 3, $ [I + \Delta t \, \theta E^{n+1}] $ is (strictly) diagonally dominant since $ \Delta t > 0 $. Then, it is an $ M $–matrix.

Note that a simple finite difference scheme generally does not give an M-matrix. The merit of the proposed method is that it is unconditionally stable in time because of the implicit nature of the time discretization. More precisely, following [10,Theorem 6 and Lemma 3], we can easily prove that the scheme (4.4) is stable and consistent, so the convergence of the scheme is ensured (see ^[23]).

5.   Numerical experiments

The goal of this section is carried out test problems in both $ 1 $ and $ 2 $ spatial dimensions to validate the numerical scheme presented in the previous section. All computations were performed in Matlab 2013. In our numerical experiments, we found that the fitted finite volume method and the finite difference method converge after $ 2 $ or $ 3 $ iterations in our algorithm (4.5). We have used for tolerance $ \epsilon = 10^{-6} $ in our iterative algorithm.

We will present two problems with exact solution and one problem without exact solution modelling cash management in finance.

Problem 1. Consider the following Merton's stochastic control problem such that $ \alpha = \alpha(t) $ is a feedback control belongs in $ [0, 1] $

subject to

where $ b^{\alpha_t}(t, x_t) = x_t\, \left[ r + \alpha_t\, (\mu-r)\right] $, $ \sigma^{\alpha_t} (t, x_t) = x_t\, \alpha_t \, \sigma $, $ 0 < p < 1 $ is coefficient of risk aversion, $ r $, $ \mu $, $ \sigma $ are constants, $ x_t\, \in\, \mathbb{R} $ and $ W_t $ Brownian motion. We assume $ \mu > r $. For this problem, using dynamic programming the corresponding HJB equation is given by

where

and $ a^\alpha (t, x) = \dfrac{1}{2}\bigg(\sigma^{\alpha}(t, x)\bigg)^2 $.

The divergence form of the HJB (5.2) is given by

The domain where we compare the solution is $ \Omega = \left[0, x_{{max}}\right] $, where Dirichlet boundary conditions is used at the boundaries. Of course the value of the boundary conditions are taken to be equal to the exact solution. The exact solution given in ^[9] is given at every $ (\tau^n, x_i) $ by

We use the following $ L^2([0, T]\times \Omega) $ norm of the absolute error

where $ v^m $ is the numerical approximation of $ v $ computed from our numerical scheme. The 3 D graphs of the implicit fitted finite volume ($ \theta = 1 $) with its corresponding exact solution are given at Figures 1 and 2 with parameters coming from ^[16]. For our computation, we have $ [0, 10] $ for computational domain with $ N = 1500 $, $ r = 0.0449 $, $ \mu = 0.0657 $, $ \sigma = 0.2537 $, $ p = 0.5255 $ and $ T = 1 $.

We compare the fitted finite volume method and the finite difference method in Table 1.

From Tables 1–3, we can observe that the implicit fitted finite volume is slightly accurate comparing to the implicit finite difference method. Both fitted finite volume and finite difference method converge with order $ 0.9 $ in time.

Problem 2. Consider the following two dimensional Merton's stochastic control model such that $ \alpha_1 = \alpha_1(t) $ and $ \alpha_2 = \alpha_2(t) $ are a feedback control in $ [0, 1] $

where

$ 0 < p < 1 $ is coefficient of risk aversion, $ r_1 $, $ \mu_1 $, $ r_2 $, $ \mu_2 \; \sigma $ are constants, $ x_t, \, y_t\, \in\, \mathbb{R} $ and $ \rho_{12} \in [0, 1) $ the correlation of the two Brownian motion. We assume that $ \mu_1 > r_1 $ and $ \mu_2 > r_2 $. For this problem, using dynamic programming the corresponding HJB equation is given by

where

and the two dimensional divergence form is given by

where $ k(t, x, y, \alpha) = A(t, x, y, \alpha)\nabla v(t, x, y)+ b(t, x, y, \alpha) v(t, x, y) $

By identification,

The two dimensional Ansatz exact solution ^[9] at $ \left(\tau^n, x_i, y_j\right) $ is given by

We use the following $ L^2([0, T]\times \Omega), \, \, \Omega = [0, x_{{\text{max}}}] \times [0, y_{{\text{max}}}] $ norm of the absolute error

where $ v^m $ is the numerical approximation of $ v $ computed from our numerical scheme.

The 3 D graphs of the implicit fitted finite volume ($ \theta = 1 $ at the final time $ T = 1 $) with its corresponding exact solution are given at Figures 3 and 4, with $ N_1 = 25 $, $ N_2 = 20 $, $ r_1 = 0.0449/2 $, $ \mu_1 = 0.0657/2 $, $ r_2 = 0.0448/2 $, $ \mu_2 = 0.0656/2 $, $ \sigma = 0.2537/2 $, $ p = 0.5255/2 $ and $ \rho_{12} = 0.9 $.

We compare the fitted finite volume method and the finite difference method in Table 4. As in dimension 1, we can observe the accuracy of the fitted scheme comparing to the finite difference scheme, thanks to the fitted technique. We can also observe that the accuracy is high when the number of time subdivisions is small. For $ \rho_{12} = 0.9 $ the fitted finite volume and the finite difference converge with order $ 1.3 $ in time.

Problem 3. We consider a optimal cash management under a stochastic volatility Model problem coming from ^[21]. We assume that the firm invests its cash in a bank account and a stock in a portfolio of value $ w_t $ at time $ t $, and the proportion of wealth invested in the stock at time $ t $ is $ u_t $. The interest rate earned in the bank account is $ r_1 $, the return from the stock at time $ t $ has two components, the cash dividend rate $ r_2 $, the capital gain rate $ R_t $. The dynamic of the capital gain rate $ R_t $ is assumed to be governed by the stochastic process

and the volatility $ \sigma_t $ with modeled by

Suppose that the firm has a demand rate $ d(t) $ for cash at time $ t $, and that the demand rate $ d(t) $ is normally distributed with mean $ 0 $ and variance $ 0.2 $. We assume that $ u_t\, \in\, [0, 1] $ and the wealth dynamics for this cash management problem is given by

The objective of the firm is to maximize the expectation of the total holdings at the terminal time $ T $. The portfolio optimization problem is given by

We assume that the two Brownian motions are correlated, that is $ d{W_{1}}_t\, d{W_{2}}_t = \rho\, dt $. For this problem of optimal cash management the analytical solution is not available, so our numerical scheme will to provide approximated solution. Using dynamic programming, the corresponding HJB equation for this optimal cash management problem (5.16) is given by

with terminal condition $ J(T, \cdot) = w_T. $ To simplify the problem, we assume that

Therefore (5.17) is equivalent to

with terminal condition $ H(T, R, \sigma) = 1 $. The HJB Eq (5.18) is a problem with two state variables $ R $ and $ \sigma $. The divergence form of the problem (5.18) is then given by

where $ k(t, R, \sigma, u) = A(t, R, \sigma, u)\nabla H(t, R, \sigma)+ b(t, R, \sigma) H(t, R, \sigma) $

By identification,

Because we have a stochastic volatility model, to solve the PDE equation, we have considered the following boundary conditions of Heston model

Because the PDE has two second derivatives in the two spatial directions, four boundary conditions are needed. This comes from the fact that the two second order derivatives give rise to two unknown integration constants. To meet this requirement, at the boundary $ \sigma = 0 $ it is considered inserting $ \sigma = 0 $ into the PDE to complete the set of four boundary conditions:

The HJB equation is solved using some parameters values in ^[21] given in the following tabular

Figure 5 shows a sample of fitted finite volume solution of the wealth rate $ H $ at the point $ (1/2, 1/2) $ from $ t = 1 $ to $ t = 10 $ with $ N_1 = 10 $, $ N_2 = 10 $, $ R_{\text{max}} = 1/2 $, $ \sigma_{\text{max}} = 1/2 $. We can estimate the mean and moment of $ H $ using Monte Carlo Method by generating many samples of $ H $. Table 6 presents the optimal investment policy from $ t = 1 $ to $ t = 10 $.

6.   Conclusions

We presented a fitted finite volume method to solve the HJB equation from stochastic optimal control problems coupled with implicit temporal discretization. It was shown that the corresponding system matrices are M-matrices, so the maximum principle is preserved for the discrete system. Numerical experiments in $ 1 $ and $ 2 $ dimensions are performed to prove that the fitted scheme can be more accurate comparing to the standard finite difference methods. We intent to extend our method in high dimension.

Acknowledgments

The first author was supported by the project African Center of Excellence in Mathematical Sciences and Applications (ACE-MSA) at IMSP Benin and the European Mathematical Society (EMS).

Conflict of interest

The authors declare no potential conflict of interests.