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

Authors: Lei Ge; Qiang Zhang
Journal: Complexity
 (2020)

---

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


---

## Body

## 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 Section2, 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 assetBt, whose price is governed by(1)dBt=rBtdt,with a constant risk-free interest rate r and a risky asset St modeled as(2)dSt=Stμvt+σvtdWtS.In (2), μvt and σvt are the return and the stochastic volatility of the stock price St, respectively. vt=σt2 is the stochastic variance of St. Empirical studies show presence of mean reversion in the stock movements [14]. Heston model [5] is selected for vt, namely,(3)dvt=κθ−vtdt+ξvtdWtv.Here,dWtv and dWtS are the increments of the Wiener processes under a probability P. The correlation between dWtv and dWtS is ρ, namely, Corrdvt/vt,dSt/St=ρ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 vt tends to θ), κ is the rate at which vt reverts to θ, and ξ is the volatility of the stock variance vt. The parameters κ,θ,ξ are positive constants and need to satisfy the Feller condition, 2κθ>ξ2, to ensure that vt is strictly positive. The risk premia is defined as(4)A=μt−rσt2≡μt−rvt.Following [1, 5, 15–17], we assume A is a constant. This means the stock excess return is proportional to the stock variance.Consider an investor who has an initial wealthw0 and needs to determine strategies for asset allocation and consumption over an investment horizon 0,T. Let wt be the investor’s wealth at time t. The strategies consist of an asset allocation rate φt and a consumption rate ct, which mean he/she allocates φtwt to the risky asset and 1−φtwt to the riskless asset at time t and consumes ctdt over the time interval t,t+dt. Thus, under the strategies φt and ct, the wealth process is governed by(5)dwt=φtwtStdSt+1−φtwtrdt−ctdt.The goal is to maximize the expected utilities over the investment horizon, namely,(6)supφt,ctE∫0Tαe−βtu1ctdt+1−αe−βTu2wT.In (6), φt and ct 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. 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. u1⋅ and u2⋅ 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 foru1⋅ and u2⋅:(7)u1ct=acct1−γ1−γ,forγ≠1,aclogct,forγ=1,u2wT=awwT1−γ1−γ,forγ≠1,awlogwT,forγ=1,where γ, ac, and aw are positive constants. Since u1⋅ and u2⋅ stand for the intertemporal consumption utility and the terminal wealth utility of the same investor, we use the same γ in u1⋅ and u2⋅. However, ac and aw can be different since c and w have different dimensions.LetVt,w,v be the value function of problem (6), which is given by(8)Vt,wt,vt=supφt,ctE∫tTαe−βt′u1ct′dt′+1−αe−βTu2wTℱt,where ℱT is the filtration associated to the stochastic processes in this problem. The terminal condition is obtained by setting t=T in (8):(9)VT,wT,vT=1−αe−βTu2wT.Based on the HJB dynamic programming procedure,V is governed by(10)0=supφ,cαe−βtu1c+Vt+rw+φwAv−cVw+κθ−vVv+12φ2w2vVww+φwρξvVwv+12ξ2vVvv,with the optimal strategies φ∗ and c∗ determined by(11)φ∗=−AVw+ρξVwvwVww,(12)c∗=Vwacαe−βt−1/γ.After substituting expressions (11) and (12) into (10), one obtains an equation for the value function V:(13)γ1−γacαe−βt1/γVw1−γ/−γ+rwVw+Vt+κθ−vVv+12ξ2vVvv−12vAVw+ρξVwv2Vww=0.Based on the terminal condition and the scaling property of (13), it is reasonable to guess that(14)Vτ,w,v=e−βT−τaww1−γ1−γfτ,vγ,where τ=T−t, (13) becomes(15)0=−fτ+12ξ2vfvv+κθ−v+1−γγAρξvfv−121−γ1−ρ2ξ2vfv2f+1−γA22γ2v+1−γrγ−βγf+αacaw1/γ,with(16)f0,v=1−α1/γ,and (11) and (12) become(17)φ∗=Aγ+ρξfvf,(18)c∗w=αacaw1/γf−1.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) as(19)−fτ+a1vfvv+a2v+a3fv+a4vfv2f+a5v+a6f+α1/γa7=0,with the initial condition f0,v=1−α1/γ, where(20)a1=12ξ2,a2=1−γγAρξ−κ,a3=κθ,a4=−121−γξ21−ρ2,a5=1−γ2γ2A2,a6=1−γrγ−βγ,a7=acaw1/γ.


