Finite Difference Methods for the Hamilton–Jacobi–Bellman Equations Arising in Regime Switching Utility Maximization

For solving the regime switching utility maximization, Fu et al. (Eur J Oper Res 233:184–192, 2014) derive a framework that reduce the coupled Hamilton–Jacobi–Bellman (HJB) equations into a sequence of decoupled HJB equations through introducing a functional operator. The aim of this paper is to develop the iterative finite difference methods (FDMs) with iteration policy to the sequence of decoupled HJB equations derived by Fu et al. (2014). The convergence of the approach is proved and in the proof a number of difficulties are overcome, which are caused by the errors from the iterative FDMs and the policy iterations. Numerical comparisons are made to show that it takes less time to solve the sequence of decoupled HJB equations than the coupled ones.


Introduction
The utility maximization is a kind of stochastic control problems. The dynamic programming approach is often applied to the optimal value function and the so-called HJB equation is derived (see the books [25,34] for the stochastic control and its applications). Since the HJB equation is a fully nonlinear PDE, the closed-form classical solution cannot be found except for some simple cases: a Black-Scholes complete market model with particular utility functions, see [6,7]. For constrained market models it has to use numerical methods to solve  the HJB equations. The standard approach to solve HJB equation by finite difference schemes is to discretize the derivatives in HJB equation and to solve the resulting finite dimensional control problem. The nonlinear discretized equations are often solved using policy iteration schemes (see e.g., [1,2,[10][11][12][13][14][15]20,21,26,27,30,31]). Among them, the work by [20,21] and [2] outlines the theory and implementation of the schemes for solving the coupled HJB equations arising in the American options under regime switching models. No detailed convergence proofs are given therein.
In this paper we propose a different way to solve the problems of terminal wealth utility maximization under regime switching models. The HJB system for the problems is composed by d coupled HJB equations, where d is the number of regime states. Fu et al. [16] introduce a functional operator to generate a sequence of value functions and show that the optimal value function is the limit of this sequence. To get the value functions in the sequence, it needs to solve d decoupled HJB equations in each iterative step. Thus the coupled HJB equations are separated by d single HJB equations and henceforth it seems to be much simpler to solve although the iterations are involved. We study the iterative FDMs with policy iterations for solving the sequence of decoupled HJB equations in [16] and prove the convergence. We use several examples to show that solving the decoupled HJB is more efficient than solving the coupled ones.
The regime switching model allows parameters of asset price dynamics to depend on a finite state Markov chain process. It provides good flexibility for characterizing macro market uncertainties while preserves analytic tractability for underlying asset price dynamics. Hamilton [17] introduces a regime switching model for nonstationary time series and business cycles. Hardy [18] applies a two-regime model to provide a good fit to monthly stock market returns. There has been active research in portfolio optimization with regime switching models. Zhang et al. [35] and Yin et al. [33] study the trading rules in a regime switching market. Zhou and Yin [36] investigate the mean-variance portfolio optimization in regime switching model. Canakoglu and Özekici [9] discuss the HARA utility maximization in a regime switching model. Honda [19], Sass and Haussmann [29], and Rieder and Bäuerle [28] solve portfolio optimization problems with partial information and regime switching drift processes. Bäuerle and Rieder [5] and Fu et al. [16] show that the value function satisfies the HJB system of fully coupled nonlinear PDEs and prove the verification theorem. For a power or logarithmic utility function, the HJB equations can be reduced to a system of linear ODEs which are then solved with matrix exponentials. For general utility functions, it seems not possible to solve the system of HJB equations analytically. Ma et al. [23] develop the dual control monte-carlo methods to compute the tight bounds of value function in regime switching utility maximization, but it is not possible to guarantee the convergence in theory and the computation of the lower bound is rather time-consuming. The iterative FDMs with iteration policy developed in this paper is proved to be convergent and numerical comparison is made to show that it is much faster than computing the lower bound using the approach in [23].
The remaining parts of the paper are arranged as follows. In Sect. 2, we introduce the utility maximization under regime switching models and the HJB system and analyze the iterative FDMs with policy iterations for the sequence of decoupled HJB equations. In Sect. 3, we carry out a variety of numerical examples to test the convergence and compare efficiency of the proposed algorithms with the existing methods. Conclusions are given in the final section. The standard FDMs with policy iterations for the coupled HJB equations are given in the appendix.

