Abstract
Based on the method of dynamic programming, this paper uses analysis methods governed by the nonlinear and inhomogeneous partial differential 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 finite difference 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 finite difference method. Lastly, experiments are conducted to show the performance of the proposed algorithm.
1. Introduction
How to optimally allocate assets and optimally consume are extremely important and difficult topics in portfolio management [1–3]. These 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–7]. This adds further complication to the problem. The 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–10], complete markets which means that the stock movement and the volatility movement are either perfectly correlated or perfectly anticorrelated [9–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 three-level Crank–Nicolson scheme and Richardson’s extrapolation technique. The Crank–Nicolson scheme has second-order 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.
This 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.
2. Financial Market and Stochastic Control
We consider a market consisting of one riskless asset , whose price is governed bywith a constant risk-free interest rate and a risky asset modeled as
In (2), and are the return and the stochastic volatility of the stock price , respectively. is the stochastic variance of . Empirical studies show presence of mean reversion in the stock movements [14]. Heston model [5] is selected for , namely,
Here, and are the increments of the Wiener processes under a probability . The correlation between and is , namely, . We assume is a constant. In (3), is the long-run average variance (i.e., as tends to infinity, the expected value of tends to ), is the rate at which reverts to , and is the volatility of the stock variance . The parameters are positive constants and need to satisfy the Feller condition, , to ensure that is strictly positive. The risk premia is defined as
Following [1, 5, 15–17], we assume is a constant. This means the stock excess return is proportional to the stock variance.
Consider an investor who has an initial wealth and needs to determine strategies for asset allocation and consumption over an investment horizon . Let be the investor’s wealth at time . The strategies consist of an asset allocation rate and a consumption rate , which mean he/she allocates to the risky asset and to the riskless asset at time and consumes over the time interval . Thus, under the strategies and , the wealth process is governed by
The goal is to maximize the expected utilities over the investment horizon, namely,
In (6), and are control variables for this optimization problem. is the expectation operator under the probability . is the subjective discount rate, namely, the time preference of the investor. The larger is, the more weight the investor puts on the present than on the future. The parameter determines the relative importance between intertemporal consumption and the terminal wealth. and 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. Therefore, we adopt the CRRA utility function for and :where , , and are positive constants. Since and stand for the intertemporal consumption utility and the terminal wealth utility of the same investor, we use the same in and . However, and can be different since and have different dimensions.
Let be the value function of problem (6), which is given bywhere is the filtration associated to the stochastic processes in this problem. The terminal condition is obtained by setting in (8):
Based on the HJB dynamic programming procedure, is governed bywith the optimal strategies and determined by
After substituting expressions (11) and (12) into (10), one obtains an equation for the value function :
Based on the terminal condition and the scaling property of (13), it is reasonable to guess thatwhere , (13) becomeswithand (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.
3. Numerical Method
In this section, we develop a numerical method for solving (15). For the sake of conciseness of our expressions, we rewrite (15) aswith the initial condition , where
3.1. Crank–Nicolson Scheme and Richardson’s Extrapolation
We use a three-level Crank–Nicolson scheme (see [18–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 . Since the boundary conditions at and at are not known, we use one-sided difference method at these two numerical boundaries. Step sizes and are used to discretize and , respectively. Thus, and . We adopt the standard notation .
The three-level Crank–Nicolson scheme involves the levels , , and . It is straightforward to discretize all linear terms in (19) with second-order errors, namely,
The nonlinear term has two factors and . We discretize the factor at level and approximate the factor as an average between at level and that at level , namely,
This discretization scheme leads to a set of linear equations. Based on the expressions given by (21) and (22), equation (19) can be discretized asfor and , where is the maximal value of andwith
Since (23) is not applicable to the boundaries at and , we used one-sided difference technique to discretize (19) and obtained the boundary equations in the following. It is straightforward to show that, at , we haveand at , we have
By substituting these expressions into (19), we havewherewith
From (23), (28), and (29), the numerical solution of (19) for is determined by the following system of linear equations:
After eliminating , , , and , (33) can be transformed into the following tridiagonal matrix form for :where
The method given by (21) and (22) is a two-step method, namely, depends on and . At the zeroth step, is given by initial condition (16), namely, , for . We now determine , the solution at the first step. Performing Taylor expansion on at giveswhere is given by (16) and , , and can be determined analytically from (19):
The details of derivations for , , and are given in Appendix A. From (36), is given bywith an error of .
Knowing , the numerical solutions of optimal portfolio and consumption rules can be obtained from (17) and (18):for , where is given by (34) and is given by
In summary, our numerical solutions for and are determined by (16) and (38), respectively, and the numerical solutions for with are determined by (34). The numerical solution of obtained by the three-level Crank–Nicolson scheme has an accuracy of .
3.2. Performing Richardson’s Extrapolation
To further improve the accuracy of the numerical method, we apply Richardson’s extrapolation technique to . We will choose proportional to . Let represent obtained by (34) with a step size . Then,where is the exact value. We perform two computations with the step sizes and , respectively. Then, we have the following two equations:
From (43) and (44), we solve and obtain an expression based on Richardson’s extrapolation technique:
After substituting into (39) and (40), we obtain the expressions for and with an accuracy of :where is given by (45) and 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 , , , , , and . Here, we choose .(i)Step 1: initialize by initial condition (16), namely, , for .(ii)Step 2: initialize by (38) for .(iii)Step 3: for , knowing and , for , can be determined from (34), which is in a tridiagonal form and can be easily and efficiently solved. has an accuracy of .(iv)Step 4: to obtain and , we substitute from Steps 1–3 into (39) and (40). This provides the numerical solutions for optimal strategies and without extrapolation, which have accuracy of .(v)Step 5: to obtain , , and , we repeat Steps 1–3 with the step size to obtain . Then, from (45), (46), and (47), we obtain , , and , all of which have accuracy of .
In the next section, we will verify the accuracy of the numerical solutions without extrapolation and those with extrapolation.
4. Validation Study of the Numerical Method
Equation (19) is an inhomogeneous equation. However, since the inhomogeneous term 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 . Let be the solution of the homogeneous equation of (19), namely, the case of . Then, satisfieswith the initial condition . Following the procedure outlined in [10], the exact solution for can be obtained. After expressing asfrom (49), and are governed byand the solutions arewhere and . From (17) and (18), we obtain the exact solutions of optimal strategies and for the case of :where is given by (52).
The exact solutions , , and for the case of 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 . Since neither the inhomogeneous term nor the constant initial condition affects the stability or the accuracy of a numerical method, the accuracy and the stability of the method remain valid for . The numerical results for are presented at the end of this section.
4.1. Numerical Validation
To set parameters for numerical validation, we use the estimation values of the parameters , and given in [4, 15, 22, 23] and the historical records of . These values are listed in Table 1.
We note that since determines the wealth scale and determines the temporal scale, without loss of generality, we choose in this study.
For the range of the state variables and , we consider . 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 in the interval (to eliminate possible influence from the numerical boundary, the authors choose in their numerical computations). Since wealth does not appear in (19), (17), and (18), its value is irrelevant in our study.
In Table 2, we show the comparison between , , and . is the exact solution of (49). When , given by (34) is the numerical solution of (49) without performing Richardson’s extrapolation and given by (45) is the numerical solution of (49) after performing Richardson’s extrapolation. The relative errors in and , namely, and , are shown in the last two columns of Table 2.
In Table 3, we show the comparison between the exact solution given by (54), the numerical solution determined from (39) without performing Richardson’s extrapolation, and determined from (46) after performing Richardson’s extrapolation. The relative errors in and , namely, and , are shown in the last two columns of Table 3.
There is no error in since both numerical and theoretical values of are zero.
In Table 4, we show the global relative errors and the computational times of and for different values of and . The global relative error in is defined as the maximum of the local relative errors between and in the domain and . The global relative error in is defined in the same way.
Table 4 confirms that our numerical solutions and have accuracy at orders of and that and have accuracy at orders of . Therefore, 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.
4.2. 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 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 and , 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 and . 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 and , the proposed algorithm is 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 and . Tables 6–9 show that the proposed algorithm before and after performing the extrapolation technique has second-order accuracy and third-order accuracy, respectively, for the extensive sets of parameters.
4.3. Evidence for Stability and Convergence
To provide the evidence for stability, we calculate the maximum relative errors of for , respectively. The 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 be the maximum relative error; then, the convergence order is given by
In Table 10, the maximum relative errors and convergence orders of for different values of and are given. It shows that the convergence orders are always equal to 2.0 as and go to zero, which numerically verifies the convergence of the proposed algorithm.
4.4. Numerical Results for
We have confirmed the accuracy of our numerical solutions for . This guarantees that our numerical solution for will also have an accuracy of . In Tables 11–13, we present the numerical solutions of , , and for , and 0.9 with step sizes . In Tables 14 and 15, we show the results for and 10 with and . All other parameter values in Tables 11–15 are the same as the ones in Table 1. All digits shown in Tables 11–15 are exact, in the sense that they do not change when we further refine the values of and . Therefore, we have provided the first set of exact numerical data for optimal asset allocation and consumption strategies.
5. 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. The 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. The 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.
Appendix
A. Derivations for , , and
From (19), we obtain
By taking the derivative of (A.1) with respect to , we obtain
By taking the derivative of (A.2) with respect to , we obtain
By setting in (A.1) and using the initial condition , we obtain
By setting in (A.2) and using the initial condition and (A.4), we obtain
By setting in (A.3) and using the initial condition , (A.4), and (A.5), we obtain
B. HARA Utility Setting
In this appendix, the portfolio selection problem under HARA utility setting is considered. We would like to show that the proposed method could also be applied to this problem under HARA utility setting.
We assume that HARA utility functions and are used for consumption utility and terminal wealth utility, respectively:
The market setting and investors’ objective are the same as before which are given by (2)–(6). Following the procedure in this paper, the value function satisfies the following equation:and optimal strategies are given by
Based on the terminal condition and the scaling property of this problem, a reasonable trial solution form for is assumed:where and need to be determined and . The initial conditions of and are
To determine and , we substitute expression (B.4) into (B.2) and obtain
By settingequation (B.6) becomes an equation involving only, namely,
Thus, we have achieved dimension reduction by removing dependence in (B.2). The solution of is obtained from (B.7) and (B.5):
Equation (B.8) has the same form as (19), which could be solved numerically by applying the proposed algorithm the same as that in Section 3. Then, the optimal strategies could be obtained:
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The work of L. Ge was supported by the Fundamental Research Funds for the Central Universities (220110004005040120 and 220110001002020043). The work of Q. Zhang was supported by the Research Grants Council of the Hong Kong Special Administrative Region, China (project CityU 11335816).