A Generalized Finite Difference Method for Solving Hamilton–Jacobi–Bellman Equations in Optimal Investment

: This paper studies the numerical algorithm of stochastic control problems in investment optimization. Investors choose the optimal investment to maximize the expected return under uncertainty. The optimality condition, the Hamilton–Jacobi–Bellman (HJB) equation, satisﬁed by the value function and obtained by the dynamic programming method, is a partial differential equation coupled with optimization. One of the major computational difﬁculties is the irregular boundary conditions presented in the HJB equation. In this paper, two mesh-free algorithms are proposed to solve two different cases of HJB equations with regular and irregular boundary conditions. The model of optimal investment under uncertainty developed by Abel is used to study the efﬁcacy of the proposed algorithms. Extensive numerical studies are conducted to test the impact of the key parameters on the numerical efﬁcacy. By comparing the numerical solution with the exact solution, the proposed numerical algorithms are validated.


Introduction
The stochastic optimal control problem of investment optimization, which involves making investment strategies under uncertainty, has always been an important issue in finance and economics.In 1970, Merton [1] conducted a pioneering study on stochastic optimal control in finance, which garnered the application of optimal control in finance widespread attention.Similar theoretical and applied analysis in capital investment includes Abel [2], Karatzas [3], Dixit and Pindyck [4], and Zhou and Li [5].
Stochastic optimal control problems are usually difficult to solve directly in terms of control and state variables.Therefore, most existing studies use indirect methods which are based on the optimality conditions derived from the problem.The indirect methods include the dynamic programming [6] and the maximum principle [7][8][9].The dynamic programming method derives the Hamilton-Jacobi-Bellman (HJB) equation satisfied by the value function that maximizes the objective function.The optimal control problem is then converted into a problem solving the partial differential equations (PDEs) coupled with optimization.The resulting PDEs are often nonlinear, and thus remain challenging to solve.By solving the corresponding HJB equations, the optimal control and value function of the control problems are derived.This approach is the most commonly used solution technique in finance [10][11][12].
The classical numerical algorithms for HJB equations are grid methods (mesh methods) such as the finite difference method (FDM).Wang et al. [13] approximated the viscosity solution of the HJB equation using the upwind finite difference.Peyrl et al. [14] proposed a successive approximation algorithm, where the partial differential equations are solved with the upwind difference scheme.Ma and Ma [15] developed the iterative finite difference methods with policy iterations to solved the sequence of decoupled HJB equations.Inoue et al. [16] investigated the convergence properties of the upwind difference scheme for HJB equations.Forsyth et al. [17] studied the optimal control problem in nonlinear option pricing, using Newton's iteration method to solve the implicit discretized control equations.Wang et al. [12] used a fully implicit time-stepping finite difference method to solve a nonlinear HJB equation derived from an asset allocation problem based on the mean-variance approach.In addition, there are the finite element method [18][19][20][21][22][23] and the finite volume method [24][25][26].
Mesh methods have been widely studied and applied in partial differential equations.However, mesh methods require a regular mesh of grids and are only applicable to regular boundaries, whereas, in dynamic stochastic economic models, the natural modeling domain may have an irregular shape [27,28].For example, in the option pricing problem, the pricing model for American backdated options can be described as a two-dimensional free boundary problem [29][30][31].De Angelis et al. [32] considered the problem of optimal entry into an irreversible investment scheme with a cost function that is non-convex for the control variables, where it is necessary to study optimal entry policies with an irregular boundary.If the state space is irregular, the special handling of region interiors and boundaries is required, which may increase some computational complexity [33].
As a consequence, mesh-free numerical methods have been proposed.Relevant research on meshless or semi-meshless methods for solving partial differential equations includes the global radial basis function method [34], the least squares configuration radial basis function method [35], the Haar wavelet collocation method [36], the Chebyshev method [37], the generalized finite difference method (hereafter GFDM) (see Benito et al. [38,39]), etc.Among them, the GFDM is one of the most promising methods.Based on Taylor series expansion and moving least squares approximation, the GFDM approximates the spatial derivatives at each discrete node as the linear combination of node values in its neighborhood.It overcomes the dependence of the traditional difference method on a regular grid and allows more flexibility in obtaining approximations for functions.Existing research work, such as nonlinear parabolic and hyperbolic partial differential equations studied by Benito et al. [40] and Ureña et al. [41], and linear elliptic partial differential equations studied by Gavete et al. [39], showed that the GFDM is very effective in solving second-order partial differential equations.
Our major contribution is to solve the HJB equation using a numerical scheme based on a meshless method, the GFDM, which is suitable for handling both regular and irregular boundaries.We propose two numerical algorithms, one for the general cases and the other for the case where the optimal control in the HJB equation can be explicitly obtained.For general cases, we propose a successive approximation algorithm which combines the GFDM discretization scheme with the optimization algorithm.For special cases with the explicit expression of optimal control, we first obtain the discretization scheme using the GFDM and then solve the discretized nonlinear equations using Newton's iterative method.We carry out experiments under regular and irregular boundaries, and verify the validity of the proposed numerical algorithms.We also test the effect of key parameters on the accuracy of the proposed algorithms including total points, weighting functions, and irregular index.
The remainder of this paper is organized as follows.Section 2 presents a review of the HJB equation and the GFDM.Section 3 develops the numerical methods to solve the HJB equation.In Section 4, an example of the optimal investment problem is applied to demonstrate the efficiency and accuracy of the numerical methods proposed.Moreover, the key parameters of the algorithms are also tested to verify the validity of the algorithms.Section 5 provides concluding remarks.

