A fitted finite volume method for stochastic optimal control Problems

In this article, we provide a numerical method based on fitted finite volume method to approximate the Hamilton-Jacobi-Bellman (HJB) equation coming from stochastic optimal control problems. The computational challenge is due to the nature of the HJB equation, which may be a second-order degenerated partial differential equation coupled with optimization. In the work, we discretize the HJB equation using the fitted finite volume method and show that matrix resulting from spatial discretization is an M-matrix. The optimization problem is solved at every time step using iterative method. Numerical results are presented to show the robustness of the fitted finite volume numerical method comparing to the standard finite difference method.


Introduction
We consider the numerical approximation of the following controlled Stochastic Differential Equation (SDE) defined in R n (n ≥ 1) by is the drift term and the d-dimensional diffusion coefficients. Note that ω t are d-dimensional independent Brownian motion on (Ω, F, (F t ) t≥0 , P), α = (α t ) t≥0 is an F-adapted process, valued in A closed convex subset of R m (m ≥ 1) and satisfying some integrability conditions and/or state constraints. Precise assumptions on b and σ to ensure the existence of the unique solution x t of (1) can be found in [8].
The existence and uniqueness of the viscosity solution of the HJB equation (5) is well known and can be found in [8]. Equation (5) is a initial value problem. There are two unknown functions in this equation, the value function v and the optimal control α. In most practical situations, (5) is not analytically solvable therefore numerical approximations are the only tools appropriate to provide reasonable approximations. Numerical approximation of HJB-equation of type (5) is therefore an active research area and has attracted a lot of attentions [20,18,11,10,15,14,16,13,9].
While 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 degenerated PDE, we need to couple with an optimization problem at every grid point and every time step. In terms of existing numerical methods, there are two basic threads of literature concerning controlled 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 [6]. A more recent approach is based on numerical methods such as finite difference method which ensure convergence to the viscosity solution of the HJB equation couple with an optimization problem at each time [12] .
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 [7], chapter 26). This degeneracy also has an adverse impact on the accuracy of our stochastic optimal control problems since its numerical resolution implies the resolution of PDE, coupled with optimization problem.
In this article, we propose a numerical scheme based on a finite volume method suitable to handle the degeneracy of the linear operator while solving numerically the HJB equation in dimension 1 and 2. The method is coupled with implicit time-stepping method for temporal discretization method and the iterative method presented in [9] for optimization problem at every time step.
More precisely, this method is based on fitted finite volume technique proposed in [1] to solve the degenerated Black Sholes equations. Note that to the best of our knowledge, such method has not been used to solve the stochastic optimal control problem (5).
The merit of the method is that it is absolutely stable in time because of the implicit nature of the time discretisation and the corresponding matrix after spatial discretization is a positive-definite Numerical simulations prove that our proposed method is more accurate that the standard method based on finite difference spatial discretization.
The rest of this article is organized as follows. In section 2, we present the finite volume method with the fitting technique for dimension 1 and 2. We will also show that the system matrix of the resulting discrete equations is an M -matrix. In section 3, we will present the temporal discretization and optimization problem in dimension 1 and 2. Numerical experiments using Matlab software will be performed in section 4 to demonstrate the accuracy of the proposed numerical method. We conclude the work at section 5 by summarizing our finding.

Spatial discretization
As we already know, the resolution of the HJB equation (5) involves a spatial discretisation, a temporal discretisation 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 equation (5) solving our stochastic optimal control problem (4). Details in this section can be found in [4], where such methods have been used to solve the degenerated Black Sholes equation for option pricing with constant coefficients.

Spatial discretization based on fitted finite volume method in dimension 1
Consider the more generalized HJB equation (5) in dimension 1 (n = 1) which can be written in the form.
∂v(x, t) ∂t where a(t, x, α) > 0, α = α(x, t) and bounded. As usual, we truncate the problem in the We also set Applying the mid-points quadrature rule to the first and the last point terms, we obtain the above dv dt l i + sup 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 = 2, · · · N 1 −1. This discussion is divided into two cases for i ≥ 1 and I 0 = (0, x 1 ).
The term a(x, t, α i ) x ∂v ∂x + b(x, t, α i ) v is approximated by solving the boundary value problem a(x, t, α i ) x ∂v ∂x Integrating (11) yields the first-order linear equations where C 1 denotes an additive constant. As in [4], the solution is given by Note that in this deduction we have assumed that b(x i+1/2 , t, α i ) = 0. By setting using the boundary conditions in (11) yields Solving the following linear system with respect to C 1 and C 2 yields yields Case II: This is the degenerated zone. The aims here is to approximate ρ at x 1/2 in the sub-interval I 0 . In this case, the following problem is considered where C 2 is an unknown constant to be determined. Following [4], integrating (18) yields Since Therefore we have By replacing ρ by its approximated value, (9) becomes for i = 0, 1, · · · , N 1 − 1 By setting τ = 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 Dirichlet boundary and final conditions, (61) is an M-matrix for any α ∈ A.