Discretization of the Decoupled HJB Equations
Consider a fixed time horizon [0, T ]. Let ( , F , P) be a complete probability space, W a standard brownian motion, α a continuous time finite state observable Markov Chain process (MCP), which are independent of each other, and let {F t } be the natural filtration generated by W and α completed with all P-null sets.
We identify the state space of {α t } as a finite set of unit vectors E := {e 1 , e 2 , . . . , e d } where e i ∈ R d is a column vectors with one in the i-th position and zeros elsewhere, j = 1, . . . , d. Denote by Q = (q i j ) d×d the generator of the Markov Chain {α t } with q i j ≥ 0 for i = j and d j=1 q i j = 0 for each j ∈ D := {1, . . . , d}. The MCP α has a semi-martingale representation.
where Q is the transpose of Q, M is a purely discontinuous square integrable Martingale with initial value zero. Assume the financial market consists of one risk-free bond and one risky stock. The bond and stock price processes B and S are assumed to follow the stochastic differential equations (SDE) where r t = rα t , μ t = μα t , σ t = σ α t and r = (r 1 , . . . , r d ) is a vector of risk-free interest rates with r i being the rate in regime i, and μ = (μ 1 , . . . , μ d ) and σ = (σ 1 , . . . , σ d ) are vectors of return and volatility rates of the risky asset. Assume all rates are positive constants. Denote by θ := (θ 1 , . . . , θ d ) the vector of market prices of risk with θ i = μ i −r i σ i for i ∈ D. Let X be the wealth process of a portfolio comprising the bond B and the stock S. The wealth process X satisfies the SDE: where π t is a progressively measurable control process and represents the proportion of wealth X t invested in risky asset S t and θ t = θα t is the market price of risk at time t.
The utility maximization problem is defined by: where U is a utility function that is continuous, increasing and concave on [0, ∞]. Stochastic control is a standard method that to solve problem (1). To do so, we define the value functions where E t,x, j is the conditional expectation operator given X t = x, α t = e j for j ∈ D and t := {π s , s ∈ [t, T ]} is the set of all admissible control strategies over [t, T ]. It is proved by [16] that for a continuous, strictly increasing and concave utility function U , the optimal value functions V (t, x, j), for j ∈ D, satisfy the following system of HJB equations, on [0, T ] × (0, +∞) and the terminal and boundary conditions are given by where q j := d =1, = j q j , and the boundary condition (5) will be specified for concrete problems in the following sections. The verification theorem is also given by [16]. Moreover, [16] define a functional operator , and claim that the sequences V (m) (t, x, j) converge to the value function V (t, x, j) as m tends to ∞. The function V (0) (t, x, j) is computed by Using dynamic programming principle and variable transformation τ , Eq. (6) leads to the following HJB equations for j ∈ D on [0, T ] × (0, +∞) with terminal and boundary conditions where φ(t, . This treatment reduces a system of fully coupled HJB equations (2) to a sequence of decoupled HJB equations (8).
We mainly study the iterative FDMs with policy iterations for the system of decoupled HJB equations (8). A grid is constructed consisting of a set of M + 1 nodes {x 0 , . . . , be the approximation to V (m+1) (τ n , x i , j). Equation (8) can be discretized by a standard finite difference method to give with V (m),n 0 For the convenience of analysis, (12) is re-written as where π ( j) n+1 ∈ arg sup where at each node, ξ ∈ {0, 1} is chosen to ensure that α n+1 ) are positive. For ease of analysis, we can also write the Eq. (13) coupled with boundary conditions (10) and (11) into the matrix form. Let Then (13) can be written as where From [3,4,13], it is known that the stability, consistency and monotonicity of the discretization can ensure the convergence to the viscosity solution. So we will analyze the stability, consistency and monotonicity of (12) or equivalent form (13) or its matrix form (15). (11) is bounded, then the fully implicit iterative FDMs (12) are stable,

Lemma 2.1 (Stability of the iterative FDMs) If boundary function φ in
and From (16), we derive that , there must exist i * and j * , such that (17) gives that So we have Combining (18) and (19), we obtain Iteratively using (20) gives that To proceed the analysis, it is convenient to denote (12) as where G j is defined by the left-hand side minus the right-hand side of (12), and denote (8) as where F j is defined by the left-hand side minus the right-hand side of (8) and , j ∈ D is the current regime state. Next we give the definitions of the upper and lower semi-continuous envelopes of function F j .

Definition 2.1
The upper and lower semi-continuous envelopes of function F j are defined respectively by where B(·, •) denotes the neighborhood with center · and size •.
We now give the definitions of the iteration-based viscosity sub-solution, the iterationbased viscosity super-solution and the iteration-based viscosity solution of (22) as follows. (7), then V (m+1) is called the iteration-based viscosity sub-solution of (22). (7), then V (m+1) is called the iteration-based viscosity super-solution of (22). (iii) If it is both a sub-solution and super-solution of (22), then we call that V (m+1) is the iteration-based viscosity solution of (22). where We expand ϕ (m+1) at node ( τ , x, j) with the Taylor series for (23) to give that Therefore, combining (23) with (24) gives that Consequently, it follows from (25) that Thus the proof is complete. (11) is bounded, then the implicit iterative FDMs (12) are monotone in the sense that for any ρ 1 , ρ 2 , ρ 3 , ρ 4 ≥ 0, it holds true that
From Lemmas 2.1, 2.2, 2.3, we know that (12) or (15) is a consistent, stable, monotone discretization. In [3,4,13,20,21,26,27], they all mention that a consistent, stable, monotone discretization converges to the viscosity solution. In this paper to prove the convergence of the iterative FDMs (12) for solving (8), we must specially deal with the operator iterations from m-th step to (m + 1)-th step. The result is presented in the following Theorem 2.1.  (7) is equal to the exact value. Using Lemmas 2.1, 2.2, 2.3, 2.4 and following the lines in [4], we can prove that the solution V (1),h,ρ of (12) converges to the unique viscosity solution V (1) of the nonlinear PDE (8) as ρ → 0 and h → 0. Without loss of generality, we assume that the solution V (m),h,ρ of (12) converges to the unique viscosity solution V (m) of the nonlinear PDE (8) as ρ → 0 and h → 0. To complete the proof of this theorem by methods of induction, we only need to prove that the theorem holds true for m + 1. Let

Theorem 2.1 (Convergence of the iterative FDMs) Assumed that the original HJB equation
where h = x and ρ = τ are the spatial and temporal mesh sizes.
To implement the iterative FDM scheme (12), we need the following algorithm of iteration policy.

4:
for k = 0, 1, 2, . . . do 5: solve 6:  Proof We will first prove that this algorithm is convergent by showing that the sequences (V( j)) k for k ≥ 1 are non-decreasing and bounded. Subtracting the equations for steps k and k + 1 on line 6 in Algorithm 1 leads to that

From Algorithm 1 (line 7), we know that
Therefore, From (14) and α n+1 ) k ) has positive diagonals, non-positive off-diagonals, and is diagonally dominant. So it is an M-matrix. Therefore, we have i.e., the sequences (V( j)) k for k ≥ 1 are non-decreasing. Now we prove that the sequences are bounded. To this end, let Since V (m+1),n ( j), D (m),n+1 ( j) and φ n+1 ( j) − φ n ( j) are bounded with infinity norm, we know that b( j) is bounded. Using the notation of b( j), the equation on line 6 in Algorithm 1 can be written as ). Since all the coefficients, α n+1 i , β n+1 i , q j , are positive, we derive that k+1 . Then we derive that This gives that V max ≤ B max . ThereforeV k+1 is bounded from above. Consequently, the non-decreasing sequencesV k+1 in Algorithm 1 are convergent. Now we prove that the solution of Algorithm 1 is unique. To this end, letV 1 andV 2 be the two solutions of Algorithm 1, i.e., where Subtracting Eq. (30) from Eq. (29) gives (2)) is an M-matrix, we haveV 1 ≥V 2 . In the same manner, we can proveV 2 ≥V 1 . Therefore, we obtain thatV 1 =V 2 .