Review of Stochastic Optimal Control Problem and HJB Equation
In this section, we will give a review of stochastic optimal control problems over an infinite time horizon.Let us consider the state variable x(s) which satisfies the following stochastic differential equation together with the infinite horizon objective functional The goal of the optimal control problem is to maximize J(t, x, u(.)), that is where we assume that u is the control and the set of admissible controls has the form with U = L 2 (0, ∞; R m ), and U ad ⊂ R m denotes a compact convex subset.The known constant α is a discount rate.w(t) is a k-dimensional standard Brownian motion on a probability space (Ω, F , P), g(x, u) : G × U → R n and σ(x, u) : G × U → R n×k denote the drift and diffusion terms, respectively, which satisfy the linear growth condition and Lipschitz condition such that the stochastic differential equation for the state variable x(t)∈ G has a unique strong solution [42], where G is some open and bounded set on R n .We define the value function of the problem as follows: and set V(x) := V(0, x).This value function gives the best value for every initial condition, given the set of admissible control.Define the Hamiltonian as It is characterized as the solution of the Hamilton-Jacobi-Bellman (HJB) equation where the operator A is defined as and a = σ(x, u)σ(x, u) T is a symmetric matrix.
The HJB Equation (6) includes two subproblems to be solved: the Hamiltonian operator to be optimized and the partial differential equation.

Review of Generalized Finite Difference Method
In this paper, the GFDM, mainly based on the Taylor series expansion and the weighted least squares method, is used to approximate the spatial derivatives.In order to derive the GFD scheme for the HJB Equation (6), we first review the basic principles of GFDM in Benito et al. [38,39].
Assuming that the solution domain for an unknown function V is Z ⊂ R 2 , we need to discretize the solution domain Z into N points, namely Z = {x 1 , x 2 , . . ., x N }.For each x c = (y c , K c ) ∈ Z, we define the subdomain M = {x c ; x 1 , x 2 , . . ., x p } ⊂ Z where the center point x c ∈ Z, and x i = (y i , K i ) ∈ Z where i = 1, . . ., p is a set of points located in the neighborhood of x c selected based on certain criteria such as a four-quadrant method or distance criterion.Define V c = V(x c ) and V i = V(x i ).By Taylor series expansion, V i expanded at the center point x c is where We only retain the Taylor expansion formula to the second-order term, so that the approximate value with second-order accuracy of V i can be obtained, and the approximate value of V i truncated at the second order is set to v i .Based on the formula (8), the weighted residual function can be defined as where w i = w(x i ) denotes the weighting function at point x i .The weighting functions are usually in the form: e −n(dist) 2 or 1 dist n , where dist is the distance between x i and x c , n ∈ N. We define The goal is to minimize the weighted residual function (9) with respect to the partial derivatives.Taking the partial derivatives of function (9) with respect to (10) and setting them equal to 0, i.e., ∂F ∂H = 0, we have ).
From (11), each equation is linear with respect to the five unknown partial derivatives.In order to facilitate the calculation, the above formulae can be rewritten as a system of matrix equations AH = b, where Let . Assuming that all the points in each support domain are different, then it is easy to prove that the solution of linear equations AH = b is unique and H can be explicitly solved as By introducing canonical base vectors e i (i = 1, 2, . . ., p), vector H can be expressed as where We can write the explicit expressions for the value of the five unknown partial derivatives as From (16), it can be seen that the spatial partial derivatives of each point can be approximated as a linear combination of the function values of each node in its support domain, which is similar to the traditional finite difference method.