### 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 v∈0,vmax. Since the boundary conditions at v=0 and at v=vmax 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. Thus, τ=nΔτ and v=mΔv. We adopt the standard notation fmn=fnΔτ,mΔv.The three-level Crank–Nicolson scheme involves the levelsn−1, n, and n+1. It is straightforward to discretize all linear terms in (19) with second-order errors, namely,(21)fτmn=fmn+1−fmn−12Δτ+OΔτ2,fvmn=12fm+1n+1−fm−1n+12Δv+fm+1n−1−fm−1n−12Δv+OΔτ2+OΔv2,fvvmn=12fm+1n+1−2fmn+1+fm−1n+1Δv2+fm+1n−1−2fmn−1+fm−1n−1Δv2+OΔτ2+OΔv2.The nonlinear termfv2/f has two factors fv/f and fv. We discretize the factor fv/f at level n and approximate the factor fv as an average between fv at level n−1 and that at level n+1, namely,(22)fv2fmn=fvfmn12fvmn+1+fvmn−1+OΔτ2=fm+1n−fm−1n4Δvfmnfm+1n+1−fm−1n+12Δv+fm+1n−1−fm−1n−12Δv+OΔτ2+OΔv2.This discretization scheme leads to a set of linear equations. Based on the expressions given by (21) and (22), equation (19) can be discretized as(23)e1mfm−1n+1+e2mfmn+1+e3mfm+1n+1=e4m+OΔτ2+OΔv2,for 0<m<M and n≥2, where M is the maximal value of m and(24)e1m=a1mΔv−λm2Δv,e2m=−2a1mΔv−1Δτ+a5mΔv+a6,e3m=a1mΔv+λm2Δv,e4m=−1Δτfmn−1−2α1/γa7−a1mΔvfm+1n−1−2fmn−1+fm−1n−1−λm2Δvfm+1n−1−fm−1n−1−a5mΔv+a6fmn−1,with(25)λm=a2mΔv+a3+a4mfm+1n−fm−1n2fmn.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(26)fv0,τ=14Δv−3f0n+1+4f1n+1−f2n+1−3f0n−1+4f1n−1−f2n−1+OΔτ2+OΔv2,fvv0,τ=12Δv22f0n+1−5f1n+1+4f2n+1−f3n+1+2f0n−1−5f1n−1+4f2n−1−f3n−1+OΔτ2+OΔv2,fv2f0,τ=−3f0n+4f1n−f2n8Δv2f0n−3f0n+1+4f1n+1−f2n+1−3f0n−1+4f1n−1−f2n−1+OΔτ2+OΔv2,and at m=M, we have(27)fvvmax,τ=14Δv3fMn+1−4fM−1n+1+fM−2n+1+3fMn−1−4fM−1n−1+fM−2n−1+OΔτ2+OΔv2,fvvvmax,τ=12Δv22fMn+1−5fM−1n+1+4fM−2n+1−fM−3n+1+2fMn−1−5fM−1n−1+4fM−2n−1−fM−3n−1+OΔτ2+OΔv2,fv2fvmax,τ=3fMn−4fM−1n+fM−2n8Δv2fMn3fMn+1−4fM−1n+1+fM−2n+1+3fMn−1−4fM−1n−1+fM−2n−1+OΔτ2+OΔv2.By substituting these expressions into (19), we have(28)d10f0n+1+d20f1n+1+d30f2n+1+d40f3n+1=d50+OΔτ2+OΔv2,(29)d4MfM−3n+1+d3MfM−2n+1+d2MfM−1n+1+d1MfMn+1=d5M+OΔτ2+OΔv2,where(30)d10=−3λ02Δv−1Δτ+a6,d20=2λ0Δv,d30=−λ02Δv,d40=0,d50=−1Δτf0n−1−2α1γa7−λ02Δv−3f0n−1+4f1n−1−f2n−1−a6f0n−1,(31)d1M=2a1MΔv+3λM2Δv−1Δτ+a5MΔv+a6,d2M=−5a1MΔv−2λMΔv,d3M=4a1MΔv+λM2Δv,d4M=−a1MΔv,d5M=−1ΔτfMn−1−2α1/γa7−a1MΔv2fMn−1−5fM−1n−1+4fM−2n−1−fM−3n−1−λM2Δv3fMn−1−4fM−1n−1+fM−2n−1−a5MΔv+a6fMn−1,with(32)λ0=a3,λM=a2MΔv+a3+a4M3fMn−4fM−1n+fM−2n2fMn.From (23), (28), and (29), the numerical solution of (19) for n≥2 is determined by the following system of linear equations:(33)d10d20d30d40e11e21e31e12e22e32⋱⋱⋱e1M−1e2M−1e3M−1d4Md3Md2Md1Mf0n+1f1n+1f2n+1⋮fM−1n+1fMn+1=d50e41e42⋮e4M−1d5M.After eliminatingd30, d40, d3M, and d4M, (33) can be transformed into the following tridiagonal matrix form for n≥2:(34)e20e30e11e21e31e12e22e32⋱⋱⋱e1M−1e2M−1e3M−1e1Me2Mf0n+1f1n+1f2n+1⋮fM−1n+1fMn+1=e40e41e42⋮e4M−1e4M,where(35)e20=d10−e32d30−e22d40e32e31e11,e30=d20−d40e32e12−e32d30−e22d40e32e31e21,e40=d50−d40e32e42−e32d30−e22d40e32e31e41,e1M=d2M−d4Me1M−2e3M−2−e1M−2d3M−e2M−2d4Me1M−2e1M−1e2M−1,e2M=d1M−e1M−2d3M−e2M−2d4Me1M−2e1M−1e3M−1,e4M=d5M−d4Me1M−2e4M−2−e1M−2d3M−e2M−2d4Me1M−2e1M−1e4M−1.The method given by (21) and (22) is a two-step method, namely, fn+1 depends on fn and fn−1. At the zeroth step, f0 is given by initial condition (16), namely, fm0=1−α1/γ, for 0≤m≤M. We now determine f1, the solution at the first step. Performing Taylor expansion on fΔτ,v at τ=0 gives(36)fΔτ,v=f0,v+fτ0,vΔτ+12fττ0,vΔτ2+16fτττ0,vΔτ3+OΔτ4,where f0,v is given by (16) and fτ0,v, fττ0,v, and fτττ0,v can be determined analytically from (19):(37)fτ0,v=a5v+a6f0,v+α1/γa7,fττ0,v=a52v2+a2a5+2a5a6v+a3a5+a62f0,v+a5v+a6α1/γa7,fτττ0,v=a53v3+3a2a52+3a6a52v2+a5a22+3a5a6a2+2a1a52+3a3a52+2a4a52+3a5a62v+a2a3a5+3a3a6a5+a63f0,v+a2v+a3a5α1/γa7+a5v+a62α1/γa7.The details of derivations forfτ0,v, fττ0,v, and fτττ0,v are given in Appendix A. From (36), fm1 is given by(38)fm1=fm0+fτ0,mΔvΔτ+12fττ0,mΔvΔτ2+16fτττ0,mΔvΔτ3,with an error of OΔτ4.Knowingf, the numerical solutions of optimal portfolio and consumption rules can be obtained from (17) and (18):(39)φmn∗=Aγ+ρξfvmnfmn,(40)cmn∗w=αacaw1/γfmn−1,for 0≤m≤M, where fmn is given by (34) and fvmn is given by(41)fvmn=−3f0n+4f1n−f2n2Δv,form=0,fm+1n−fm−1n2Δv,for0<m<M,3fMn−4fM−1n+fM−2n2Δv,form=M.In summary, our numerical solutions forfm0 and fm1 are determined by (16) and (38), respectively, and the numerical solutions for fmn with n≥2 are determined by (34). The numerical solution of f obtained by the three-level Crank–Nicolson scheme has an accuracy of OΔτ2+OΔv2.


### 3.2. Performing Richardson’s Extrapolation

To further improve the accuracy of the numerical method, we apply Richardson’s extrapolation technique tof. We will choose Δv proportional to Δτ. Let fτn,vm,Δτ represent fmn obtained by (34) with a step size Δτ. Then,(42)fτn,vm,Δτ=fexactτ,v+C1Δτ2+C2Δτ3+OΔτ4,where fexact is the exact value. We perform two computations with the step sizes Δτ and Δτ/2, respectively. Then, we have the following two equations:(43)fτn,vm,Δτ=fexactτ,v+C1Δτ2+C2Δτ3+OΔτ4,(44)fτ2n,v2m,Δτ2=fexactτ,v+C122Δτ2+C223Δτ3+OΔτ4.From (43) and (44), we solve fexact and obtain an expression based on Richardson’s extrapolation technique:(45)fextpl.nΔτ,mΔv=43fτ2n,v2m,Δτ2−13fτn,vm,Δτ=fexact+OΔτ3.After substitutingfextpl. into (39) and (40), we obtain the expressions for φ∗ and c∗/w with an accuracy of OΔτ3:(46)φextpl.∗=Aγ+ρξfextpl.vfextpl.,(47)cextpl.∗w=αacaw1/γfextpl.−1,where fextpl. is given by (45) and fextpl.v is given by(48)fextpl.vmn=16Δv−9fextpl.mn+16fextpl.m+1n−8fextpl.m+2n+fextpl.m+4n,form≤1,16Δvfextpl.m−2n−4fextpl.m−1n+4fextpl.m+1n−fextpl.m+2n,for1<m<M−1,16Δv9fextpl.mn−16fextpl.m−1n+8fextpl.m−2n−fextpl.m−4n,form≥M−1.For the purpose of giving a quick understanding of our method, Figure1 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∗, fextpl., φextpl.∗, and cextpl.∗. Here, we choose Δv=Δτ.(i)
Step 1: initializef by initial condition (16), namely, fm0=1−α1/γ, for 0≤m≤M.(ii)
Step 2: initializefm1 by (38) for 0≤m≤M.(iii)
Step 3: forn>2, knowing fmn−1 and fmn, for 0≤m≤M, fmn+1 can be determined from (34), which is in a tridiagonal form and can be easily and efficiently solved. fmn has an accuracy of OΔτ2.(iv)
Step 4: to obtainφ∗ and c∗, we substitute fmn from Steps 1–3 into (39) and (40). This provides the numerical solutions for optimal strategies φm∗n and cm∗n without extrapolation, which have accuracy of OΔτ2.(v)
Step 5: to obtainfextpl., φextpl.∗, and cextpl.∗, we repeat Steps 1–3 with the step size Δτ/2 to obtain fτ2n,v2m,Δτ/2. Then, from (45), (46), and (47), we obtain fextpl., φextpl.∗, and cextpl.∗, all of which have accuracy of OΔτ3.Figure 1
Flowchart of the algorithm.In the next section, we will verify the accuracy of the numerical solutions without extrapolation and those with extrapolation.


## 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 v∈0,vmax. Since the boundary conditions at v=0 and at v=vmax 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. Thus, τ=nΔτ and v=mΔv. We adopt the standard notation fmn=fnΔτ,mΔv.The three-level Crank–Nicolson scheme involves the levelsn−1, n, and n+1. It is straightforward to discretize all linear terms in (19) with second-order errors, namely,(21)fτmn=fmn+1−fmn−12Δτ+OΔτ2,fvmn=12fm+1n+1−fm−1n+12Δv+fm+1n−1−fm−1n−12Δv+OΔτ2+OΔv2,fvvmn=12fm+1n+1−2fmn+1+fm−1n+1Δv2+fm+1n−1−2fmn−1+fm−1n−1Δv2+OΔτ2+OΔv2.The nonlinear termfv2/f has two factors fv/f and fv. We discretize the factor fv/f at level n and approximate the factor fv as an average between fv at level n−1 and that at level n+1, namely,(22)fv2fmn=fvfmn12fvmn+1+fvmn−1+OΔτ2=fm+1n−fm−1n4Δvfmnfm+1n+1−fm−1n+12Δv+fm+1n−1−fm−1n−12Δv+OΔτ2+OΔv2.This discretization scheme leads to a set of linear equations. Based on the expressions given by (21) and (22), equation (19) can be discretized as(23)e1mfm−1n+1+e2mfmn+1+e3mfm+1n+1=e4m+OΔτ2+OΔv2,for 0<m<M and n≥2, where M is the maximal value of m and(24)e1m=a1mΔv−λm2Δv,e2m=−2a1mΔv−1Δτ+a5mΔv+a6,e3m=a1mΔv+λm2Δv,e4m=−1Δτfmn−1−2α1/γa7−a1mΔvfm+1n−1−2fmn−1+fm−1n−1−λm2Δvfm+1n−1−fm−1n−1−a5mΔv+a6fmn−1,with(25)λm=a2mΔv+a3+a4mfm+1n−fm−1n2fmn.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(26)fv0,τ=14Δv−3f0n+1+4f1n+1−f2n+1−3f0n−1+4f1n−1−f2n−1+OΔτ2+OΔv2,fvv0,τ=12Δv22f0n+1−5f1n+1+4f2n+1−f3n+1+2f0n−1−5f1n−1+4f2n−1−f3n−1+OΔτ2+OΔv2,fv2f0,τ=−3f0n+4f1n−f2n8Δv2f0n−3f0n+1+4f1n+1−f2n+1−3f0n−1+4f1n−1−f2n−1+OΔτ2+OΔv2,and at m=M, we have(27)fvvmax,τ=14Δv3fMn+1−4fM−1n+1+fM−2n+1+3fMn−1−4fM−1n−1+fM−2n−1+OΔτ2+OΔv2,fvvvmax,τ=12Δv22fMn+1−5fM−1n+1+4fM−2n+1−fM−3n+1+2fMn−1−5fM−1n−1+4fM−2n−1−fM−3n−1+OΔτ2+OΔv2,fv2fvmax,τ=3fMn−4fM−1n+fM−2n8Δv2fMn3fMn+1−4fM−1n+1+fM−2n+1+3fMn−1−4fM−1n−1+fM−2n−1+OΔτ2+OΔv2.By substituting these expressions into (19), we have(28)d10f0n+1+d20f1n+1+d30f2n+1+d40f3n+1=d50+OΔτ2+OΔv2,(29)d4MfM−3n+1+d3MfM−2n+1+d2MfM−1n+1+d1MfMn+1=d5M+OΔτ2+OΔv2,where(30)d10=−3λ02Δv−1Δτ+a6,d20=2λ0Δv,d30=−λ02Δv,d40=0,d50=−1Δτf0n−1−2α1γa7−λ02Δv−3f0n−1+4f1n−1−f2n−1−a6f0n−1,(31)d1M=2a1MΔv+3λM2Δv−1Δτ+a5MΔv+a6,d2M=−5a1MΔv−2λMΔv,d3M=4a1MΔv+λM2Δv,d4M=−a1MΔv,d5M=−1ΔτfMn−1−2α1/γa7−a1MΔv2fMn−1−5fM−1n−1+4fM−2n−1−fM−3n−1−λM2Δv3fMn−1−4fM−1n−1+fM−2n−1−a5MΔv+a6fMn−1,with(32)λ0=a3,λM=a2MΔv+a3+a4M3fMn−4fM−1n+fM−2n2fMn.From (23), (28), and (29), the numerical solution of (19) for n≥2 is determined by the following system of linear equations:(33)d10d20d30d40e11e21e31e12e22e32⋱⋱⋱e1M−1e2M−1e3M−1d4Md3Md2Md1Mf0n+1f1n+1f2n+1⋮fM−1n+1fMn+1=d50e41e42⋮e4M−1d5M.After eliminatingd30, d40, d3M, and d4M, (33) can be transformed into the following tridiagonal matrix form for n≥2:(34)e20e30e11e21e31e12e22e32⋱⋱⋱e1M−1e2M−1e3M−1e1Me2Mf0n+1f1n+1f2n+1⋮fM−1n+1fMn+1=e40e41e42⋮e4M−1e4M,where(35)e20=d10−e32d30−e22d40e32e31e11,e30=d20−d40e32e12−e32d30−e22d40e32e31e21,e40=d50−d40e32e42−e32d30−e22d40e32e31e41,e1M=d2M−d4Me1M−2e3M−2−e1M−2d3M−e2M−2d4Me1M−2e1M−1e2M−1,e2M=d1M−e1M−2d3M−e2M−2d4Me1M−2e1M−1e3M−1,e4M=d5M−d4Me1M−2e4M−2−e1M−2d3M−e2M−2d4Me1M−2e1M−1e4M−1.The method given by (21) and (22) is a two-step method, namely, fn+1 depends on fn and fn−1. At the zeroth step, f0 is given by initial condition (16), namely, fm0=1−α1/γ, for 0≤m≤M. We now determine f1, the solution at the first step. Performing Taylor expansion on fΔτ,v at τ=0 gives(36)fΔτ,v=f0,v+fτ0,vΔτ+12fττ0,vΔτ2+16fτττ0,vΔτ3+OΔτ4,where f0,v is given by (16) and fτ0,v, fττ0,v, and fτττ0,v can be determined analytically from (19):(37)fτ0,v=a5v+a6f0,v+α1/γa7,fττ0,v=a52v2+a2a5+2a5a6v+a3a5+a62f0,v+a5v+a6α1/γa7,fτττ0,v=a53v3+3a2a52+3a6a52v2+a5a22+3a5a6a2+2a1a52+3a3a52+2a4a52+3a5a62v+a2a3a5+3a3a6a5+a63f0,v+a2v+a3a5α1/γa7+a5v+a62α1/γa7.The details of derivations forfτ0,v, fττ0,v, and fτττ0,v are given in Appendix A. From (36), fm1 is given by(38)fm1=fm0+fτ0,mΔvΔτ+12fττ0,mΔvΔτ2+16fτττ0,mΔvΔτ3,with an error of OΔτ4.Knowingf, the numerical solutions of optimal portfolio and consumption rules can be obtained from (17) and (18):(39)φmn∗=Aγ+ρξfvmnfmn,(40)cmn∗w=αacaw1/γfmn−1,for 0≤m≤M, where fmn is given by (34) and fvmn is given by(41)fvmn=−3f0n+4f1n−f2n2Δv,form=0,fm+1n−fm−1n2Δv,for0<m<M,3fMn−4fM−1n+fM−2n2Δv,form=M.In summary, our numerical solutions forfm0 and fm1 are determined by (16) and (38), respectively, and the numerical solutions for fmn with n≥2 are determined by (34). The numerical solution of f obtained by the three-level Crank–Nicolson scheme has an accuracy of OΔτ2+OΔv2.