Numerical Examples
In this section, we solve several examples using the iterative FDMs with policy iterations for the sequence of decoupled HJB equations and the coupled ones for power, non-HARA and Yaari utility functions. The iterative FDMs with policy iterations for the coupled HJB equations stem from [20], which solve the American option pricing under regime switching. But the scheme has to be modified for the utility maximization, since the policy for the utility maximization is different from that for American options. For convenience to the readers, we provide the scheme in the appendix. Moreover the boundary conditions to the HJB system are constructed. The boundary conditions (5) are constructed by where the second equality is calculated by [8]. The construction of the boundary condition is motivated by that we allocate all of the wealth measured by the utility to the risk-free bond over [t, T ]. To verify the boundary condition (32) is correct, we compare it with the exact ones where the expression of a(t, j) is given by [16].
In Table 1, FDM-D-HJB denotes the iterative FDMs with policy iterations for solving the decoupled HJB and FDM-C-HJB for the coupled ones, N , M are the number of time and space mesh. The benchmark value is calculated by explicit formula given by (see [16]). The benchmark values are 2.19913 and 2.08313 respectively for the current regime state being 1 and 2. The numerics in Table 1 show that the approach is convergent and FDM-D-HJB takes much less time than FDM-C-HJB. Figures 1 and 2 test the convergence rates of FDM-C-HJB and FDM-D-HJB for space and time, respectively. The absolute value of the slope for the log-scale plots is just the convergence rate. Figures 1 and 2 respectively show that the convergence rate for space is about 2 and for time approximately 1. Example 3.2 For comparison, we use example in [23], which considers the non-HARA utility function . In Table 2, we observe that the value computed by the FDM-D-HJB and FDM-C-HJB is between the lower and upper bound, which shows that the computation is correct. Also we see that the FDM-D-HJB is more efficient than the FDM-C-HJB. Example 3. 3 We consider the Yaari utility function: U (x) = min(x, H ) with H = 2. Since the second derivative of the Yaari utility function is 0, which leads to a degenerate equation, we shall use the smoothing technique which is similar to [32]. Define the approximate utility function where a, b are positive constants. The riskless interest rates, return and volatility rates of risky asset are given by, The initial wealth at time t = 0 is x = 1, the investment period T = 1, the boundary condition given by (32) and the smoothing parameter ε = 10 −6 . Since the threshold H of the Yaari utility function is given by 2, it is reasonable to take x max = 2.
In Table 3, the number of time meshes is taken as N = 4000 and the number of space uniform meshes M = 200. When τ is close to 0, the value function is close to the utility function whose first-order derivative is discontinuous and therefore more meshes points are needed to improve the accuracy. Based on this observation, the nonuniform meshes are designed as follows: The number of space uniform meshes is taken as M = 200 for τ ∈ [0, T /2] and M = 100 for τ ∈ [T /2, T ]. The numerics in Table 3 show that for the Yaari utility, the values fall between the lower and upper bounds and still the FDM-D-HJB takes less time than FDM-C-HJB. Moreover, the average computational time with nonuniform meshes is about 52% of that with uniform meshes for FDM-D-HJB and 57% for FDM-C-HJB under almost the same accuracy. This shows that the nonuniform meshes can be used to improve the accuracy of the algorithm.