Numerical Scheme of HJB Equations
In this section, we present the numerical methods of the HJB equation with the Dirichlet boundary condition arising from the corresponding infinite horizon optimal control problem We propose two numerical algorithms for different cases of HJB equations: 1.
For the general case of the HJB equation coupled with optimization, we propose a successive approximation algorithm which combines the GFDM discretization scheme with the optimization algorithm.2.
For the special case where one can explicitly express the optimal control u by maximizing the Hamiltonian, we propose an algorithm combining the GFDM and Newton's iterative method.

The General Case
Assuming that the state variable is x = (y, K), we consider the HJB Equation ( 17) in two-dimensional space, which can be written in the general form We will give a discretization scheme which is based on the GFDM for the spatial variables in the HJB Equation (18).Assuming that the number of interior points in the solution domain is M, then for c = 1, 2, . . ., M, we can obtain the discrete equation at each point x c as with the boundary value V j at the boundary points x j for j = M + 1, M + 2, . . ., N, Using the formulae (16), and denoting the approximate value of V c by v c (c = 1, 2, . . ., M), we have Denote the left-hand side of Equation ( 21) as F c (v), where v = (v 1 , v 2 , . . ., v M ) T .Then, each discrete point x c corresponds to a discrete equation The Hamiltonian in Equation ( 21) can be solved using the standard optimization tools such as "fminbnd" in Matlab software which is based on a golden section search and parabolic interpolation method.After solving the optimization for every discrete point, it turns out that the equations are only related to unknown discrete value functions.Therefore, we need to solve a system of equations consisting of M discrete linear algebra equations with M unknown variables v c , that is Using the GFDM discretization scheme and the optimization method given above, we propose an algorithm alternately updating the control u and the value function v. Algorithm 1 summarizes the general successive approximation algorithm which is an iterative method.

Algorithm 1 Successive Approximation Algorithm
Input: Initial control law ũ0 ; given tolerance ε Output: Approximation of control law u, value function v Solve the boundary value problem using GFDM scheme in Section 3.1 for the given control u k c , i.e., solve equations Compute u k+1 by solving the optimization for every discrete point x c :

The Special Case
For some special cases, one can explicitly express the optimal control u by maximizing the Hamiltonian in the form û = û(x, DV, D 2 V).
Substituting ( 24) into (18), we have for x ∈ Z.That is, the HJB equation is transformed into a partial differential equation without the optimization.Consequently, we can obtain the optimal control and optimal value by directly solving the partial differential equation.However, quite often, the obtained equations are nonlinear partial differential equations, presenting computational challenges.We will use Newton's iteration method to solve the system of nonlinear generalized finite difference equations.Denote the left-hand side of Equation ( 25) by F c (v). Similarly to the discretization scheme in the general case, we can obtain the discrete equations for each discrete point The idea of the iterative method is to treat the system as the zero of the operator Define D as the differential operator.The differentiation of F is then denoted by DF .Algorithm 2 summarizes how to find the zero of F using Newton's iterative method.

Algorithm 2 Newton Iterations
Input: Initial guess ṽ0 ; max number of iterations maxIter; given tolerance ε Output: Approximation of v Initialize k = 0; v 0 = ṽ0 while k < maxIter do Obtain F c (v) for every discrete point using Formula ( 25) and ( 16)

Optimal Investment Problem
The case study presents an optimal investment model under uncertainty in Abel [2].Consider the problem of a firm making optimal investment decisions under uncertainty.The evolution of the firm's capital stock K(t) over time satisfies the following differential equation where u(t) is the control variable denoting the amount of capital investment, and δ ≥ 0 denotes the depreciation rate of the capital stock.The output price y(t) satisfies the following stochastic process where w(t) is a standard Brownian motion on a probability space (Ω, F , P), and the coefficients µ and σ positively represent the growth rate and volatility, respectively.The output technology is assumed to be of the Cobb-Douglas form.The elements of production considered are the capital stock K(t) and the labor L(t).Assuming returns to scale are constant, then the form of the Cobb-Douglas production function is The firm pays a fixed wage w, and the labor can be instantaneously at no cost, so the firm chooses L(t) to maximize the instantaneous operating income at time t.Take the partial derivative of formula (30) with respect to L(t), and set the partial derivative to be 0; it then yields Substituting formula (31) into (30), the maximized total profit at time t becomes where c and γ are constants, and γ = 1 1−ξ , c = γ −γ (γ − 1) γ−1 w 1−γ .The firm will incur certain adjustment costs in its operation.Assuming that the convex adjustment cost generated by capital is h(u(t)), the actual profit at time t is The firm's goal is to choose the capital investment u(t) that maximizes the expected discount cash flow at a given discount rate α over an infinite time horizon.This problem is a standard stochastic optimal control problem, and the value function is where U represents the set of all admissible controls.
Assume that the convex adjustment cost is h(u(t)) = u(t) 2 2 .Using the dynamic programming principle and Ito's formula, we can obtain that the value function V satisfies the following HJB equation In order to solve the Equation (35) using numerical methods, a Dirichlet condition is imposed at the boundary ∂Z: V(y, K) = B(y, K).Then, the problem to be solved is where with

