Loading [MathJax]/jax/output/SVG/jax.js
Research article

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

  • Received: 13 October 2020 Accepted: 28 December 2020 Published: 12 January 2021
  • MSC : 65M75

  • 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 in one and two dimensional domain. The computational challenge is due to the nature of the HJB equation, which may be a second-order degenerate partial differential equation coupled with optimization. For such problems, standard scheme such as finite difference losses its monotonicity and therefore the convergence toward the viscosity solution may not be guarantee. In the work, we discretize the HJB equation using the fitted finite volume method, which has for main feature to tackle the degeneracy of the equation. The time discretisation is performed using the Implicit Euler method, which is unconditionally stable. We show that matrices resulting from spatial discretization and temporal discretization are M-matrices. 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.

    Citation: Christelle Dleuna Nyoumbi, Antoine Tambue. A fitted finite volume method for stochastic optimal control problems in finance[J]. AIMS Mathematics, 2021, 6(4): 3053-3079. doi: 10.3934/math.2021186

    Related Papers:

    [1] Zuliang Lu, Ruixiang Xu, Chunjuan Hou, Lu Xing . A priori error estimates of finite volume element method for bilinear parabolic optimal control problem. AIMS Mathematics, 2023, 8(8): 19374-19390. doi: 10.3934/math.2023988
    [2] Zui-Cha Deng, Fan-Li Liu, Liu Yang . Numerical simulations for initial value inversion problem in a two-dimensional degenerate parabolic equation. AIMS Mathematics, 2021, 6(4): 3080-3104. doi: 10.3934/math.2021187
    [3] Dennis Llemit, Jose Maria Escaner IV . Value functions in a regime switching jump diffusion with delay market model. AIMS Mathematics, 2021, 6(10): 11595-11609. doi: 10.3934/math.2021673
    [4] Tiantian Zhang, Wenwen Xu, Xindong Li, Yan Wang . Multipoint flux mixed finite element method for parabolic optimal control problems. AIMS Mathematics, 2022, 7(9): 17461-17474. doi: 10.3934/math.2022962
    [5] Yu Xu, Youjun Deng, Dong Wei . Numerical solution of forward and inverse problems of heat conduction in multi-layered media. AIMS Mathematics, 2025, 10(3): 6144-6167. doi: 10.3934/math.2025280
    [6] Shuangbing Guo, Xiliang Lu, Zhiyue Zhang . Finite element method for an eigenvalue optimization problem of the Schrödinger operator. AIMS Mathematics, 2022, 7(4): 5049-5071. doi: 10.3934/math.2022281
    [7] Lidiya Melnikova, Valeriy Rozenberg . One dynamical input reconstruction problem: tuning of solving algorithm via numerical experiments. AIMS Mathematics, 2019, 4(3): 699-713. doi: 10.3934/math.2019.3.699
    [8] Zuliang Lu, Fei Cai, Ruixiang Xu, Chunjuan Hou, Xiankui Wu, Yin Yang . A posteriori error estimates of hp spectral element method for parabolic optimal control problems. AIMS Mathematics, 2022, 7(4): 5220-5240. doi: 10.3934/math.2022291
    [9] Yuelong Tang . Error estimates of mixed finite elements combined with Crank-Nicolson scheme for parabolic control problems. AIMS Mathematics, 2023, 8(5): 12506-12519. doi: 10.3934/math.2023628
    [10] Jiali Wu, Maoning Tang, Qingxin Meng . A stochastic linear-quadratic optimal control problem with jumps in an infinite horizon. AIMS Mathematics, 2023, 8(2): 4042-4078. doi: 10.3934/math.2023202
  • 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 in one and two dimensional domain. The computational challenge is due to the nature of the HJB equation, which may be a second-order degenerate partial differential equation coupled with optimization. For such problems, standard scheme such as finite difference losses its monotonicity and therefore the convergence toward the viscosity solution may not be guarantee. In the work, we discretize the HJB equation using the fitted finite volume method, which has for main feature to tackle the degeneracy of the equation. The time discretisation is performed using the Implicit Euler method, which is unconditionally stable. We show that matrices resulting from spatial discretization and temporal discretization are M-matrices. 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.


    We consider the numerical approximation of the following controlled Stochastic Differential Equation (SDE) defined in Rn (n1) by

    dxs=b(s,xs,αs)ds+σ(s,xs,αs)dWs,s(t,T]xt=x (1.1)

    where

    b:[0,T]×Rn×ARn(t,xt,αt))b(t,xt,αt)

    is the drift term and

    σ:[0,T]×Rn×ARn×d(t,xt,αt))σ(t,xt,αt)

    the d-dimensional diffusion coefficients. Note that Wt are d-dimensional independent Brownian motion on (Ω,F,(Ft)t0,P), α=(αt)t0 is an F-adapted process, valued in A compact convex subset of Rm(m1) is the set of admissible controls satisfying some integrability conditions and/or state constraints. Precise assumptions on b and σ to ensure the existence of the uniqueness solution xt of (1.1) can be found in [9].

    Given a function g:RnR and f:[0,T]×Rn×AR, the value function is defined by

    v(t,x)=supαAE[Ttf(s,xs,α)ds+g(xT)],(t,x)[0,T]×Rn. (1.2)

    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)

    {vt(t,x)+supαA[Lαv(t,x)+f(t,x,α)]=0on [0,T)×Rnv(T,x)=g(x),xRn (1.3)

    where

    Lαv(t,x)=ni=1(b(t,x,α))iv(t,x)xi+ni,j=1(a(t,x,α))i,j2v(t,x)xixj, (1.4)

    where a(t,x,α)=(12(σ(t,x,α))(σ(t,x,α))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 Rn, often it is restricted to a bounded region

    U=(p1,q1)×(p1,q1)××(pn,qn), (1.5)

    with pk and qk are constants for k=1,2,,n for computational reasons. This initial value problem with Rn in (1.3) being replaced by 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 α. 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.

    In this part, we will present results about existence and uniqueness of the value function. Firstly, assume that b and σ are continuous and for every αA, b(,,α) and σ(,,α) are in C1([0,T]×Rn). To ensure the existence and uniqueness of the value function solution, we make the following assumptions.

    Assumption 1. We assume that there exists C0 such that for all αA, x,yRn and t,s[0,T],

    |b(t,x,α)b(s,y,α)|C(|xy|+|ts|)|σ(t,x,α)σ(s,y,α)|C(|xy|+|ts|) (2.1)

    Assumption 2. We assume that there exists C0 such that for every αA, xRn and t,s[0,T],

    |b(t,x,α)|C(1+|x|)|σ(t,x,α)|C(1+|x|). (2.2)

    Let T>0, URn be an open bounded set and S(n,R) the set of symmetric matrix n×n. We set O=(0,T)×U and O the parabolic boundary of O defined by : O=U×(0,T)(¯U)×{T}.

    Assumption 3. Let f:¯O×AR and g:¯OR be two continuous functions such that there exists C0 such that, for all (t,x),(s,y)¯O, αA

    |f(t,x,α)f(s,y,α)|C(|xy|+|ts|)|f(t,x,α)|C(1+|x|)|g(t,x)|C(1+|x|). (2.3)

    The value function v from ¯OR is now given by

    v(t,x,α)=supαAE[τtf(s,xs,α(s))ds+g(τ,xτ)], (2.4)

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

    {vt(t,x)+supαA[Lαv(t,x)+f(t,x,α)]=0on Ov(T,x)=g(x)on O, (2.5)

    where

    Lαv(t,x)=ni=1(b(t,x,α))iv(t,x)xi+ni,j=1(a(t,x,α))i,j2v(t,x)xixj, (2.6)

    a(t,x,α)=(12(σ(t,x,α))(σ(t,x,α))T)i,j the symmetric positive semi definite diffusion coefficient matrix.

    Let us consider the following equation without final condition

    vt(t,x)+supαA[Lαv(t,x)+f(t,x,α)]=0on O (2.7)

    We denote

    USC(ˉO)={v:¯OR|vupper semi-continuous on¯O},LSC(¯O)={v:¯OR|vlower semi-continuous on¯O}. (2.8)

    Definition 1. [17] A function uUSC(¯O) is a viscosity subsolution of Eq (2.7) if and only if u is such that for every test function ρC(¯O), uρ has a strict local maximum at (t,x)O with u(t,x)=ρ(t,x), implies

    tρ(t,x)+supαA[Lαρ(t,x)+f(t,x,α)]0, (2.9)

    A function uLSC(¯O) is a viscosity supersolution of Eq (2.7) if and only if u is such that for every test function ρC(¯O), uρ has a strict local minimum at (t,x)O with u(t,x)=ρ(t,x), implies

    tρ(t,x)+supαA[Lαρ(t,x)+f(t,x,α)]0, (2.10)

    A function uC(¯O) 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.

    vt(t,x)+supαA[Lαv(t,x)+f(t,x,α)]=0on O (2.11)
    v(T,x)=g(x)on O, (2.12)

    Definition 2. A function vUSC(¯O) is a viscosity subsolution of (2.11) if v is a viscosity subsolution in the sense of the definition (1) and vg.

    Likewise, a function vLSC(¯O) is a viscosity supersolution of (2.11) if v is a viscosity subsolution in the sense of the definition (1) and vg on O.

    A function v on C(¯O) 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.

    {vt(t,x)+supαA[Lαv(t,x)+f(t,x,α)]=0on  Ov(T,x)=g(x),onO. (2.13)

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

    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.

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

    v(t,x)t+supαA[x(a(t,x,α)x2v(t,x)x+b(t,x,α)xv(t,x))+c(t,x,α)v(t,x)]=0 (3.1)

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

    As usual, we truncate the problem in the finite interval I=[0,xmax]. Let the interval I=[0,xmax] be divided into N1 sub-intervals Ii:=(xi,xi+1),i=0N11 with 0=x0<x1<<xN1=xmax. We also set xi+1/2=xi+xi+12 and xi1/2=xi1+xi2 for each i=1N11. These mid-points form a second partition of [0,xmax] if we define x1/2=x0 and xN1+1/2=xmax. Integrating both size of (3.1) over Ji=(xi1/2,xi+1/2) and taking αi=α(xi,t), we have

    xi+1/2xi1/2vtdx+xi+1/2xi1/2supαiA[xx(a(t,x,αi)xvx+b(t,x,αi)v)+c(t,x,αi)v]dx=0 (3.2)

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

    vi(t)tli+supαiA[[xi+1/2ρ(v)|xi+1/2xi1/2ρ(v)|xi1/2]+c(t,xi,αi)vili]=0, (3.3)

    for i=1,2,N11, where li=xi+1/2xi1/2 is the length of Ji. Note that vi denotes the nodal approximation to v(t,xi) and ρ(v) is the flux associated with v defined by

    ρ(v):=a(t,x,αi)xvx+b(t,x,αi)v. (3.4)

    Clearly, we now need to derive approximation of the flux defined above at the mid-point xi+1/2, of the interval Ii for i=0,1,N11. This discussion is divided into two cases for i1 and i=0 on I0=[0,x1].

    CaseI_: Approximation of ρ at xi+1/2 for i1.

    The term (a(t,x,αi)xvx+b(t,x,αi)v) is approximated by solving the two point boundary value problem

    (a(t,xi+1/2,αi)xvx+b(t,xi+1/2,αi)v)=0,xIiv(xi)=vi(t),v(xi+1)=vi+1(t). (3.5)

    Integrating (3.5) yields the first-order linear equations

    ρi(v)(t)=a(t,xi+1/2,αi)xvx+b(t,xi+1/2,αi)v=C1 (3.6)

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

    v(t)=C1b(t,xi+1/2,αi)+C2xb(t,xi+1/2,αi)a(t,xi+1/2,αi). (3.7)

    Note that in this deduction we have assumed that b(t,xi+1/2,αi)0. By setting βi(t)=b(t,xi+1/2,αi)a(t,xi+1/2,αi), using the boundary conditions in (3.5) yields

    vi(t)=C1b(t,xi+1/2,αi)+C2xβi(t)iandvi+1(t)=C1b(t,xi+1/2,αi)+C2xβi(t)i+1. (3.8)

    Solving the following linear system with respect to C1 and C2 yields

    {vi(t)=C1b(t,xi+1/2,αi)+C2xβi(t)ivi+1(t)=C1b(t,xi+1/2,αi)+C2xβi(t)i+1 (3.9)

    yields

    ρi(v)(t)=C1=b(t,xi+1/2,αi)(xβi(t)i+1vi+1(t)xβi(t)ivi(t))xβi(t)i+1xβi(t)i (3.10)

    ρi(v)(t) provides an approximation to the ρ(v)(t) at xi+1/2. Similarly the approximation of ρ(v)(t) at xi1/2 is given by

    ρ(v)(t)|xi1/2=b(t,xi1/2,αi)(xβi1(t)ivi(t)xβi1(t)i1vi1(t))xβi1(t)ixβi1(t)i1 (3.11)

    CaseII_: This is the degenerated zone. The aims here is to approximate ρ at x1/2 in the sub-interval I0=[0,x1] for i=0. In this case, the following problem is considered

    (a(t,x1/2,α1)xvx+b(t,x1/2,α1)v)=C2in[0,x1]v(0)=v0(t),v(x1)=v1(t) (3.12)

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

    ρ0(v)|1/2(t)=a(t,x1/2,αi)x1/2vx+b(t,x1/2,α1)v=b(t,x1/2,α1)v0(t)+C2x1/2. (3.13)

    Since x1/2=x1+x02 with x0=0, we have C2x1=(a(x1/2,t,α1)+b(t,x1/2,α1))(v1(t)v0(t)). Therefore we have

    ρ0(v)|1/2(t)=12[(a(t,x1/2,α1)+b(t,x1/2,α1))v1(t)(a(t,x1/2,α1)b(x1/2,t,α1))v0(t)]. (3.14)

    By replacing ρ by its approximated value 3.10, 3.11 and 3.14, (3.3) becomes for i=0,1,,N11

    dvi(t)dt+supαiA1li[xi+1/2ρi(v)(t)|xi+1/2xi1/2ρi(v)(t)|xi1/2+ci(t,αi)vi(t)li]=0 (3.15)

    By setting τ=Tt and including the boundary conditions, we have the following system of Ordinary Differential Equation (ODE) coupled with optimisation problem.

    {dv(τ)dτ+supαAN11[A(α,τ)v(τ)+G(α,τ)]=0v(0) given , (3.16)

    which can be rewritten as

    {dv(τ)dτ+infαAN11[E(α,τ)v(τ)+F(α,τ)]=0v(0) given,  (3.17)

    where AN11=A×A××AN11, v(τ)=(v1(τ),,vN11(τ)) and F(α,τ)=(F1(α1,τ),,FN11(αN11,τ)) includes all Dirichlet boundary and final conditions, A(α,τ)=E(α,τ) and G(α,τ)=F(α,τ) are defined as for i=1,,N11

    Ei,i+1(αi,τ)=xi+1/2bi+1/2(τ,αi)xβi(τ)i+1li(xβi(τ)i+1xβi(τ)i), (3.18)
    Ei,i(αi,τ)=(xi+1/2bi+1/2(τ,αi)xβi(τ)ili(xβi(τ)i+1xβi(τ)i)+xi1/2bi1/2(τ,αi)xβi1(τ)ili(xβi1(τ)ixβi1(τ)i1)ci(τ,αi)), (3.19)
    Ei,i1(αi,τ)=xi1/2bi1/2(τ,αi)xβi1(τ)i1li(xβi1(τ)ixβi1(τ)i1), (3.20)
    E1,1(α1,τ)=x1+1/2b1+1/2(τ,α1)xβ1(τ)1l1(xβ1(τ)2xβ1(τ)1)+14l1x1(a1/2(τ,α1)+b1/2(τ,α1))c1(τ,α1) (3.21)
    E1,2(α1,τ)=x1+1/2b1+1/2(τ,α1)xβ1(τ)2l1(xβ1(τ)2xβ1(τ)1) (3.22)

    G(α,τ)=[14l1x1(a1/2(τ,α1)b1/2(τ,α1))v000xN11/2bN11/2(τ,αN11)xβN11(τ)N1lN11(xβN11(τ)N1xβN11(τ)N11)vN1].

    Theorem 2. Assume that c given in (3.1) is negative and let h=max1iN1li. If h is relatively small then the matrix E(α,τ) in the system (3.17) is an M-matrix for any αAN11.

    Proof. Let us show that E(α,τ) has positive diagonal, non-positive off diagonal, and is diagonally dominant. We first note that

    bi+1/2(τ,α)xβi(τ)i+1xβi(τ)i=ai+1/2(τ,α)βi(τ)xβi(τ)i+1xβi(τ)i>0, (3.23)

    for i=1,,N11, and all bi+1/2(τ,α)0,bi1/2(τ,α)0, with ai+1/2(τ,α)>0 and ai1/2(τ,α)>0.

    This also holds when bi+1/2(τ,α)0 and bi1/2(τ,α)0, that is

    limbi+1/2(τ,α)0bi+1/2(τ,α)xβi(τ)i+1xβi(τ)i=bi+1/2(τ,α)eβi(τ)ln(xi+1)eβi(τ)ln(xi)=bi+1/2(τ,α)βi(τ)ln(xi+1)βi(τ)ln(xi)=ai+1/2(τ,α)(lnxi+1xi)1>0,limbi1/2(τ,α)0bi1/2(τ)xβi1(τ)ixβi1(τ)i1=bi1/2(τ,α)eβi1(τ)ln(xi)eβi1(τ)ln(xi1)=bi1/2(τ,α)βi1(τ)ln(xi)βi1(τ)ln(xi1)=ai1/2(τ,α)(lnxixi1)1>0. (3.24)

    Using the definition of E(α,τ) given above, we see that

    Ei,i0,Ei,i+10,Ei,i10i=2,,N11,
    |Ei,i||Ei,i1|+|Ei,i+1|

    because xβi(τ)i+1xβi(τ)i+xβi(τ)1iβi(τ)h, xβi1(τ)i1xβi1(τ)ixβi1(τ)1iβi1(τ)h and

    |Ei,i||Ei,i+1||Ei,i1|=(bi+1/2(τ)xβi(τ)i+1xβi(τ)i)>0(hβiβixβi1i)00+(bi1/2(τ,α)xβi1(τ)ixβi1(τ)i1)>0(hβi1βi1xβi11i)00ci(τ,αi).

    Note that for i=1, we have E1,10 if a1/2(τ,α1)+b1/2(τ,α1), are nonnegative and c1(τ,α1)<0. So E(α,τ) is diagonally dominant and is therefore an M-matrix.

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

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

    v(t,x,y)t+supαA[(k(t,x,y,α))+c(t,x,y,α)v(t,x,y)]=0, (3.25)

    where k(t,x,y,α)=A(t,x,y,α)v(t,x,y)+bv(t,x,y) is the flux,

    b=(xb1(t,x,y,α),yb2(t,x,y,α))TandA=(a11a12a21a22), (3.26)

    We will assume that a21=a12. We also define the following coefficients, which will help us to build our scheme smoothly a11(t,x,y,α)=a(t,x,y,α)x2,a22(t,x,y,α)=¯a(t,x,y,α)y2 and a12=a21=d1(t,x,y,α)xy. Although this initial value problem (3.25) is defined on the unbounded region R2, we often restrict our consideration to a bounded region. As usual the two dimensional domain is truncated to Ix×Iy, where Ix=[0,xmax] and Iy=[0,ymax] be divided into N1 and N2 sub-intervals:

    Ixi:=(xi,xi+1),Iyj:=(yj,yj+1),i=0,,N11,j=0,,N21

    with 0=x0<x1<<xN1=xmax and 0=y0<y1<<yN2=ymax. This defines a mesh on Ix×Iy with all the mesh lines perpendicular to one of the axes.

    We also set

    xi+1/2=xi+xi+12,xi1/2=xi1+xi2,yj+1/2=yj+yj+12,yj1/2=yj1+yj2,

    for each i=1,,N11 and each j=1,,N21. We denote N=(N11)×(N21). These mid-points form a second partition of Ix×Iy if we define x1/2=x0, xN1+1/2=xmax, y1/2=y0 and yN2+1/2=ymax. For each i=0,1,,N1 and j=0,1,,N2, we set hxi=xi+1/2xi1/2 and hyj=yj+1/2yj1/2.

    We now discuss the finite volume method for (3.25). Integrating both size of (3.25) over Ri,j=[xi1/2,xi+1/2]×[yj1/2,yj+1/2], we have

    xi+1/2xi1/2yj+1/2yj1/2v(t,x,y)tdxdy+xi+1/2xi1/2yj+1/2yj1/2supαA[(k(v(t,x,y,α)))+c(t,x,y,α)v(t,x,y)]dxdy=0, (3.27)

    for i=1,2,,N11, j=1,2,,N21. Applying the mid-points quadrature rule to the first and the last point terms, we obtain the above

    dvi,j(t)dtli,j+supαi,jA[Ri,j(k(v(t,x,y,αi,j))dxdy+ci,j(t,αi,j)vi,j(t)li,j]=0 (3.28)

    for i=1,2,N11, j=1,2,N21 where li,j=(xi+1/2xi1/2)×(yj+1/2yj1/2) is the length of Ri,j, and vi,j(t) denotes the nodal approximation to v(t,xi,yj). We now consider the approximation of the middle term in (3.28). Let n denote the unit vector outward-normal to Ri,j. By Ostrogradski Theorem, integrating by parts and using the definition of flux k(v), we have

    Ri,j(k(v))dxdy=Ri,jk(v(t,x,y,αi,j))nds=(xi+1/2,yj+1/2)(xi+1/2,yj1/2)(a11vx+a12vy+xb1v)dy(xi1/2,yj+1/2)(xi1/2,yj1/2)(a11vx+a12vy+xb1v)dy+(xi+1/2,yj+1/2)(xi1/2,yj+1/2)(a21vx+a22vy+yb2v)dx(xi+1/2,yj1/2)(xi1/2,yj1/2)(a21vx+a22vy+yb2v)dx. (3.29)

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

    (xi+1/2,yj+1/2)(xi+1/2,yj1/2)(a11vx+a12vy+xb1v)dy(a11vx+a12vy+xb1v)|(xi+1/2,yj)hyj. (3.30)

    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 (xi+1/2,yj), of the interval Ixi for i=0,1,N11. This discussion is divided into two cases for i1 and i=0 on I0=[0,x1]. This is really an extension of the one dimensional fitted finite volume presented in the previous section.

    CaseI_: For i1.

    Remember that a11(t,x,y,α)=a(t,x,y,α)x2, we approximate the term (a11vx+xb1v) by solving the following two points boundary value problem

    (a(t,xi+1/2,yj,αi,j)xvx+b1(t,xi+1/2,yj,αi,j)v)=0,v(t,xi,yj)=vi,j(t),v(t,xi+1,yj)=vi+1,j(t). (3.31)

    Integrating (3.31) yields the first-order linear equations

    a(t,xi+1/2,yj,αi,j)xvx+b1(t,xi+1/2,yj,αi,j)v=C1 (3.32)

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

    C1=b1i+1/2,j(t,αi,j)(xβi,j(t)i+1vi+1,jxβi,j(t)ivi,j)xβi,j(t)i+1xβi,j(t)i. (3.33)

    Therefore,

    a11vx+a12vy+xb1vxi+1/2(b1i+1/2,j(t,αi,j)(xβi,j(t)i+1vi+1,jxβi,j(t)ivi,j)xβi,j(t)i+1xβi,j(t)i+d1yvy), (3.34)

    where b1(t,xi+1/2,yj,αi,j)=b1i+1/2,j(t,αi,j), a(t,xi+1/2,yj,αi,j)=ai+1/2,j(t,αi,j), βi,j(t)=b1i+1/2,j(t,αi,j)ai+1/2,j(t,αi,j) and a12=a21=d1(t,x,y,α)xy. Finally, we use the forward difference,

    vy|(xi+1/2,yj)vi,j+1vi,jhyj

    Finally,

    [a11vx+a12vy+xb1v](xi+1/2,yj)hyjxi+1/2×(b1i+1/2,j(t,αi,j)(xβi,j(t)i+1vi+1,jxβi,j(t)ivi,j)xβi,j(t)i+1xβi,j(t)i+d1i,j(t,αi,j)yjvi,j+1vi,jhyj)hyj. (3.35)

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

    [a11vx+a12vy+xb1v](xi1/2,yj)hyjxi1/2×(b1i1/2,j(t,αi,j)(xβi1,j(t)ivi,jxβi1,j(t)i1vi1,j)xβi1,j(t)ixβi1,j(t)i1+d1i,j(t,αi,j)yjvi,j+1vi,jhyj)hyj. (3.36)

    CaseII_: For j1.

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

    (xi+1/2,yj+1/2)(xi1/2,yj+1/2)(a21vx+a22vy+yb2v)dx(a21vx+a22vy+yb2v)|(yj+1/2,xi)hxi. (3.37)

    Following the first case of this section, we have

    [a21vx+a22vx2+yb2v](xi,yj+1/2)hxiyj+1/2×(b2i,j+1/2(t,αi,j)(yj+1ˉβi,j(t)vi,j+1yjˉβi,j(t)vi,j)yj+1ˉβi,j(t)yjˉβi,j(t)+d1i,j(t,αi,j)xivi+1,jvi,jhxi)hxi.

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

    [a21vx+a22vy+yb2v](xi,yj1/2)hxiyj1/2×(b2i,j1/2(t,αi,j)(yˉβi,j1(t)jvi,jyˉβi,j1(t)j1vi,j1)yˉβi,j1(t)jyˉβi,j1(t)j1+d1i,j(t,αi,j)xivi+1,jvi,jhxi)hxi, (3.38)

    for j=1,,N21, where ˉβi,j(t)=b2i,j+1/2(t,αi,j)ˉai,j+1/2(t,αi,j) with a22(t,x,y,α)=ˉa(t,x,y,α)y2.

    CaseIII_: Approximation of the flux at I0. Note that the analysis in case I does not apply to the approximation of the flux on [0,x1] because (3.31) is degenerated. Therefore, we reconsider the following form

    (a(t,x1/2,yj,α1,j)xvx+b1(t,x1/2,y,α1,j)v)C2in[0,x1]v(x0,yj)=v0,j,v(x1,yj)=v1,j, (3.39)

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

    [a11vx+a12vy+xb1v](x1/2,yj)hyjx1/2(12[(ax1/2,j(t,α1,j)+b1x1/2,j(t,α))v1,j(ax1/2,j(t,α1,j)b1x1/2,j(t,α1,j))v0,j]+d11,j(t,α1,j)yjv1,j+1v1,jhyj)hyj. (3.40)

    CaseIV_: Approximation of the flux at J0.

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

    [a21vx+a22vy+yb2v](xi,y1/2)hxiy1/2(12[(ˉai,y1/2(t,αi,1)+b2i,y1/2(t,α))vi,1(ˉai,y1/2(t,αi,1)b2i,y1/2(t,αi,1))vi,0]+d1i,1(t,αi,1)xivi+1,1vi,1hxi)hxi. (3.41)

    By replacing the flux by his value for i=1,,N11 and j=1,,N21, Eq (3.28) becomes

    dvi,jdt+supαi,jA1li,j[xi+1/2(b1i+1/2,j(t,α)(xβi,j(t)i+1vi+1,jxβi,j(t)ivi,j)xβi,j(t)i+1xβi,j(t)i)hyj+xi+1/2(d1i,j(t,αi,j)yjvi,j+1vi,jhyj)hyjxi1/2(b1i1/2,j(t,αi,j)(xβi1,jivi,jxβi1,j(t)i1vi1,j)xβi1,j(t)ixβi1,j(t)i1+d1i,j(t,αi,j)yjvi,j+1vi,jhyj)hyj+yj+1/2(b2i,j+1/2(t,αi,j)(yj+1ˉβi,j(t)vi,j+1yjˉβi,j(t)vi,j)yj+1ˉβi,j(t)yjˉβi,j(t)+d1i,j(t,αi,j)xivi+1,jvi,jhxi)hxiyj1/2(b2i,j1/2(t,αi,j)(yjˉβi,j1(t)vi,jyj1ˉβi,j1(t)vi,j1)yjˉβi,j1(t)yj1ˉβi,j1(t)+d1i,j(t,αi,j)xivi+1,jvi,jhxi)hxi+ci,j(t,αi,j)vi,jli,j]=0 (3.42)

    By setting τ=Tt and including the boundary conditions, we have the following system

    {supαAN[ei,ji1,j(τ,α)vi1,j+ei,ji,j(τ,α)vi,j+ei,ji+1,j(τ,α)vi+1,j+ei,ji,j1(τ,α)vi,j1+ei,ji,j+1(τ,α)vi,j+1]dvi,jdτ=0,withvi,j(0)given,

    where for i=1,,N11, j=1,,N21 and N=(N11)×(N21), AN=A×A××A(N11)×(N21). We have

    e1,j0,j=14l1,jhyjx1(ax1/2,j(τ,α1,j)b1x1/2,j(τ,α1,j))v0,je1,j1,j=14l1,jhyjx1(ax1/2,j(τ,α1,j)+b1x1/2,j(τ,α1,j))12c1,j(τ,α1,j)+d11,j(τ,α1,j)xihyjl1,j+x1+1/2hyjb11+1/2,j(τ,α1,j)xβ1,j(τ)1l1,j(xβ1,j(τ)2xβ1,j(τ)1)e1,j2,j=d11,j(τ,α1,j)xihyjl1,jx1+1/2hyjb11+1/2,j(τ,α1,j)xβ1,j(τ)2l1,j(xβ1,j(τ)2xβ1,j(τ)1)ei,1i,0=14li,1hxiy1(ˉai,y1/2(τ,αi,1)b2i,y1/2(τ,αi,1))vi,0ei,1i,1=14li,1hxiy1(ˉai,y1/2(τ,αi,1)+b2i,y1/2(τ,αi,1)12ci,1(τ,αi,1)+d1i,1(τ,αi,1)yjhxili,1+y1+1/2hxib2i,1+1/2(τ,αi,1)yˉβi,1(τ)1li,1(yˉβi,1(τ)2yˉβi,1(τ)1)ei,1i,2=d1i,1(τ,αi,1)yjhxili,1y1+1/2hxib2i,1+1/2(τ,αi,1)y2ˉβi,1(τ)li,1(y2ˉβi,1(τ)y1ˉβi,1(τ))ei,ji+1,j=d1i,j(τ,αi,j)xihyjli,jxi+1/2hyjb1i+1/2,j(τ,αi,j)xβi,j(τ)i+1li,j(xβi,j(τ)i+1xβi,j(τ)i)ei,ji1,j=xi1/2hyjb1i1/2,j(τ,αi,j)xβi1,j(τ)i1li,j(xβi1,j(τ)ixβi1,j(τ)i1)ei,ji,j+1=d1i,j(τ,αi,j))yjhxili,jyj+1/2hxib2i,j+1/2(τ,αi,j))yj+1ˉβi,j(τ)li,j(yj+1ˉβi,j(τ)yjˉβi,j(τ))ei,ji,j1=yj1/2hxib2i,j1/2(τ,αi,j))yj1ˉβi,j1(τ)li,j(yjˉβi,j1(τ)yj1ˉβi,j1(τ)),ei,ji,j=d1i,j(τ,αi,j))xihyjli,j+xi+1/2hyjb1i+1/2,j(τ,αi,j))xβi,j(τ)ili,j(xβi,j(τ)i+1xβi,j(τ)i)+xi1/2hyjb1i1/2,j(τ,αi,j))xβi1,j(τ)ili,j(xβi1,j(τ)ixβi1,j(τ)i1)ci,j(τ,αi,j))d1i,j(τ,αi,j))yjhxili,j+yj+1/2hxib2i,j+1/2(τ,αi,j))yjˉβi,j(τ)li,j(yj+1ˉβi,j(τ)yjˉβi,j(τ))+yj1/2hxib2i,j1/2(τ,αi,j))yjˉβi,j1(τ)li,j(yjˉβi,j1(τ)yj1ˉβi,j1(τ)). (3.43)

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

    {dv(τ)dτ+infαAN[E(τ,α)v(τ)+F(τ,α)]=0,withv(0)given, (3.44)

    or

    {dv(τ)dτ=supαAN[A(τ,α)v(τ)+G(τ,α)]withv(0)given, (3.45)

    where A(τ,α)=E(τ,α), E(τ,α)(I,J)=(ei,ji,j), i,i=1,,N11, j,j=1,,N21, n1=N11,n2=N21;I:=I(i,j)=i+(j1)n1 and J:=J(i,j)=i+(j1)n1, v=(v1,1,,v1,N21,,vN11,1,,vN11,N21) and G(τ,α)=F(τ,α) includes boundary condition.

    Theorem 3. Assume that the coefficients of A given by (3.26) are positive, c<0 and let h=max1iN11jN2{hxi,hyj}. If h is relatively small then the matrix E(τ,α)=(ei,ji,j)i=1,,N11j=1,,N21, in (3.44) is an M-matrix for any αAN.

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

    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

    dv(τ)dτ=supαAN[A(τ,α)v(τ)+G(τ,α)],v(0)given. (4.1)

    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,m1. We denote v(τn)vn, An(α)=A(τn,α) and Gn(α)=G(τn,α).

    For θ[12,1], following [11], the θ-Method approximation of (4.1) in time is given by

    vn+1vn=ΔtsupαAN(θ[An+1(α)vn+1+Gn+1(α)]+(1θ)[An(α)vn+Gn(α)]),

    this also can be written as

    infαAN([I+ΔtθEn+1]vn+1+Fn+1(α)+[I+ΔtθEn]vn+Fn(α))=0. (4.2)

    As we can notice, to find the unknown vn+1 we need also to solve an optimization problem. Let

    αn+1(argsupαAN{θΔt[An+1(α)vn+1+Gn+1(α)]+(1θ)Δt[An(α)vn+Gn(α)]}). (4.3)

    Then, the unknown vn+1 is solution of the following equation

    [IθΔtAn+1(αn+1)]vn+1=[I+(1θ)ΔtAn(αn+1)]vn+[θΔtGn+1(αn+1)+(1θ)ΔtGn(αn+1)], (4.4)

    Note that when θ=12 the time-stepping scheme becomes the Crank-Nicolson scheme and when θ=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 (vn+1)0=vn,

    2. Let ˆvk=(vn+1)k,

    3. For k=0,1,2 until convergence (ˆvk+1ˆvkϵ, given tolerance) solve

    αki (4.5)
    (argsupαAN{θΔt[An+1(α)ˆvk+Gn+1(α)]i+(1θ)Δt[An(α)vn+Gn(α)]i})αk=(αk)i (4.6)
    [IθΔtAn+1(αk)]ˆvk+1=[I+(1θ)ΔtAn(αk)]vn (4.7)
    +[θΔtGn+1(αk)+(1θ)ΔtGn(αk)],

    4. Let kl being the last iteration in step 3, set vn+1:=ˆvkl, αn+1:=αkl.

    Indeed to find the solution at time τn+1, we use the solution at time τ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+ΔtθEn+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,,m1, the system matrix [I+ΔtθEn+1] in (4.2) is an M–matrix for each αAN.

    Proof. The proof is obvious. Indeed as in Theorem 2 and 3, [I+ΔtθEn+1] is (strictly) diagonally dominant since Δ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]).

    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 ϵ=106 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 α=α(t) is a feedback control belongs in [0,1]

    v(t,x)=maxα[0,1]E{1pxp(T)}. (5.1)

    subject to

    dxt=bαt(t,xt)dt+σαt(t,xt)dWt,

    where bαt(t,xt)=xt[r+αt(μr)], σαt(t,xt)=xtαtσ, 0<p<1 is coefficient of risk aversion, r, μ, σ are constants, xtR and Wt Brownian motion. We assume μ>r. For this problem, using dynamic programming the corresponding HJB equation is given by

    {vt(t,x)+supα[0,1][Lαv(t,x)]=0on [0,T)×R+v(T,x)=xpp,xR+ (5.2)

    where

    Lαv(t,x)=(bα(t,x))v(t,x)x+(aα(t,x))2v(t,x)x2, (5.3)

    and aα(t,x)=12(σα(t,x))2.

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

    v(t,x)t+supα[0,1][x(a(t,x,α)x2v(t,x)x+b(t,x,α)xv(t,x))+c(t,x,α)v(t,x)]=0,where (5.4)
    a(t,x,α)=12σ2α2b(t,x,α)=r+(μr)ασ2α2c(t,x,α)=(r+α(μr)σ2α2).

    The domain where we compare the solution is Ω=[0,xmax], 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 (τn,xi) by

    v(τn,xi)=exp(p×(n×ΔtT)×ρ)×(xi)pp, (5.5)
    ρ=(r+(μr)2σ2(1p)+12(p1)σ2((μr)σ2(1p))2),0<p<1 (5.6)

    We use the following L2([0,T]×Ω) norm of the absolute error

     vmvL2[[0,T]×Ω]=(m1n=0N11i=1(τn+1τn)×li×(vniv(τn,xi,))2)1/2, (5.7)

    where vm is the numerical approximation of v computed from our numerical scheme. The 3 D graphs of the implicit fitted finite volume (θ=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, μ=0.0657, σ=0.2537, p=0.5255 and T=1.

    Figure 1.  Implicit fitted finite volume.
    Figure 2.  Ansatz analytical solution.

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

    Table 1.  Comparison of the implicit fitted finite volume method and implicit finite difference method. We have taken [0,10] for computational domain with N=1500, r=0.0449, μ=0.0657, σ=0.2537, p=0.5255 and T=1. By fitting the data from this table, we found that the convergence orders in time are 0.91658 and 0.91438 respectively for the fitted finite volume method and the finite difference method.
    Time subdivision 200 150 100 50
    Error of fitted finite volume method 3.34 E-01 6.81 E-01 1.01 E-00 1.33 E-00
    Error of finite difference method 3.37 E-01 6.89 E-01 1.02 E-00 1.34 E-00

     | Show Table
    DownLoad: CSV

    From Tables 13, 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.

    Table 2.  Comparison of the implicit fitted finite volume method and implicit finite difference method. We have taken [0,5] for computational domain with N=1500, r=0.0449, μ=0.0657, σ=0.2537, p=0.5255 and T=1.
    Time subdivision 200 150 100 50
    Error of fitted finite volume method 1.64 E-01 3.35 E-01 4.97 E-01 6.53 E-01
    Error of finite difference method 1.66 E-01 3.38 E-01 5.03 E-01 6.61 E-01

     | Show Table
    DownLoad: CSV
    Table 3.  Comparison of the implicit fitted finite volume method and implicit finite difference method. We have taken [0,5] for computational domain with N=2000, r=0.0449, μ=0.0657, σ=0.2537, p=0.5255 and T=1.5. By fitting the data from this table, we found that the convergence orders in time are 0.89522 and 0.8948 respectively for the fitted finite volume method and the finite difference method.
    Time subdivision 200 150 100 50
    Error of fitted finite volume method 2.97 E-01 5.99 E-01 8.83 E-01 1.15 E-01
    Error of finite difference method 3.00 E-01 6.05 E-01 8.91 E-01 1.17 E-01

     | Show Table
    DownLoad: CSV

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

    v(t,x,y)=max(α1,α2)[0,1]×[0,1]E{1pxp(T)×1pyp(T)}, (5.8)
    subjectto{dxt=b1α1t(t,xt)dt+σα1t(t,xt)dW1tdyt=b2α2t(t,yt)dt+σα2t(t,yt)dW2t (5.9)

    where

    b1α1t(t,xt)=xt[r1+α1t(μ1r1)],b2α2t(t,yt)=yt[r2+α2t(μ2r2)],σα1t(t,xt)=xtα1tσ,σα2t(t,yt)=ytα2tσ,

    0<p<1 is coefficient of risk aversion, r1, μ1, r2, μ2σ are constants, xt,ytR and ρ12[0,1) the correlation of the two Brownian motion. We assume that μ1>r1 and μ2>r2. For this problem, using dynamic programming the corresponding HJB equation is given by

    {vt(t,x,y)+sup(α1,α2)[0,1]×[0,1][Lα1,α2v(t,x,y)]=0on [0,T)×R+×R+v(T,x,y)=xpp×ypp,x,yR+ (5.10)

    where

    Lαv(t,x,y)=(b1α1(t,x))v(t,x,y)x+(b2α2(t,y))v(t,x,y)y+12(σα1(t,x))22v(t,x,y)x2+12(σα2(t,y))22v(t,x,y)y2+(σα1(t,x))(σα2(t,y))ρ122v(t,x,y)xy,

    and the two dimensional divergence form is given by

    v(t,x,y)t+supα[0,1]×[0,1][(k(t,x,y,α))+c(t,x,y,α)v(t,x,y)]=0, (5.11)

    where k(t,x,y,α)=A(t,x,y,α)v(t,x,y)+b(t,x,y,α)v(t,x,y)

    A=(a11a12a21a22),
    a11(t,x,y,α)=12σ2α21x2,a22(t,x,y,α)=12σ2α22y2,a12(t,x,y,α)=a21(x,y,t,α)=12σ2α1α2ρ12xy.

    By identification,

    a(t,x,y,α)=12σ2α21ˉa(t,x,y,α)=12σ2α22b1(t,x,y,α)=r1+α1(μ1r1)12σ2α1α2ρ12σ2α21,b2(t,x,y,α)=r2+α2(μ2r2)12σ2α1α2ρ12σ2α22,c(t,x,y,α)=[r1+(μ1r1)α1][r2+(μ2r2)α2]+σ2(α21+α22+α1α2ρ12),d1(t,x,y,α)=12σ2α1α2ρ12.

    The two dimensional Ansatz exact solution [9] at (τn,xi,yj) is given by

    v(τn,xi,yj)=exp(p×(n×ΔtT)×ρ)×(xi)pp×(yj)pp,ρ=supα1,α2[0,1]×[0,1][r1+r2+(μ1r1)α1+(μ2r2)α2+12σ2α21(p1)+12σ2α22(p1)+σ2α1α2ρ12p],0<p<1/2.

    We use the following L2([0,T]×Ω),Ω=[0,xmax]×[0,ymax] norm of the absolute error

    vmvL2[Ω×[0,T]]=(m1n=0N11i=1N21j=1(τn+1τn)×hxi×hyj×(vni,jv(τn,xi,xj,))2)1/2, (5.12)

    where vm 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 are given at Figures 3 and 4, with N1=25, N2=20, r1=0.0449/2, μ1=0.0657/2, r2=0.0448/2, μ2=0.0656/2, σ=0.2537/2, p=0.5255/2 and ρ12=0.9.

    Figure 3.  Implicit fitted finite volume at finite time T = 1.
    Figure 4.  Ansatz analytical solution at finite time T = 1.

    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 ρ12=0.9 the fitted finite volume and the finite difference converge with order 1.3 in time.

    Table 4.  Errors table for fitted finite volume method and finite difference method in dimension 2 with ρ12=0.9. By fitting the data from this table, we found that the convergence order in time is 1.3 for both the fitted finite volume method and the finite difference method.
    Time subdivision 200 150 100 50
    Error of fitted finite volume method 2.90 E-03 6.50 E-03 1.15 E-02 2.09 E-02
    Error of finite difference method 3.01 E-03 6.52 E-02 1.16 E-02 2.10 E-01

     | Show Table
    DownLoad: CSV

    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 wt at time t, and the proportion of wealth invested in the stock at time t is ut. The interest rate earned in the bank account is r1, the return from the stock at time t has two components, the cash dividend rate r2, the capital gain rate Rt. The dynamic of the capital gain rate Rt is assumed to be governed by the stochastic process

    dRt=[β1Rt+f]dt+σtdW1t, (5.13)

    and the volatility σt with modeled by

    dσt=ασtdt+βσtdW2t. (5.14)

    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 ut[0,1] and the wealth dynamics for this cash management problem is given by

    dwt=wtutr2dt+wt(1ut)r1dt+wtRtdtd(t)wtdt. (5.15)

    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

    J(T,w,R,σ)=maxu[0,1]E{wT}. (5.16)
    subjectto{dwt=wtutr2dt+wt(1ut)r1dt+wtRtdtd(t)wtdt,dRt=[β1Rt+f]dt+σtdW1tdσt=ασtdt+βσtdW2t

    We assume that the two Brownian motions are correlated, that is dW1tdW2t=ρ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

    Jt+maxu[0,1]{(f+β1R)JR+w(ur2+(1u)r1+uRd(t)Jw+1/2(σ2JRR+β2σ2Jσσ+2ρβσ2JσR)+ασJσ}=0, (5.17)

    with terminal condition J(T,)=wT. To simplify the problem, we assume that

    J(t,w,R,σ)=wH(t,R,σ).

    Therefore (5.17) is equivalent to

    Ht+maxu[0,1]{(f+β1R)HR+(ur2+(1u)r1+uRd(t))H+1/2(σ2JRR+β2σ2Hσσ+2ρβσ2HσR)+(ασ)Hσ}=0 (5.18)

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

    H(t,R,σ)t+supu[0,1][(k(t,R,σ,u))+c(t,R,σ,u)H(t,R,σ)]=0, (5.19)

    where k(t,R,σ,u)=A(t,R,σ,u)H(t,R,σ)+b(t,R,σ)H(t,R,σ)

    A=(a11a12a21a22),
    a11=12σ2,a22=12β2σ2,a12=a21=12σ2ρβ.

    By identification,

    a(t,R,σ)=σ22R2,ˉa(t,R,σ)=β22b1(t,R,σ)=fR+β1σρβR,b2(t,R,σ)=αβ2c(t,R,σ,u)=ur2+(1u)r1+uRd(t)β1α+β2.

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

    H(t,0,σ)=0,H(t,R,σmax)=R,HR(t,Rmax,σ)=1. (5.20)

    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 σ=0 it is considered inserting σ=0 into the PDE to complete the set of four boundary conditions:

    Ht(t,R,0)+maxu[0,1]{(f+β1R)HR(t,R,0)+(ur2+(1u)r1+uRd)H(R,0,t)}=0 (5.21)

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

    Table 5.  Parameters values for the cash management problem.
    f β1 β α r1 r2 ρ
    0.12 0.96 0.3 0.85 0.024 0.01 0.5

     | Show Table
    DownLoad: CSV

    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 N1=10, N2=10, Rmax=1/2, σ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.

    Figure 5.  A fitted finite volume sample solution of the wealth rate H at the point 1/2,1/2.
    Table 6.  Optimal choice using fitted scheme and finite difference.
    Time 1 2 3 4 5 6 7 8 9 10
    Optimal choice using fitted scheme 1 1 1 1 1 1 1 1 1 1
    Optimal choice using finite difference 1 0 1 1 1 1 0 0 1 1

     | Show Table
    DownLoad: CSV

    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.

    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).

    The authors declare no potential conflict of interests.



    [1] D. N. DE G. Allen, R. V. Southwell, Relaxation methods applied to determine the motion, in two dimensions, of a viscous fluid past a fixed cylinder, Q. J. Mech. Appl. Math., 8 (1955), 129–145. doi: 10.1093/qjmam/8.2.129
    [2] C. Benezet, J. Chassagneux, C. Reisinger, A numerical scheme for the quantile hedging problem. Available from: https://arXiv.org/abs/1902.11228v1.
    [3] V. Henderson, K. Kladvko, M. Monoyios, C. Reisinger, Executive stock option exercise with full and partial information on a drift change point. Available from: arXiv: 1709.10141v4.
    [4] C. S. Huang, C. H. Hung, S. Wang, A fitted finite volume method for the valuation of options on assets with stochastic volatilities, Computing, 77 (2006), 297–320. doi: 10.1007/s00607-006-0164-4
    [5] N. V. Krylov, Controlled Diffusion Processes, Applications of Mathematics, Springer-Verlag, New York, 1980.
    [6] S. Wang, A Novel fitted finite volume method for the Black-Scholes equation governing option pricing, IMA J. Numer. Anal., 24 (2004), 699–720. doi: 10.1093/imanum/24.4.699
    [7] P. Forsyth, G. Labahn, Numerical methods for controlled Hamilton-Jacobi-Bellman PDEs in finance, J. Comput. Finance, 11 (2007), 1–43.
    [8] P. Wilmott, The Best of Wilmott 1: Incorporating the Quantitative Finance Review, John Wiley & Sons, 2005.
    [9] H. Pham, Optimisation et contrôle stochastique appliqués à la finance, Mathématiques et applications, Springer-Verlag, New York, 2000.
    [10] L. Angermann, S. Wang, Convergence of a fitted finite volume method for the penalized Black–Scholes equation governing European and American Option pricing, Numerische Math., 106 (2007), 1–40. doi: 10.1007/s00211-006-0057-7
    [11] H. Peyrl, F. Herzog, H. P. Geering, Numerical solution of the Hamilton-Jacobi-Bellman equation for stochastic optimal control problems, WSEAS international conference on Dynamical systems and control, 2005.
    [12] N. V. Krylov, On the rate of convergence of finite-difference approximations for Bellman's equations with variable coefficients, Probab. Theory Relat. Fields, 117 (2000), 1–16. doi: 10.1007/s004400050264
    [13] N. V. Krylov, The rate of convergence of finite-difference approximations for Bellman equations with Lipschitz coefficients, Appl. Math. Optim., 52 (2005), 365–399. doi: 10.1007/s00245-005-0832-3
    [14] E. R. Jakobsen, On the rate of convergence of approximations schemes for Bellman equations associated with optimal stopping time problems, Math. Models Methods Appl. Sci., 13 (2003), 613–644. doi: 10.1142/S0218202503002660
    [15] M. G. Crandall, P. L. Lions, Viscosity solutions of Hamilton-Jacobi equations, Trans. Am. Math. Soc., 277 (1983), 1–42. doi: 10.1090/S0002-9947-1983-0690039-8
    [16] J. Holth, Merton's portfolio problem, constant fraction investment strategy and frequency of portfolio rebalancing, Master Thesis, University of Oslo, 2011. Available from: http://hdl.handle.net/10852/10798.
    [17] Iain Smears, Hamilton-Jacobi-Bellman Equations, Analysis and Numerical Analysis, Research report, 2011. Available from: http://fourier.dur.ac.uk/Ug/projects/highlights/PR4/Smears_HJB_report.pdf
    [18] N. Krylov, Approximating value functions for controlled degenerate diffusion processes by using piece-wise constant policies, Electron. J. Probab., 4 (1999), 1–19. doi: 10.1214/ECP.v4-999
    [19] M. G. Crandall, H. Ishii, P. L. Lions, User's guide to viscosity solutions of second order partial differential equations, Bull. Am. Math. Soc., 27 (1992), 1–67. doi: 10.1090/S0273-0979-1992-00266-5
    [20] C. S. Huang, S. Wang, K. L. Teo, On application of an alternating direction method to Hamilton-Jacobi-Bellman equations, J. Comput. Appl. Math., 27 (2004), 153–166.
    [21] N. Song, W. K. Ching, T. K. Siu, C. K. F. Yiu, On optimal cash management under a stochastic volatility model, East Asian J. Appl. Math., 3 (2013), 81–92. doi: 10.4208/eajam.070313.220413a
    [22] K. Debrabant, Espen R. Jakobsen, Semi-Lagrangian schemes for linear and fully non-linear diffusion equations, Math. Comput., 82 (2013), 1433–1462.
    [23] G. Barles, P. E. Souganidis, Convergence of approximation schemesfor fully nonlinear second order equations, Asymptot. Anal., 4 (1991), 271–283. doi: 10.3233/ASY-1991-4305
  • This article has been cited by:

    1. Mohammad Shirazian, Optimal control design for linear time-varying systems by interpolated variational iteration method, 2022, 1077-5463, 107754632211277, 10.1177/10775463221127739
    2. Christelle Dleuna Nyoumbi, Antoine Tambue, A Novel High Dimensional Fitted Scheme for Stochastic Optimal Control Problems, 2023, 61, 0927-7099, 1, 10.1007/s10614-021-10197-4
    3. Zifu Fan, Youpeng Tao, Wei Zhang, Kexin Fan, Jiaojiao Cheng, Research on open and shared data from government-enterprise cooperation based on a stochastic differential game, 2023, 8, 2473-6988, 4726, 10.3934/math.2023234
    4. Jiamian Lin, Xi Li, SingRu (Celine) Hoe, Zhongfeng Yan, A Generalized Finite Difference Method for Solving Hamilton–Jacobi–Bellman Equations in Optimal Investment, 2023, 11, 2227-7390, 2346, 10.3390/math11102346
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2992) PDF downloads(172) Cited by(4)

Figures and Tables

Figures(5)  /  Tables(6)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog