Numerical Solutions to Optimal Portfolio Selection and Consumption Strategies under Stochastic Volatility

Based on the method of dynamic programming, this paper uses analysis methods governed by the nonlinear and inhomogeneous partial diﬀerential equation to study modern portfolio management problems with stochastic volatility, incomplete markets, limited investment scope, and constant relative risk aversion (CRRA). In this paper, a three-level Crank–Nicolson ﬁnite diﬀerence scheme is used to determine numerical solutions under this general setting. One of the main contributions of this paper is to apply this three-level technology to solve the portfolio selection problem. In addition, we have used a technique to deal with the nonlinear term, which is another novelty in performing the Crank–Nicolson algorithm. The Crank–Nicolson algorithm has also been extended to third-order accuracy by performing Richardson’s extrapolation. The accuracy of the proposed algorithm is much higher than the traditional ﬁnite diﬀerence method. Lastly, experiments are conducted to show the performance of the proposed algorithm.


Introduction
How to optimally allocate assets and optimally consume are extremely important and difficult topics in portfolio management [1][2][3]. ese topics are important not only for theoretical consideration but also for applications in the financial industry. Early studies usually assumed the volatility of the risky asset to be a constant. However, in recent years, researchers found that volatility should be modeled as stochastic rather than deterministic [4][5][6][7]. is adds further complication to the problem. e optimal asset allocation and optimal consumption strategies are governed by the Hamilton-Jacobi-Bellman (HJB) equation. Due to the nonlinearity and inhomogeneity of this partial differential equation, no exact solution has been found. Furthermore, even numerical solutions are not available. In this paper, we present an accurate and efficient numerical method for solving this equation and generate the first set of accurate numerical solutions for this problem.
Due to the importance of portfolio selection under stochastic volatilities, several important theoretical works have been carried out, and exact solutions have been obtained under certain special settings, such as no consumption [8][9][10], complete markets which means that the stock movement and the volatility movement are either perfectly correlated or perfectly anticorrelated [9][10][11][12], or when investors have unit elasticity of intertemporal substitution of consumption [13].
In this paper, we consider this optimal stochastic control problem under a general setting: stochastic volatility, incomplete markets, finite investment horizons, and CRRA utility. Our numerical method combines a threelevel Crank-Nicolson scheme and Richardson's extrapolation technique. e Crank-Nicolson scheme has secondorder accuracy in terms of discretization error, and Richardson's extrapolation technique further improves the accuracy. We verify that our numerical method is accurate and efficient. is paper is organized as follows. In Section 2, we describe the model for financial market, the stochastic control optimization procedure, and the governing HJB equation for the optimal asset allocation and consumption strategies. In Section 3, we present our numerical method for solving the HJB equation. In Section 4, we verify the accuracy and the efficiency of our numerical method and present accurate numerical solutions for the optimal asset allocation strategy and the optimal consumption strategy. In the last section, we present our conclusions.