Numerical Results
From (35), we are able to obtain the explicit expression of the optimal control given by û = ∂V ∂K .Therefore, we are able to use both algorithms proposed in Section 3 to obtain the numerical approximation of V and u in problem (36).The effectiveness of the proposed schemes is verified in the cases of regular and node distribution, and the influence of related parameters on numerical accuracy is explored.
We first set the basic parameters involved in the problem (34) as follows: The errors are evaluated using the following formulae: where N represents the total number of points considered in the solution domain, approx(i) and exact(i) represent the approximate solution and exact solution at point x i , respectively, and exact max represents the maximum value of the exact values of the points of all support domains.The analytical solution of Equation ( 35) is The threshold ε used to terminate the iteration is set as 10 −4 , and the maximum number of iterations maxIter is set as 20.All examples in this paper are implemented using Matlab R2017b on a 2.10 GHz desktop computer with 12 GB of memory.
The default weighting function is chosen as e − dist 2 , and the points that constitute the support domain are selected based on the distance criterion (see Figure 1).According to [43], the degree of the irregularity of node clouds will also have an impact on the accuracy of the results.Therefore, under the irregular node distribution, the influence of the irregular index of the cloud of nodes on the numerical solution is mainly considered.Following [43], we define the irregularity index of the clouds of nodes (IIC) as follows.The distance from the node x i (i = 1, 2, . . ., p) to the center node x k in each support domain is Taking each point x k in the solution domain as the central node, the average distance of the support domain where it is located is where p represents the number of nodes in each support domain except the central node.
Then, the average distance from all support domains in the entire solution domain is where M represents the number of interior points in the solution domain.Thus, the irregularity index of clouds of nodes can be defined as To discuss the convergence of the proposed algorithms, we evaluate the order of convergence using the following formula: where and h k denote the exact value of the function V, approximate solution of the function V, the L 2 norm of the error, and the step size of the k-th set of points, respectively.

The Regular Node Distributions
Consider the regular node distribution given in Figure 2. We set the solution domain as a half-unit side square area [0.5, 1] × [1.5, 2], and the total number of nodes N are set as 64, 225, 361 and 529, respectively.The errors of the value function V and the optimal control u are listed in Table 1.It can be seen that our proposed algorithms perform well over the regular domain, and the greater the number of nodes, the smaller the errors.Since we use the exact form of optimal control in conducting the Algorithm 2, it is not surprising that the errors of u obtained by Algorithm 2 are smaller than those obtained by Algorithm 1.The CPU time in seconds and the order of convergence of V are listed in Table 2.The results show that our algorithms are stable and are close to the second-order convergence.
In addition, we implement the standard FDM [44] to solve our problem, where the central-difference approximation for the first-and the second-order derivatives are adopted, and then , the results obtained by the GFDM and by the FDM are compared in Table 3.It is observed that the results from the GFDM has a greater level of accuracy than the results obtained by the FDM.
We then investigate the stability property of proposed algorithms.Table 4 gives the global error and the mean relative error under different numbers of nodes and weighting functions with the implementation of Algorithm 1.It shows the consistent accuracy across various N and weighting functions.The standard deviation of the global errors of the value function V and that of the optimal control u are 5.2468 × 10 −6 and 1.2410 × 10 −3 , respectively.Thus, Algorithm 1 is stable in the numerical approximation.Since Algorithm 2 differs from Algorithm 1 in replacing the optimization search with the Newton iteration due to the known expression in the optimal control, the stability of Algorithm 2 is inferred.Therefore, both algorithms are efficient and robust.
The convergence histories using Algorithm 1 are shown in Figure 3. Figures 4 and 5 show the comparison of the exact solution and numerical solutions of V and u, respectively, with N = 529.In Figure 3, we see that Algorithm 1 converges to the solution of the problem after a few iterations, which demonstrates the fast convergence of the algorithm.
From Table 1 as well as Figures 4 and 5, we observe that Algorithms 1 and 2 perform equally well, so we can conclude that both algorithms can be used.However, we do find that Algorithm 1 takes a longer computational time than Algorithm 2, because Algorithm 1 requires solving the optimization in each iteration as shown in Table 2.