Conclusions
In this paper, we extend the finite difference methods (FDMs) with policy iterations to the HJB system arising from regime switching utility maximization problems. The coupled HJB equations and the sequence of decoupled HJB equations derived by [16] are solved respectively by the standard and iterative FDMs with policy iterations. Numerical examples for power, non-HARA, and Yaari utilities are conducted to exhibit the accuracy and efficiency of the approach and show that solving the sequence of decoupled HJB equations is more efficient than the coupled one. The convergence of the approach is proved and some new techniques (e.g., introducing of the iteration-based viscosity solution) are used to overcome the difficulties caused by the errors from the iterative FDMs for solving of the sequences of HJB equations and the policy iterations. In the future it will be interesting to study the numerical methods for the high-dimensional HJB equations arising in the utility maximization based on multiple stochastic factors. To avoid the curse of the dimensionality, it worth to study the radial basis function methods and the kernel-based methods which are proposed by [22] and [24] for solving the linear partial differential equations arising in option pricing.

Compliance with Ethical Standards
Conflict of interest The authors declared that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
Let V n i ( j) be the approximation of V (τ n , x i , j). Equation (34) can be discretized by a standard FDM which can be re-written as where π ( j) n+1 ∈ arg sup At each node, in order to ensure α n+1 i (π Define matrix operator A(π n+1 ) by