## 3.2. Performing Richardson’s Extrapolation

To further improve the accuracy of the numerical method, we apply Richardson’s extrapolation technique tof. We will choose Δv proportional to Δτ. Let fτn,vm,Δτ represent fmn obtained by (34) with a step size Δτ. Then,(42)fτn,vm,Δτ=fexactτ,v+C1Δτ2+C2Δτ3+OΔτ4,where fexact is the exact value. We perform two computations with the step sizes Δτ and Δτ/2, respectively. Then, we have the following two equations:(43)fτn,vm,Δτ=fexactτ,v+C1Δτ2+C2Δτ3+OΔτ4,(44)fτ2n,v2m,Δτ2=fexactτ,v+C122Δτ2+C223Δτ3+OΔτ4.From (43) and (44), we solve fexact and obtain an expression based on Richardson’s extrapolation technique:(45)fextpl.nΔτ,mΔv=43fτ2n,v2m,Δτ2−13fτn,vm,Δτ=fexact+OΔτ3.After substitutingfextpl. into (39) and (40), we obtain the expressions for φ∗ and c∗/w with an accuracy of OΔτ3:(46)φextpl.∗=Aγ+ρξfextpl.vfextpl.,(47)cextpl.∗w=αacaw1/γfextpl.−1,where fextpl. is given by (45) and fextpl.v is given by(48)fextpl.vmn=16Δv−9fextpl.mn+16fextpl.m+1n−8fextpl.m+2n+fextpl.m+4n,form≤1,16Δvfextpl.m−2n−4fextpl.m−1n+4fextpl.m+1n−fextpl.m+2n,for1<m<M−1,16Δv9fextpl.mn−16fextpl.m−1n+8fextpl.m−2n−fextpl.m−4n,form≥M−1.For the purpose of giving a quick understanding of our method, Figure1 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∗, fextpl., φextpl.∗, and cextpl.∗. Here, we choose Δv=Δτ.(i)
Step 1: initializef by initial condition (16), namely, fm0=1−α1/γ, for 0≤m≤M.(ii)
Step 2: initializefm1 by (38) for 0≤m≤M.(iii)
Step 3: forn>2, knowing fmn−1 and fmn, for 0≤m≤M, fmn+1 can be determined from (34), which is in a tridiagonal form and can be easily and efficiently solved. fmn has an accuracy of OΔτ2.(iv)
Step 4: to obtainφ∗ and c∗, we substitute fmn from Steps 1–3 into (39) and (40). This provides the numerical solutions for optimal strategies φm∗n and cm∗n without extrapolation, which have accuracy of OΔτ2.(v)
Step 5: to obtainfextpl., φextpl.∗, and cextpl.∗, we repeat Steps 1–3 with the step size Δτ/2 to obtain fτ2n,v2m,Δτ/2. Then, from (45), (46), and (47), we obtain fextpl., φextpl.∗, and cextpl.∗, all of which have accuracy of OΔτ3.Figure 1
Flowchart of the algorithm.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 α1/γa7 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. Then, f^ satisfies(49)−f^τ+a1vf^vv+a2v+a3f^v+a4vf^v2f^+a5v+a6f^=0,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(50)f^τ,v=eh1τv+h2τ,from (49), h1τ and h2τ are governed by(51)h1′τ−a1+a4h1τ2−a2h1τ−a5=0,withh10=0,h2′τ−a3h1τ−a6=0,withh20=0,and the solutions are(52)h1τ=2a5eΔτ−1a2+Δ−a2−ΔeΔτ,forΔ>0,−2a5τa2τ−2,forΔ=0,2a5−Δcot−Δτ/2−a2/−Δ,forΔ<0,(53)h2τ=−a3a8loga2+Δe−Δτ−a2+Δ2Δ+a6−a3Δ2a8−a2a32a8τ,forΔ>0,−a3a8log1−a2τ2+a6−a2a32a8τ,forΔ=0,−a3a8logcos−Δτ2−a2−Δsin−Δτ2+a6−a2a32a8τ,forΔ<0,where a8=a1+a4 and Δ=a22−4a5a8. From (17) and (18), we obtain the exact solutions of optimal strategies φ^∗ and c∗^ for the case of α=0:(54)φ^exact∗=Aγ+ρξh1τ,(55)c^exact∗w=0,where h1τ is given by (52).The exact solutionsf^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/γa7 nor the constant initial condition 1−α1/γ affects the stability or the accuracy of a numerical method, the accuracy and the stability of the method remain valid for α≠0. The numerical results for α≠0 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κ,θ,ξ,ρ,A,β, and γ given in [4, 15, 22, 23] and the historical records of r. These values are listed in Table 1.Table 1
Values for the parametersκ,θ,ξ,ρ,A,r,β,γ,ac, and aw.
ParameterκθξρArβγacawValue1.60480.04640.3796−76.70%1.551%0.06211We note that sinceaw determines the wealth scale and aw/ac determines the temporal scale, without loss of generality, we choose ac=aw=1 in this study.For the range of the state variablest 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=vt in the interval 0.1,0.8 (to eliminate possible influence from the numerical boundary, the authors choose vmax=2 in their numerical computations). Since wealth w does not appear in (19), (17), and (18), its value is irrelevant in our study.In Table2, 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. The 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.Table 2
Comparison between the exact solutionf^exact, the numerical solution without Richardson’s extrapolation f^num., and the numerical solution after Richardson’s extrapolation f^extpl.. The relative errors in numerical solutions are shown in the last two columns. Here, the numerical solutions are obtained with Δτ=Δv=0.01.
τσ=vf^exactf^num.f^extpl.f^num.−f^exact/f^exactf^extpl.−f^exact/f^exact0.10.10.99612021120.99612029040.99612021127.9×10−82.1×10−120.10.40.99193788440.99193768090.99193788442.1×10−76.2×10−120.10.80.97867207880.97867068450.97867207891.4×10−68.1×10−1110.10.95693653320.95693675140.95693653302.3×10−71.3×10−1010.40.93395165570.93395110520.93395165575.9×10−71.4×10−1110.80.86404514140.86404216230.86404514153.4×10−68.9×10−11100.10.60625318630.60625317020.60625318392.6×10−84.0×10−9100.40.58706961780.58706961360.58706961597.1×10−93.2×10−9100.80.52966763030.52966766620.52966762886.8×10−82.7×10−91000.10.00617630000.00617629860.00617629972.2×10−74.2×10−81000.40.00598086420.00598086300.00598086402.0×10−74.1×10−81000.80.00539607210.00539607140.00539607181.3×10−74.0×10−8In Table3, we show the comparison between the exact solution φ^exact∗ given by (54), the numerical solution φ^num.∗ determined from (39) without performing Richardson’s extrapolation, and φ^extpl.∗ determined from (46) after performing Richardson’s extrapolation. The relative errors in φ^num.∗ and φ^extpl.∗, namely, φ^num.∗−φ^exact∗/φ^exact∗ and φ^extpl.∗−φ^exact∗/φ^exact∗, are shown in the last two columns of Table 3.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.∗. The relative errors in numerical solutions are shown in the last two columns. Here, the numerical solutions are obtained with Δτ=Δv=0.01.
τσ=vφ^exact∗φ^num.∗φ^extpl.∗φ^num.∗−φ^exact∗/φ^exact∗φ^extpl.∗−φ^exact∗/φ^exact∗0.10.10.78316676090.78316727060.78316676096.5×10−71.5×10−110.10.40.78316676090.78316735680.78316676097.6×10−72.8×10−110.10.80.78316676090.78316764870.78316676081.1×10−69.1×10−1110.10.82219087610.82219246870.82219087411.9×10−62.4×10−910.40.82219087610.82219251740.82219087602.0×10−66.6×10−1110.80.82219087610.82219274480.82219087602.3×10−65.8×10−11100.10.83741215490.83741222330.83741214788.2×10−88.6×10−9100.40.83741215490.83741215890.83741215444.8×10−96.4×10−10100.80.83741215490.83741215650.83741215481.9×10−91.6×10−101000.10.83741219690.83741226520.83741218978.2×10−88.6×10−91000.40.83741219690.83741220090.83741219644.7×10−96.4×10−101000.80.83741219690.83741219850.83741219681.9×10−91.6×10−10There is no error inc∗ since both numerical and theoretical values of c∗ are zero.In Table4, we show the global relative errors and the computational times of f^num. and f^extpl. for different values of Δτ and Δv. The 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. The global relative error in f^extpl. is defined in the same way.Table 4
Maximum relative errors and computational times off^num. and f^extpl. in the domain 0≤τ≤100 and 0≤v=σ≤0.8.
Step sizef^num.f^extpl.ΔτΔvMax. rel. err.Comp. timeMax. rel. err.Comp. time0.020.021.5×10−50.393.3×10−70.840.010.013.8×10−61.524.2×10−83.310.0050.0059.6×10−76.025.2×10−913.180.00250.00252.4×10−724.026.3×10−1052.54Table4 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. 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 Table1. Extensive sets are given in Table 5.Table 5
Values for the parametersκ,θ,ξ,ρ,A,r,β,γ,ac, and aw.
ParameterκθξρArβγacawSet 12.50.04640.3796−76.70%1.551%0.06211Set 21.60480.080.3796−76.70%1.551%0.06211Set 31.60480.04640.6−76.70%1.551%0.06211Set 41.60480.04640.3796−50%1.551%0.06211Set 51.60480.04640.3796−76.70%1.251%0.06211Set 61.60480.04640.3796−76.70%1.555%0.06211Set 71.60480.04640.3796−76.70%1.551%0.15211Set 81.60480.04640.3796−76.70%1.551%0.06511Set 91.60480.04640.3796−76.70%1.551%0.06221Set 101.60480.04640.3796−76.70%1.551%0.06212By 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.Table 6
Maximum relative error forf^num. and φ^num. with Δτ=0.01andΔv=0.01.
SetMaximum relative errorf^num.φ^num.Set 15.8×10−63.9×10−6Set 23.6×10−62.7×10−6Set 33.6×10−63.9×10−6Set 44.0×10−61.8×10−6Set 52.4×10−62.0×10−6Set 64.0×10−62.8×10−6Set 74.2×10−62.9×10−6Set 82.1×10−61.4×10−6Set 93.8×10−62.7×10−6Set 103.8×10−62.7×10−6Table 7
Maximum relative error forf^extpl. and φ^extpl. with Δτ=0.01andΔv=0.01.
SetMaximum relative errorf^extpl.φ^extpl.Set 11.2×10−79.3×10−8Set 25.0×10−84.3×10−8Set 39.5×10−85.9×10−8Set 45.9×10−83.1×10−8Set 53.3×10−83.1×10−8Set 65.6×10−84.5×10−8Set 75.8×10−84.7×10−8Set 82.5×10−82.02×10−8Set 95.4×10−84.4×10−8Set 105.4×10−84.4×10−8By taking half of the step size above, namely, choosingΔτ=0.005 and Δv=0.005, 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 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 third-order accuracy, respectively, for the extensive sets of parameters.Table 8
Maximum relative error forf^num. and φ^num. with Δτ=0.005andΔv=0.005.
SetMaximum relative errorf^num.φ^num.Set 11.4×10−69.7×10−7Set 29.1×10−76.7×10−7Set 39.0×10−79.8×10−7Set 41.0×10−64.6×10−7Set 56.0×10−75.0×10−7Set 610.0×10−76.9×10−7Set 71.0×10−67.2×10−7Set 85.1×10−73.5×10−7Set 99.6×10−76.7×10−7Set 109.6×10−76.7×10−7Table 9
Maximum relative error forf^extpl. and φ^extpl. with Δτ=0.005andΔv=0.005.
SetMaximum relative errorf^extpl.φ^extpl.Set 11.6×10−81.2×10−8Set 26.4×10−97.5×10−9Set 32.1×10−87.6×10−9Set 47.6×10−94.0×10−9Set 54.2×10−94.0×10−9Set 67.1×10−95.8×10−9Set 77.4×10−96.0×10−9Set 83.2×10−92.5×10−9Set 96.9×10−95.6×10−9Set 106.9×10−95.6×10−9