The Irregular Node Distributions
We conduct experiments under two irregular node distributions with the total points N = 361 (we view the rectangular boundaries used by traditional mesh algorithms as regular boundaries, and those with non-rectangular boundaries as irregular ones).The shapes of the nodes are shown in Figure 6.The solution domains of both clouds of nodes are given as 1.
Shape a: The accuracy of the numerical solutions under different shapes is shown in Table 5.The L 2 difference between the two successive iterations obtained by Algorithm 1 is shown in Figure 7. Similarly to the case of regular nodes, Algorithm 1 converges to the solution of the problem after several iterations.We will then test the influence of the irregular index of the clouds of points (IIC) and a key parameter, the weighting functions, on the accuracy of the numerical solution using Algorithm 2.
For each support domain, p = 15 nodes were selected.Table 6 shows the numerical results for different values of IIC.It can be seen that the IIC obtained for the circular region is greater, and the error of the results is larger.From the basic principle of GFDM, it is known that the approximated value of the central node of each support domain is the weighted sum of the function values of the rest of the nodes, and the nodes closer to the central node have a greater impact on the function approximation, so the nodes that are closer to the central node should be given higher weights.Based on this feature, this paper uses two different weighting functions: the exponential function e −n (dist) 2 and the potential function 1  dist n for comparison.Table 7 gives the numerical results under different weighting functions.Both forms of weighting functions, the exponential function and the potential function, yield small errors.Moreover, the exponential weighting function produces more stable results with respect to the changes in n than the potential function.In sum, our two mesh-free numerical methods are robust in solving HJB equations (both general and some special form) with regular and irregular boundary conditions.The effectiveness and the accuracy of the proposed methods are evidenced by comprehensive numerical analyses.

Conclusions
In this paper, we present two algorithms for computing different cases of HJB equations with boundary conditions.The first algorithm, combining the mesh-free GFD method and successive approximation method, is proposed to solve the general case of HJB equations.The second algorithm, combining the mesh-free GFD method and Newton's iterative method, is proposed to solve some special case of HJB equations.Numerical examples show that both algorithms are efficient and have high accuracy over regular and irregular domains.We also tested the effect of key parameters on the accuracy of the proposed algorithm under regular node distributions and irregular node distributions, respectively, including the total points, weighting functions and irregular index.
Due to the meshless nature of the generalized finite difference method, the proposed algorithms can be easily applied for solving nonlinear PDEs over complicated and realistic boundaries.Additionally, it can also be extended to solve HJB-FP systems in the mean field type control problem, which will be investigated in our future study.

Figure 1 .
Figure 1.Criterion for selecting points in the support domain-distance criterion.

Figure 3 .
Figure 3. L 2 difference between two successive iterations with 361 and 529 regular points using Algorithm 1.

Figure 4 .Figure 5 .
Figure 4.The exact solution and numerical solution of value function V (N = 529).

Figure 7 .
Figure 7. L 2 difference between two successive iterations with irregular points using Algorithm 1.

Figures 8 -
show the exact solutions and numerical solutions of V and u under different shapes, respectively.It can be seen that the results are also accurate under irregular points as in regular ones.

Figure 8 .
Figure 8.The exact solution and numerical solution of the value function V under shape a.

Figure 9 .
Figure 9.The exact solution and numerical solution of the optimal control u under shape a.

Figure 10 .
Figure 10.The exact solution and numerical solution of the value function V under shape b.

Figure 11 .
Figure 11.The exact solution and numerical solution of the optimal control u under shape b.

Table 1 .
Numerical results of both algorithms under regular node distributions.

Table 2 .
Comparison of both algorithms on the order of convergence and computation time (s).

Table 3 .
Comparison of the numerical errors with the FDM and the GFDM using Algorithm 1.

Table 4 .
Robustness analysis of Algorithm 1 for the number of nodes and weighting functions.

Table 5 .
Numerical errors of both algorithms under irregular node distributions.

Table 6 .
Numerical errors of Algorithm 2 under irregular nodes versus IIC.