Financial Market and Stochastic Control
We consider a market consisting of one riskless asset B t , whose price is governed by with a constant risk-free interest rate r and a risky asset S t modeled as In (2), μ(v t ) and σ(v t ) are the return and the stochastic volatility of the stock price S t , respectively. v t � σ 2 t is the stochastic variance of S t . Empirical studies show presence of mean reversion in the stock movements [14]. Heston model [5] is selected for v t , namely, Here, dW v t and dW S t are the increments of the Wiener processes under a probability P. e correlation between dW v t and dW S t is ρ, namely, Corr((dv t /v t ), (dS t /S t )) � ρdt. We assume ρ is a constant. In (3), θ is the long-run average variance (i.e., as t tends to infinity, the expected value of v t tends to θ), κ is the rate at which v t reverts to θ, and ξ is the volatility of the stock variance v t . e parameters κ, θ, ξ are positive constants and need to satisfy the Feller condition, 2κθ > ξ 2 , to ensure that v t is strictly positive. e risk premia is defined as Following [1,5,[15][16][17], we assume A is a constant. is means the stock excess return is proportional to the stock variance.
Consider an investor who has an initial wealth w 0 and needs to determine strategies for asset allocation and consumption over an investment horizon [0, T]. Let w t be the investor's wealth at time t. e strategies consist of an asset allocation rate φ t and a consumption rate c t , which mean he/she allocates φ t w t to the risky asset and (1 − φ t )w t to the riskless asset at time t and consumes c t dt over the time interval [t, t + dt]. us, under the strategies φ t and c t , the wealth process is governed by e goal is to maximize the expected utilities over the investment horizon, namely, In (6), φ t and c t are control variables for this optimization problem. E is the expectation operator under the probability P. β is the subjective discount rate, namely, the time preference of the investor. e larger β is, the more weight the investor puts on the present than on the future. e parameter α determines the relative importance between intertemporal consumption and the terminal wealth. u 1 (·) and u 2 (·) are the investor's utility functions which measure the investor's degree of satisfaction with the outcomes from intertemporal consumption and terminal wealth, respectively.
CRRA utility functions have been widely adopted for modeling investors' behavior. erefore, we adopt the CRRA utility function for u 1 (·) and u 2 (·): with the optimal strategies φ * and c * determined by After substituting expressions (11) and (12) into (10), one obtains an equation for the value function V: Based on the terminal condition and the scaling property of (13), it is reasonable to guess that where τ � T − t, (13) becomes with and (11) and (12) become Equation (15) is a nonlinear and inhomogeneous partial differential equation. Since no closed-form solution is available for this equation, numerical computation plays a critical role for studying this important practical problem in modern finance. However, there are even no numerical solutions available in the literature.

Numerical Method
In this section, we develop a numerical method for solving (15). For the sake of conciseness of our expressions, we rewrite (15) as with the initial condition

Crank-Nicolson Scheme and Richardson's Extrapolation.
We use a three-level Crank-Nicolson scheme (see [18][19][20][21]) of second-order accuracy to solve the nonlinear and inhomogeneous partial differential equation given by (19) and use Richardson's extrapolation technique for further improving accuracy. Numerically, one can only solve (19) over a finite domain v ∈ [0, v max ]. Since the boundary conditions at v � 0 and at v � v max are not known, we use one-sided difference method at these two numerical boundaries.
Step sizes Δτ and Δv are used to discretize τ and v, respectively. us, τ � nΔτ and v � mΔv. We adopt the standard notation f n m � f(nΔτ, mΔv). e three-level Crank-Nicolson scheme involves the levels n − 1, n, and n + 1. It is straightforward to discretize all linear terms in (19) with second-order errors, namely, Complexity e nonlinear term f 2 v /f has two factors f v /f and f v . We discretize the factor f v /f at level n and approximate the factor f v as an average between f v at level n − 1 and that at level n + 1, namely, is discretization scheme leads to a set of linear equations. Based on the expressions given by (21) and (22), equation (19) can be discretized as for 0 < m < M and n ≥ 2, where M is the maximal value of m and with Since (23) is not applicable to the boundaries at m � 0 and m � M, we used one-sided difference technique to discretize (19) and obtained the boundary equations in the following. It is straightforward to show that, at m � 0, we have and at m � M, we have By substituting these expressions into (19), we have with From (23), (28), and (29), the numerical solution of (19) for n ≥ 2 is determined by the following system of linear equations: Complexity After eliminating d 3 (0), d 4 (0), d 3 (M), and d 4 (M), (33) can be transformed into the following tridiagonal matrix form for n ≥ 2: where e method given by (21) and (22) is a two-step method, namely, f n+1 depends on f n and f n− 1 . At the zeroth step, f 0 is given by initial condition (16), We now determine f 1 , the solution at the first step. Performing Taylor expansion on f(Δτ, v) at τ � 0 gives where f(0, v) is given by (16) and f τ (0, v), f ττ (0, v), and f τττ (0, v) can be determined analytically from (19): 3 5 v 3 + 3a 2 a 2 5 + 3a 6 a 2 5 v 2 + a 5 a 2 2 + 3a 5 a 6 a 2 + 2a 1 a 2 5 + 3a 3 a 2 5 + 2a 4 a 2 5 + 3a 5 a 2 6 v + a 2 a 3 a 5 + 3a 3 a 6 a 5 + a 3 6 Complexity e details of derivations for f τ (0, v), f ττ (0, v), and f τττ (0, v) are given in Appendix A. From (36), f 1 m is given by with an error of O(Δτ 4 ).
Knowing f, the numerical solutions of optimal portfolio and consumption rules can be obtained from (17) and (18): In summary, our numerical solutions for f 0 m and f 1 m are determined by (16)

Performing Richardson's Extrapolation.
To further improve the accuracy of the numerical method, we apply Richardson's extrapolation technique to f. We will choose Δv proportional to Δτ. Let f(τ n , v m , Δτ) represent f n m obtained by (34) with a step size Δτ. en, where f exact is the exact value. We perform two computations with the step sizes Δτ and Δτ/2, respectively. en, we have the following two equations: From (43) and (44), we solve f exact and obtain an expression based on Richardson's extrapolation technique: After substituting f extpl. into (39) and (40), we obtain the expressions for φ * and c * /w with an accuracy of O(Δτ 3 ): where f extpl. is given by (45) and (f extpl. ) v is given by For the purpose of giving a quick understanding of our method, Figure 1 presents a flowchart of the algorithm for solving (19) We also summarize the procedure in words in the following for obtaining the numerical solutions of f, φ * , c * , f extpl. , φ * extpl. , and c * extpl. . Here, we choose Δv � Δτ. (v) Step 5: to obtain f extpl. , φ * extpl. , and c * extpl. , we repeat Steps 1-3 with the step size Δτ/2 to obtain f(τ 2n , v 2m , (Δτ/2)). en, from (45), (46), and (47), we obtain f extpl. , φ * extpl. , and c * extpl. , all of which have accuracy of O(Δτ 3 ).
In the next section, we will verify the accuracy of the numerical solutions without extrapolation and those with extrapolation.

Validation Study of the Numerical Method
Equation (19) is an inhomogeneous equation. However, since the inhomogeneous term α (1/c) a 7 affects neither the stability nor the accuracy of the three-level Crank-Nicolson method, it is sufficient to conduct validation studies for the corresponding homogeneous equation, namely, for the case of α � 0. Let f be the solution of the homogeneous equation of (19), namely, the case of α � 0. en, f satisfies with the initial condition f(0, v) � 1. Following the procedure outlined in [10], the exact solution for f can be obtained. After expressing f(τ, v) as from (49), h 1 (τ) and h 2 (τ) are governed by  Complexity where a 8 � a 1 + a 4 and Δ � a 2 2 − 4a 5 a 8 . From (17) and (18), we obtain the exact solutions of optimal strategies φ * and c * for the case of α � 0: where h 1 (τ) is given by (52). e exact solutions f exact , φ * exact , and c * exact for the case of α � 0 given by (50), (54), and (55) offer a benchmark for testing the accuracy of our numerical solutions. We show that our numerical solutions are accurate and efficient for α � 0. Since neither the inhomogeneous term α (1/c) a 7 nor the constant initial condition (1 − α) (1/c) affects the stability or the accuracy of a numerical method, the accuracy and the stability of the method remain valid for α ≠ 0. e numerical results for α ≠ 0 are presented at the end of this section.

Numerical Validation.
To set parameters for numerical validation, we use the estimation values of the parameters κ, θ, ξ, ρ, A, β, and c given in [4,15,22,23] and the historical records of r. ese values are listed in Table 1.
We note that since a w determines the wealth scale and a w /a c determines the temporal scale, without loss of generality, we choose a c � a w � 1 in this study.
For the range of the state variables t and v, we consider T ≤ 100. Based on the historical records of the Chicago Board Options Exchange Volatility Index, a popular measure of the implied volatility of S&P 500 index options, we examine the numerical solutions for the instantaneous volatility σ t � �� v t √ in the interval [0.1, 0.8] (to eliminate possible influence from the numerical boundary, the authors choose v max � 2 in their numerical computations). Since wealth w does not appear in (19), (17), and (18), its value is irrelevant in our study.
In Table 2, we show the comparison between f exact , f num. , and f extpl. . f exact is the exact solution of (49). When α � 0, f num. given by (34) is the numerical solution of (49) without performing Richardson's extrapolation and f extpl. given by (45) is the numerical solution of (49) after performing Richardson's extrapolation. e relative errors in f num. and f extpl. , namely, | (f num. − f exact )/f exact | and | (f extpl. − f exact )/f exact |, are shown in the last two columns of Table 2.
In Table 3, we show the comparison between the exact solution φ * exact given by (54), the numerical solution φ * num.
ere is no error in c * since both numerical and theoretical values of c * are zero.
In Table 4, we show the global relative errors and the computational times of f num. and f extpl. for different values of Δτ and Δv. e global relative error in f num. is defined as the maximum of the local relative errors between f exact and f num. in the domain 0 ≤ τ ≤ 100 and 0 ≤ � v √ � σ ≤ 0.8. e global relative error in f extpl. is defined in the same way. Table 4 confirms that our numerical solutions f num. and φ * num. have accuracy at orders of Δτ 2 and that f extpl. and φ * extpl. have accuracy at orders of Δτ 3 . erefore, the extrapolation technique does improve the accuracy. For the same order of accuracy, the application of Richardson's extrapolation technique significantly saves the computational time.

Extensive Sets of Parameters.
In this section, some extensive sets of parameters are used to further validate the proposed method. Since this system contains several Complexity parameters, we vary only one of the parameters at a time and keep other parameters at their benchmark values as shown in Table 1. Extensive sets are given in Table 5.
By choosing the step size Δτ � 0.01 and Δv � 0.01, the proposed algorithm is conducted, and the relative errors are recorded both before and after performing the extrapolation technique. In Tables 6 and 7, we show the maximum relative errors of the proposed algorithm before and after performing the extrapolation technique, respectively, in the domain 0 ≤ τ ≤ 100 and 0 ≤ � v √ � σ ≤ 0.8. It can be found that the extrapolation technique does improve one-order accuracy for the ten sets of parameters in Table 5.
By taking half of the step size above, namely, choosing Δτ � 0.005 and Δv � 0.005, the proposed algorithm is  κ, θ, ξ, ρ, A, r, β, c, a c , and a Table 3: Comparison between the exact solution φ * exact , the numerical solution without Richardson's extrapolation φ * num. , and the numerical solution after Richardson's extrapolation φ * extpl. . e relative errors in numerical solutions are shown in the last two columns. Here, the numerical solutions are obtained with Δτ � Δv � 0.01.       9.6 × 10 − 7 6.7 × 10 − 7 Set 10 9.6 × 10 − 7 6.7 × 10 − 7 conducted, and the relative errors are recorded both before and after performing the extrapolation technique. In Tables 8 and 9, we show the maximum relative errors of the proposed algorithm before and after performing the extrapolation technique, respectively, in the domain 0 ≤ τ ≤ 100 and 0 ≤ � v √ � σ ≤ 0.8. Tables 6-9 show that the proposed algorithm before and after performing the extrapolation technique has second-order accuracy and thirdorder accuracy, respectively, for the extensive sets of parameters.

Evidence for Stability and Convergence.
To provide the evidence for stability, we calculate the maximum relative errors of f num. for Δτ � Δv � 0.04, 0.02, 0.01, 0.005, 0.0025, respectively. e results are shown in Figure 2. One can see that, as the step size goes to zero, the maximum relative errors tend to zero, which numerically verifies the stability of the proposed algorithm.
To provide the evidence for convergence, let e(Δτ, Δv) be the maximum relative error; then, the convergence order is given by convergence order � log 2 e(Δτ, Δv) e((Δτ/2), (Δv/2)) . (56) In Table 10, the maximum relative errors and convergence orders of f num. for different values of Δτ and Δv are given. It shows that the convergence orders are always equal to 2.0 as Δτ and Δv go to zero, which numerically verifies the convergence of the proposed algorithm. 7.1×10 − 9 5.8×10 − 9 Set 7 7.4×10 − 9 6.0×10 − 9 Set 8 3.2×10 − 9 2.5×10 − 9 Set 9 6.9×10 − 9 5.6×10 − 9 Set 10 6.9×10    Tables 11-15 are exact, in the sense that they do not change when we further refine the values of Δτ and Δv. erefore, we have provided the first set of exact numerical   Table 13: Numerical solutions for f, φ * , and c * /w after performing Richardson's extrapolation. Here, α � 0.9 and Δτ � Δv � 0.001. Other parameter values are the same as those in Table 1. data for optimal asset allocation and consumption strategies.

Conclusion
In this paper, we study the portfolio selection problem in a general setting under CRRA utility functions: stochastic volatility, incomplete markets, finite investment horizons, and consumption choice. To the best of our knowledge, no explicit solution or numerical result is available in the literature for this setting. We present an accurate and efficient numerical method for optimal asset allocation and optimal consumption strategies. e optimal strategies are depended on a solution to a nonlinear and inhomogeneous partial differential equation which is derived from the portfolio selection problem. A three-level Crank-Nicolson finite difference scheme, which has second-order accuracy, is used to determine numerical solutions. In addition, we have used a technique to deal with the nonlinear term, which is one of our main contributions. We believe that the technique to deal with the nonlinear term could be applied to other similar numerical problems. e Crank-Nicolson algorithm also has been extended to third-order accuracy by performing Richardson's extrapolation. Some experiments are conducted to verify the performance of the proposed algorithm. Based on this algorithm, we present the first set of accurate numerical solutions of optimal strategies. Since the portfolio selection problem under stochastic volatility is an important issue in modern finance, the proposed algorithm will be useful for further theoretical research and for applications in the financial industry.  Table 15: Numerical solutions for f, φ * , and c * /w after performing Richardson's extrapolation. Here, α � 0.5, c � 10, and Δτ � Δv � 0.001.
Other parameter values are the same as those in Table 1.