### 4.3. Evidence for Stability and Convergence

To provide the evidence for stability, we calculate the maximum relative errors off^num. for Δτ=Δv=0.04,0.02,0.01,0.005,0.0025, 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.Figure 2
Maximum relative errors for decreasing step size.To provide the evidence for convergence, leteΔτ,Δv be the maximum relative error; then, the convergence order is given by(56)convergence order=log2eΔτ,ΔveΔτ/2,Δv/2.In Table10, 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.Table 10
Maximum relative errors and convergence orders off^num. in the domain 0≤τ≤100 and 0≤v=σ≤0.8.
Step sizeMaximum relative errorsConvergence orderΔτΔv0.020.021.53×10−52.00.010.013.83×10−62.00.0050.0059.58×10−72.00.00250.00252.40×10−72.0


### 4.4. Numerical Results forα≠0

We have confirmed the accuracy of our numerical solutions forα=0. This guarantees that our numerical solution fextpl. for α≠0 will also have an accuracy of OΔτ3. In Tables 11–13, we present the numerical solutions of fextpl., φextpl.∗, and cextpl.∗/w for α=0.1,0.5, and 0.9 with step sizes Δτ=Δv=0.001. In Tables 14 and 15, we show the results for γ=1 and 10 with α=0.5 and Δτ=Δv=0.001. 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 Δv. Therefore, we have provided the first set of exact numerical data for optimal asset allocation and consumption strategies.Table 11
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.1 and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.10.9765644460.7830377440.3238165870.10.40.9725289550.7830373940.3251602580.10.80.9597282710.7830362680.32949718710.11.2174435710.8175056160.25974737010.41.1910905830.8174484280.26549430510.81.1107573570.8172604090.284695630100.13.0734957700.8330103060.102888629100.42.9830468230.8329476710.106008315100.82.7117159650.8327395490.1166153721000.16.2276929850.8352401450.0507776741000.46.0374350620.8352066940.0523778331000.85.4674510730.8350951010.057838243Table 12
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.5 and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.10.7749377740.7828032100.9124691100.10.40.7718288440.7828022690.9161445400.10.80.7619663280.7827992380.92800266310.11.3689741000.8128728950.51652312610.41.3425532930.8127817070.52668805410.81.2617877680.8124844520.560400726100.16.0151811760.8323826420.117553696100.45.8400613900.8323127310.121078656100.85.3145274220.8320806890.1330516761000.113.9168104220.8352387810.0508095431000.413.4916571570.8352053110.0524106691000.812.2179612210.8350936530.057874368Table 13
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.9 and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.10.4096863870.7822441532.3156329550.10.40.4081604550.7822420302.3242900850.10.80.4033185340.7822352052.35219365810.11.2314519090.8082916300.77037786910.41.2105426960.8081887210.78368429410.81.1464004990.8078556960.827532175100.17.6867898130.8321316570.123417359100.47.4639749660.8320589530.127101618100.86.7952020830.8318177430.1396107561000.118.6674543490.8352383260.0508201751000.418.0971749900.8352048500.0524216241000.816.3887024200.8350931700.057886419Table 14
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.5, γ=1, and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.10.5468592821.5500000000.9143119940.10.40.5468592821.5500000000.9143119940.10.80.5468592821.5500000000.91431199410.10.9561778201.5500000000.52291528810.40.9561778201.5500000000.52291528810.80.9561778201.5500000000.522915288100.14.0343088511.5500000000.123936966100.44.0343088511.5500000000.123936966100.84.0343088511.5500000000.1239369661000.18.3139164411.5500000000.0601401281000.48.3139164411.5500000000.0601401281000.88.3139164411.5500000000.060140128Table 15
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.5, γ=10, and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.11.0247308490.1578338050.9105151780.10.41.0232359120.1578336810.9118454320.10.81.0184671920.1578332820.91611492210.11.8407898580.1695443640.50686556510.41.8270546360.1695304530.51067602110.81.7838731460.1694856200.523037747100.19.1822547040.1785759880.101612624100.49.0714269140.1785644610.102854049100.88.7260391460.1785270170.1069251441000.138.5460734750.1801360060.0242056561000.438.0501550490.1801328720.0245211351000.836.5060889640.1801226620.025558284