Spatial discretization based on fitted finite volume method in dimension 2
Here we consider the following two dimensional problem where We will assume that a 21 = a 12 . We also also define the following coefficients, which will help us to build our scheme smoothly a 11 (x, y, t, α) = a(x, y, t, α) x 2 , a 22 (x, y, t, α) = a(x, y, t, α)y 2 and a 12 = a 21 = d 1 (x, y, t, α)xy.
As usual the two dimensional domain is truncated to I x × I y , where I x = [0, x max ] and I y = [0, y max ] be divided into N 1 and N 2 sub-intervals: with 0 = x 0 < x 1 < · · · · · · < x N 1 = x max and 0 = y 0 < y 1 < · · · · · · < y N 2 = y max . This defines a mesh on I x × I y with all the mesh lines perpendicular to one of the axes.
We also set and y N 2 +1/2 = y max . For each i = 0, 1, · · · , N 1 and j = 0, 1, · · · , N 2 , we set h We now discuss the finite volume method for (31). Integrating both size of (31) over R i,j = for i = 1, 2, · · · , N 1 − 1, j = 1, 2, · · · , N 2 − 1. Applying the mid-points quadrature rule to the first and the last point terms, we obtain the above and v i,j (t) denotes the nodal approximation to v(x i , y j , t). We now consider the approximation of the middle term in (33). Let n denote the unit vector outward-normal to ∂R i,j .
By Ostrogradski Theorem, integrating by parts and using the definition of flux k, we have To achieve this, it is clear that we now need to derive approximations of the k(v(x, y, t, α i,j )) · n defined above at the mid-point x i+1/2 , y j , of the interval I x i for i = 0, 1, · · · N 1 − 1. This discussion is divided into two cases for i ≥ 1 and i ∈ I 0 = (0, x 1 ). This is really an extension of the one dimensional fitted finite volume presented in the previous section.
Case I: For i ≥ 2.
Remember that a 11 (x, y, t, α) = a(x, y, t, α) x 2 , we approximate the term a 11 ∂v ∂x + x b 1 v by solv-ing the following two points boundary value problem Integrating (36) 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 β i,j (t) = b 1i+1/2,j (t, α i,j ) a i+1/2,j (t, α i,j ) and a 12 = a 21 = d 1 (x, y, t, α) x y. Finally, we use the forward difference, Finally, Simillary, the second term in (34) can be approximated by Case II: For j ≥ 2.
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 (36) is degenerated. Therefore, we reconsider the following form where C 2 is an unknown constant to be determined. Integrating (45), we can deduce that Case IV: Approximation of the flux at J 0 .

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. We will present it in one and two dimensional. Let us re-consider the differential equation coupled with optimization problem given in (22) or (62) by v(0) given, For temporal discretization, we use a constant time step ∆t > 0, of course variable time steps can be used. The temporal grid points given by ∆t = τ n+1 − τ n for n = 1, 2, . . . m − 1. We denote v(τ n ) ≈ v n , A n (α) = A(τ n , α) and G n (α) = G(τ n , α).
For θ ∈ 1 2 , 1 , following [9], the θ-Method approximation of (63) in time is given by As we can notice, to find the unknown v n+1 we need also to solve an optimization. Let Then, the unknown v n+1 is solution of the following equation Note that, for θ = 1 2 , we have the Crank Nickolson scheme and for θ = 1 we have the Implicit scheme. Unfortunately (64)-(65) are nonlinear and coupled and we need to iterate at every time step. The following iterative scheme close to the one in [9] is used.

Numerical experiments
The goal of this section is carried out on test problems in both 1 and 2 space dimensions to validate the numerical scheme presented in the previous section. All computations were performed in Matlab 2013 using the estimate parameters coming from [17] and [9]. We will present two problems with exact solution and one problem without exact solution modelling cash management in finance.
By identification, The two dimensional Antsaz exact solution [8] at (τ n , x i , y j ) is given by where v m is the numerical approximation of v computed from our numerical scheme. The 3 D graphs of the Implicit Fitted Finite Volume ( θ = 1 at the final time T = 1) with its corresponding exact solution is given at Figure 3 and Figure 2, with N 1 = 50, N 2 = 45, r 1 = 0.0449/2, µ 1 = 0.0657/2,    r 2 = 0.044/2, µ 2 = 0.065/2, σ = 0.2537/2 and p = 0.5255/2. We compare the fitted finite volume method and the finite difference method in Table 2. Again, we can observe the accuracy of the fitted scheme comparing to the finite difference scheme.
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 σ t with modeled by dσ t = α σ t dt + β σ t dW 2t .
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 ∈ [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 subject to We assume that the two Brownian motions are correlated, that is dW 1t dW 2t = ρ dt. For this problem of optimal Cash Management the analytical solution is not available, so our numerical scheme will to provide approximated solution. The corresponding HJB equation for this optimal cash management problem (85) is given by 1] {(f + β 1 R) J R + w (u r 2 + (1 − u) r 1 + u R − d(t)) J w + (86) 1/2 σ 2 J RR + β 2 σ 2 J σσ + 2 ρ β σ 2 J σ R + α σ J σ = 0, with terminal condition J(·, T ) = w T . To simplify the problem, we assume that J(w, R, σ, t) = wH(R, σ, t).

Conclusion
We presented a fitted finite volume method to solve the HJB equation from stochastic optimal control problems coupled with implicit temporal discretization. The optimization problem is solved at every time step using iterative method. It was shown that the corresponding system matrix is an M-matrix, so the maximum principle is preserved for the discrete system. Numerical experiments in 1 and 2 dimensions are performed to prove the accuracy of the fitted finite volume method comparing to the standard finite difference methods.