## 4.1. Numerical Validation

To set parameters for numerical validation, we use the estimation values of the parametersκ,θ,ξ,ρ,A,β, and γ given in [4, 15, 22, 23] and the historical records of r. These values are listed in Table 1.Table 1
Values for the parametersκ,θ,ξ,ρ,A,r,β,γ,ac, and aw.
ParameterκθξρArβγacawValue1.60480.04640.3796−76.70%1.551%0.06211We note that sinceaw determines the wealth scale and aw/ac determines the temporal scale, without loss of generality, we choose ac=aw=1 in this study.For the range of the state variablest 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=vt in the interval 0.1,0.8 (to eliminate possible influence from the numerical boundary, the authors choose vmax=2 in their numerical computations). Since wealth w does not appear in (19), (17), and (18), its value is irrelevant in our study.In Table2, 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. The 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.Table 2
Comparison between the exact solutionf^exact, the numerical solution without Richardson’s extrapolation f^num., and the numerical solution after Richardson’s extrapolation f^extpl.. The relative errors in numerical solutions are shown in the last two columns. Here, the numerical solutions are obtained with Δτ=Δv=0.01.
τσ=vf^exactf^num.f^extpl.f^num.−f^exact/f^exactf^extpl.−f^exact/f^exact0.10.10.99612021120.99612029040.99612021127.9×10−82.1×10−120.10.40.99193788440.99193768090.99193788442.1×10−76.2×10−120.10.80.97867207880.97867068450.97867207891.4×10−68.1×10−1110.10.95693653320.95693675140.95693653302.3×10−71.3×10−1010.40.93395165570.93395110520.93395165575.9×10−71.4×10−1110.80.86404514140.86404216230.86404514153.4×10−68.9×10−11100.10.60625318630.60625317020.60625318392.6×10−84.0×10−9100.40.58706961780.58706961360.58706961597.1×10−93.2×10−9100.80.52966763030.52966766620.52966762886.8×10−82.7×10−91000.10.00617630000.00617629860.00617629972.2×10−74.2×10−81000.40.00598086420.00598086300.00598086402.0×10−74.1×10−81000.80.00539607210.00539607140.00539607181.3×10−74.0×10−8In Table3, we show the comparison between the exact solution φ^exact∗ given by (54), the numerical solution φ^num.∗ determined from (39) without performing Richardson’s extrapolation, and φ^extpl.∗ determined from (46) after performing Richardson’s extrapolation. The relative errors in φ^num.∗ and φ^extpl.∗, namely, φ^num.∗−φ^exact∗/φ^exact∗ and φ^extpl.∗−φ^exact∗/φ^exact∗, are shown in the last two columns of Table 3.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.∗. The relative errors in numerical solutions are shown in the last two columns. Here, the numerical solutions are obtained with Δτ=Δv=0.01.
τσ=vφ^exact∗φ^num.∗φ^extpl.∗φ^num.∗−φ^exact∗/φ^exact∗φ^extpl.∗−φ^exact∗/φ^exact∗0.10.10.78316676090.78316727060.78316676096.5×10−71.5×10−110.10.40.78316676090.78316735680.78316676097.6×10−72.8×10−110.10.80.78316676090.78316764870.78316676081.1×10−69.1×10−1110.10.82219087610.82219246870.82219087411.9×10−62.4×10−910.40.82219087610.82219251740.82219087602.0×10−66.6×10−1110.80.82219087610.82219274480.82219087602.3×10−65.8×10−11100.10.83741215490.83741222330.83741214788.2×10−88.6×10−9100.40.83741215490.83741215890.83741215444.8×10−96.4×10−10100.80.83741215490.83741215650.83741215481.9×10−91.6×10−101000.10.83741219690.83741226520.83741218978.2×10−88.6×10−91000.40.83741219690.83741220090.83741219644.7×10−96.4×10−101000.80.83741219690.83741219850.83741219681.9×10−91.6×10−10There is no error inc∗ since both numerical and theoretical values of c∗ are zero.In Table4, we show the global relative errors and the computational times of f^num. and f^extpl. for different values of Δτ and Δv. The 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. The global relative error in f^extpl. is defined in the same way.Table 4
Maximum relative errors and computational times off^num. and f^extpl. in the domain 0≤τ≤100 and 0≤v=σ≤0.8.
Step sizef^num.f^extpl.ΔτΔvMax. rel. err.Comp. timeMax. rel. err.Comp. time0.020.021.5×10−50.393.3×10−70.840.010.013.8×10−61.524.2×10−83.310.0050.0059.6×10−76.025.2×10−913.180.00250.00252.4×10−724.026.3×10−1052.54Table4 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. 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 Table1. Extensive sets are given in Table 5.Table 5
Values for the parametersκ,θ,ξ,ρ,A,r,β,γ,ac, and aw.
ParameterκθξρArβγacawSet 12.50.04640.3796−76.70%1.551%0.06211Set 21.60480.080.3796−76.70%1.551%0.06211Set 31.60480.04640.6−76.70%1.551%0.06211Set 41.60480.04640.3796−50%1.551%0.06211Set 51.60480.04640.3796−76.70%1.251%0.06211Set 61.60480.04640.3796−76.70%1.555%0.06211Set 71.60480.04640.3796−76.70%1.551%0.15211Set 81.60480.04640.3796−76.70%1.551%0.06511Set 91.60480.04640.3796−76.70%1.551%0.06221Set 101.60480.04640.3796−76.70%1.551%0.06212By 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.Table 6
Maximum relative error forf^num. and φ^num. with Δτ=0.01andΔv=0.01.
SetMaximum relative errorf^num.φ^num.Set 15.8×10−63.9×10−6Set 23.6×10−62.7×10−6Set 33.6×10−63.9×10−6Set 44.0×10−61.8×10−6Set 52.4×10−62.0×10−6Set 64.0×10−62.8×10−6Set 74.2×10−62.9×10−6Set 82.1×10−61.4×10−6Set 93.8×10−62.7×10−6Set 103.8×10−62.7×10−6Table 7
Maximum relative error forf^extpl. and φ^extpl. with Δτ=0.01andΔv=0.01.
SetMaximum relative errorf^extpl.φ^extpl.Set 11.2×10−79.3×10−8Set 25.0×10−84.3×10−8Set 39.5×10−85.9×10−8Set 45.9×10−83.1×10−8Set 53.3×10−83.1×10−8Set 65.6×10−84.5×10−8Set 75.8×10−84.7×10−8Set 82.5×10−82.02×10−8Set 95.4×10−84.4×10−8Set 105.4×10−84.4×10−8By taking half of the step size above, namely, choosingΔτ=0.005 and Δv=0.005, 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 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 third-order accuracy, respectively, for the extensive sets of parameters.Table 8
Maximum relative error forf^num. and φ^num. with Δτ=0.005andΔv=0.005.
SetMaximum relative errorf^num.φ^num.Set 11.4×10−69.7×10−7Set 29.1×10−76.7×10−7Set 39.0×10−79.8×10−7Set 41.0×10−64.6×10−7Set 56.0×10−75.0×10−7Set 610.0×10−76.9×10−7Set 71.0×10−67.2×10−7Set 85.1×10−73.5×10−7Set 99.6×10−76.7×10−7Set 109.6×10−76.7×10−7Table 9
Maximum relative error forf^extpl. and φ^extpl. with Δτ=0.005andΔv=0.005.
SetMaximum relative errorf^extpl.φ^extpl.Set 11.6×10−81.2×10−8Set 26.4×10−97.5×10−9Set 32.1×10−87.6×10−9Set 47.6×10−94.0×10−9Set 54.2×10−94.0×10−9Set 67.1×10−95.8×10−9Set 77.4×10−96.0×10−9Set 83.2×10−92.5×10−9Set 96.9×10−95.6×10−9Set 106.9×10−95.6×10−9


## 4.3. Evidence for Stability and Convergence

To provide the evidence for stability, we calculate the maximum relative errors off^num. for Δτ=Δv=0.04,0.02,0.01,0.005,0.0025, 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.Figure 2
Maximum relative errors for decreasing step size.To provide the evidence for convergence, leteΔτ,Δv be the maximum relative error; then, the convergence order is given by(56)convergence order=log2eΔτ,ΔveΔτ/2,Δv/2.In Table10, 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.Table 10
Maximum relative errors and convergence orders off^num. in the domain 0≤τ≤100 and 0≤v=σ≤0.8.
Step sizeMaximum relative errorsConvergence orderΔτΔv0.020.021.53×10−52.00.010.013.83×10−62.00.0050.0059.58×10−72.00.00250.00252.40×10−72.0


## 4.4. Numerical Results forα≠0

We have confirmed the accuracy of our numerical solutions forα=0. This guarantees that our numerical solution fextpl. for α≠0 will also have an accuracy of OΔτ3. In Tables 11–13, we present the numerical solutions of fextpl., φextpl.∗, and cextpl.∗/w for α=0.1,0.5, and 0.9 with step sizes Δτ=Δv=0.001. In Tables 14 and 15, we show the results for γ=1 and 10 with α=0.5 and Δτ=Δv=0.001. 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 Δv. Therefore, we have provided the first set of exact numerical data for optimal asset allocation and consumption strategies.Table 11
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.1 and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.10.9765644460.7830377440.3238165870.10.40.9725289550.7830373940.3251602580.10.80.9597282710.7830362680.32949718710.11.2174435710.8175056160.25974737010.41.1910905830.8174484280.26549430510.81.1107573570.8172604090.284695630100.13.0734957700.8330103060.102888629100.42.9830468230.8329476710.106008315100.82.7117159650.8327395490.1166153721000.16.2276929850.8352401450.0507776741000.46.0374350620.8352066940.0523778331000.85.4674510730.8350951010.057838243Table 12
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.5 and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.10.7749377740.7828032100.9124691100.10.40.7718288440.7828022690.9161445400.10.80.7619663280.7827992380.92800266310.11.3689741000.8128728950.51652312610.41.3425532930.8127817070.52668805410.81.2617877680.8124844520.560400726100.16.0151811760.8323826420.117553696100.45.8400613900.8323127310.121078656100.85.3145274220.8320806890.1330516761000.113.9168104220.8352387810.0508095431000.413.4916571570.8352053110.0524106691000.812.2179612210.8350936530.057874368Table 13
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.9 and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.10.4096863870.7822441532.3156329550.10.40.4081604550.7822420302.3242900850.10.80.4033185340.7822352052.35219365810.11.2314519090.8082916300.77037786910.41.2105426960.8081887210.78368429410.81.1464004990.8078556960.827532175100.17.6867898130.8321316570.123417359100.47.4639749660.8320589530.127101618100.86.7952020830.8318177430.1396107561000.118.6674543490.8352383260.0508201751000.418.0971749900.8352048500.0524216241000.816.3887024200.8350931700.057886419Table 14
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.5, γ=1, and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.10.5468592821.5500000000.9143119940.10.40.5468592821.5500000000.9143119940.10.80.5468592821.5500000000.91431199410.10.9561778201.5500000000.52291528810.40.9561778201.5500000000.52291528810.80.9561778201.5500000000.522915288100.14.0343088511.5500000000.123936966100.44.0343088511.5500000000.123936966100.84.0343088511.5500000000.1239369661000.18.3139164411.5500000000.0601401281000.48.3139164411.5500000000.0601401281000.88.3139164411.5500000000.060140128Table 15
Numerical solutions forf,φ∗,andc∗/w after performing Richardson’s extrapolation. Here, α=0.5, γ=10, and Δτ=Δv=0.001. Other parameter values are the same as those in Table 1.
τσ=vfextpl.φextpl.∗cextpl.∗/w0.10.11.0247308490.1578338050.9105151780.10.41.0232359120.1578336810.9118454320.10.81.0184671920.1578332820.91611492210.11.8407898580.1695443640.50686556510.41.8270546360.1695304530.51067602110.81.7838731460.1694856200.523037747100.19.1822547040.1785759880.101612624100.49.0714269140.1785644610.102854049100.88.7260391460.1785270170.1069251441000.138.5460734750.1801360060.0242056561000.438.0501550490.1801328720.0245211351000.836.5060889640.1801226620.025558284


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



---

*Source: 9548060-2020-07-30.xml*