Optimal investment strategy post retirement without ruin possibility: A numerical algorithm

Journal of Computational and Applied Mathematics 363 (2020) 325–336

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Optimal investment strategy post retirement without ruin possibility: A numerical algorithm Hassan Dadashi Dept. Math. Institute for Advanced Studies in Basic Sciences (IASBS), Zanjan, P.O. Box 45195-1159, Iran

article

info

Article history: Received 1 August 2017 Received in revised form 21 October 2018 MSC: 60J70 93E20 65N06 Keywords: Defined contribution plan Decumulation phase HJB equation Minimum guarantee Finite difference method

a b s t r a c t We give a numerical algorithm based on the policy iteration method to find approximations of solution of an HJB equation arising from the portfolio optimization problem post retirement for a member of a defined contribution plan. In our numerical algorithm, we apply a fully implicit finite difference scheme to the HJB equation on a bounded domain. In formulating the optimal control problem, the loss function is written using a target function during the decumulation phase and at the terminal time. Following Di Giacinto et al. (2014), a minimum guarantee for the final annuity is considered. Finally, we give the simulation results for the final annuity, optimal investment strategy and optimal wealth process, for different levels of risk aversion and different loss functions. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Since in defined contribution (DC) plans the financial risks are borne mostly by the plan members, studying the portfolio optimization problem for individual members along with considering their own risk aversion is very important in accumulation (or preretirement) and decumulation (or postretirement) phase of the plan. At the retirement time the retiree faces to the annuity risk, the risk that exists at the time of low bond yield rate. In some DC plans, the retiree has the option to postpone the annuity purchase. Therefore, while drawing from the fund for their consumption, invests in the risky assets to get greater wealth for annuitization. In this work, focusing on the decumulation phase, we fix the annuitization time and investigate the optimal investment strategy. We employ the loss function applied in [1,2] and [3], in which a target-based quadratic loss function is defined totally by the wealth process. In fact, a target for the final wealth, or equivalently for the final annuity, at the terminal time is designed which consequently yields a target for the wealth amount at any prior time. Gerrard et al. [2] study the portfolio optimization problem post retirement when the consumption rate is fixed. In [3], they consider the consumption rate as a control variable too. Di Giacinto et al. [1] add to this framework a constraint on the final wealth to have a minimum guarantee on the final annuity. They study the HJB equation, obtained via the dynamic programming principle, by using the corresponding dual equation. They find a closed form for the investment strategy in the case that the loss function does not contain the running cost term. In this paper, assuming the same constraint on the wealth process, we consider the general form of the loss function, the loss function with running cost term, and give a numerical algorithm to get the approximations of the optimal strategies. Our simulation results show that considering the running cost in the E-mail address: [email protected]. https://doi.org/10.1016/j.cam.2019.02.027 0377-0427/© 2019 Elsevier B.V. All rights reserved.

326

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

value function formulation yields higher mean for the final annuity and keeps the wealth process relatively away from the minimum guarantee. In our numerical algorithm, we apply the fully implicit time stepping finite difference method. Then, in each time step, the policy iteration method is employed to tackle the nonlinearity appears in the algebraic equation. See [4,5] and [6] for some applications of iteration policy method in approximating the solutions of optimal control problems. It seems that this method is a powerful method with relatively low computational costs in approximating the solutions of HJB equations. Our aim in future works is to apply it in more developed frameworks and market models. Besides the different constraints on the control variables and the wealth dynamics, the different utility or loss functions that are employed in formulating the value functions yield quite different optimal control problems. We mention here a few related works. Milevsky and Young [7] consider the utility functions defined on the consumption and the final annuity and study the optimal annuitization problem for time dependent mortality function in all or nothing risky investment framework. In the framework of all or nothing risky investment and for different utility functions defined on the consumption and the final annuity, Stabile [8] studies the optimal investment/consumption problem and investigates the optimal time for purchasing the annuity subject to a constant force of mortality. Blake et al. [9] compare purchasing the annuity at the retirement age with distribution programs involving differing exposures to equities during retirement. Emms and Haberman [10] investigate the optimal annuitization time and determine the optimal investment strategy by minimizing the expected loss function. Albrecht and Maurer [11] compare the immediate annuitization and the income drawdown option. Moreover, they determine the probability of running out of money before the uncertain date of death. Gerrard et al. [12], by considering a target for the consumptions in decumulation phase and a target for the final wealth, investigate the optimal annuitization time together with the optimal consumption and investment rate. Vigna [13] shows that the target based approach yields portfolios that are on the efficient frontier in the mean–variance setting. The article is organized as follows. In the following section, the market model, the loss function and also the set of admissible strategies are specified. In the third section, after identifying the value function, we introduce the safety level for the wealth process. Furthermore, the HJB equation obtained via the dynamic programming principle is written and a change of variable that transforms the domain to a rectangle is introduced. In Section 4, the main achievement of this work, the numerical algorithm in finding approximations of solution of the HJB equation is presented. Finally in Section 5, the simulation results for the final annuity, optimal investment strategy and optimal wealth process are revealed. 2. The market model We consider a market consisting of a risky and a riskless asset with the following dynamics: dSt = St (µdt + σ dBt ), dAt = rAt dt , where B(·) is a Brownian motion on the filtered probability space (Ω , F , F, P). So the risky asset follows the geometric Brownian motion with constant volatility σ and expected return µ = r + σ β , in which r is the constant interest rate and β is the Sharpe ratio of the asset. Let the consumption rate during the decumulation phase, which is denoted by the time interval [0, T ], be the constant b0 . Moreover, at any time t, π (t) and 1 − π (t) denote the proportions of the fund’s portfolio that are invested in the risky and in the riskless asset, respectively. So, the dynamics of the wealth process (or the fund value) are given by dXt = {[π (t)(µ−r) + r ]Xt − b0 }dt + σ π (t)Xt dBt ,

(2.1)

X (0) = x0 . The control variable π is supposed to be chosen from the interval [0, 10]. This means that the short-selling is not permitted, but the portfolio manager can borrow from the money market up to ten times the current fund value. The above model can be extended in several directions. Among many other models, Boulier et al. [14] consider a model with stochastic interest rate. Han and Hung [15] equip the market model with the inflation rate which has important consequences on the optimal investment strategy. Hainaut and Deelstra [16] investigate the optimal time for annuity purchase in a market model with jump-diffusion dynamics when the economic utility function is replaced by the expected present value operator. We assume that the retirement age is a0 and the annuity purchase is postponed until the age a1 = a0 + T . The main purpose of a retiree in postponing the annuity purchase is to reach a desired annuity. Let F be the target for the final wealth by which the retiree can purchase the desired annuity at the terminal time T . In other words, if aa1 is the actuarial value of the unitary lifetime annuity at the age a1 , then F /aa1 would be the desired annuity. At any time t ∈ [0, T ], the wealth amount F (t) =

b0 r

+ (F −

b0 r

)e−r(T −t) ,

(2.2)

guarantees reaching the desired target F at the terminal time T, by investing the whole portfolio during the time interval [t , T ] in the riskless asset, the null strategy. Actually, this can be proved easily by applying the null strategy in the dynamics

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

327

(2.1), see [1, Lemma 3.2]. Therefore, the amounts {F (t), 0 ≤ t ≤ T } are regarded as the targets during the decumulation phase. We consider the square of distances of the wealth process from the targets {F (t), 0 ≤ t ≤ T } as one term of loss function, which is called the running cost term, and allot the weight κ to this term. Moreover, since the annuity purchase is at the terminal time T , we consider the square of the distance of the final accumulated wealth from the target F , as another term of the loss function. So, supposing the mortality rate δ to be constant and independent of the asset dynamics, the loss function is written as T

∫

κ e−ρ t (F (t) − X (t))2 dt + e−ρ T (F − X (T ))2 ,

(2.3)

0

where ρ = r + δ is the interest rate plus the mortality rate. For any 0 ≤ t ≤ T , consider the filter Ft := (Fst )s∈[t ,T ] where Fst is the σ -algebra generated by the random variables (B(u) − B(t))u∈[t ,s] . Then at any time t ≥ 0, the strategy π (·) is chosen from L2 (Ω × [t , T ]; [0, 10]), which is supposed to be Ft -progressively measurable. It should be noted that for the given process π (·), Eq. (2.1) has a unique strong solution, see [17, Section 5.6.C]. We denote this solution by X (·; t , x, π (·)). Since the loss function is defined totally by the wealth process, designing a minimum guarantee for the accumulated wealth is computationally tractable. If we consider S as the minimum guarantee at the terminal time then the set of admissible investment strategies reduces to

˜ ad (t) := {π ∈ L2 (Ω × [t , T ]; [0, 10])|π (·) is Ft − prog .meas., Π X (T ; t , x, π (·)) ≥ S a.s.}

(2.4)

3. The HJB equation

˜ ad (t) For any (t , x) ∈ [0, T ] × R+ , we define the following objective functional on Π ˜ , x; π (·)) := E[κ J(t

T

∫

e−ρ s (F (s) − X (s; t , x, π (·)))2 ds t

+ e−ρ T (F − X (T ; t , x, π (·)))2 ].

(3.1)

Our goal is to find the admissible strategies that minimize this functional. To solve this stochastic optimal control problem via the dynamic programming, we define the value function V˜ (t , x) :=

inf

˜ ad (t) π (·)∈Π

˜ , x; π (·)). J(t

(3.2)

The minimum guarantee constraint on the wealth amount at time T indicates some properties of the optimal strategy. Let us set for any t ∈ [0, T ] b0 b0 S(t) := − ( − S)e−r(T −t) . r r Then, the curve {S(t), 0 ≤ t ≤ T } is like an attractor for the wealth process. Actually, if for some 0 ≤ t < T the wealth process hits this curve, then the null strategy, π (s) = 0, s ≥ t, will be the optimal strategy and therefore, due to the wealth process dynamics, the optimal wealth process will remain on this curve from the time t onwards. Intuitively, if the accumulated wealth at any time t, X (t), becomes smaller than S(t), then by applying any admissible strategy there will be a positive probability of remaining under this curve until the terminal time T . For a rigorous proof of this statement see [1]. So, the assumption X (T ) ≥ S in (2.4) can be replaced by X (s) ≥ S(s), for all t ≤ s ≤ T ,

˜ ad (t) = {π ∈ L2 (Ω × [t , T ]; [0, 10])|π (·) is Ft − prog .meas., Π X (s; t , x, π (·)) ≥ S(s), t ≤ s ≤ T a.s.}. Moreover, as it was argued in Section 2, when the accumulated wealth equals to the target F (t), at any 0 ≤ t ≤ T , then the null strategy will keep the wealth process equals to the target function F (·) from the time t onwards, which yields zero cost. This means that the curve {F (t), 0 ≤ t ≤ T } is an attraction for the wealth process too. It should be noted that, as it is proved in [2], the wealth process will never hit the upper border, {F (t), 0 ≤ t ≤ T }, if the initial wealth x0 is less than F (0). Similarly, it can be proved that if x0 > S(0) then the optimal wealth process will not theoretically hit the bottom border. However, due to the approximation error made by the machine, the wealth process hits the bottom border in some realizations of our simulations. These observations indicate that if the initial value lies in [S(0), F (0)] then the wealth process will remain in a domain bounded at the top and at the bottom by the curves {F (t), 0 ≤ t ≤ T } and {S(t), 0 ≤ t ≤ T }, respectively. So, the optimal control problem leads to an HJB equation on a bounded domain. Moreover, notice that the problem without any constraint on the wealth process which yields an HJB equation on an unbounded domain, has explicit solution, see [2]. Due to the loss function formulation, (2.3), and the wealth process dynamics, (2.1), and by applying the dynamic programming principle we obtain the following HJB equation inf {

π∈R

∂ V˜ + A˜V˜ (t , x) + κ e−ρ t (F − x)2 } = 0, ∂t

(3.3)

328

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

where A˜ is the generator of the diffusion process (2.1), A˜ = ((π (t)[µ − r ] + r)x − b0 )

∂ 1 ∂2 + σ 2 π 2 (t)x2 2 . ∂x 2 ∂x

Additionally, the definition of V˜ implies the following Dirichlet boundary conditions V˜ (T , x) = e−ρ T (F − x)2 , V˜ (t , F (t)) = 0, V˜ (t , S(t)) = (F − S) ( 2

κ 2r

x ∈ [S , F ], t ∈ [0, T ], −2r(T −t)

(1 − e

)+e

−ρ T

),

(3.4)

t ∈ [0, T ].

Remark 3.1. Assuming the coefficients µ and σ to be time dependent deterministic functions, the structure of the generator of diffusion process X , and consequently the structure of the associated HJB equation, does not change. So, this assumption does not make great difficulties to our problem. However, for the simplicity of notations we assume that they are constants. Since we are going to apply the policy iteration method in finding an approximation of the solution, we consider the original HJB equation in contrast to Di Giacinto et al. [1], who study the corresponding dual equation. To apply the finite difference method more conveniently to the above PDE, we convert the domain of the equation to a rectangle, firstly introduced in [1]. Consider the diffeomorphism L : [0, T ] × [S , F ] → C , where C := {(t , x)|t ∈ [0, T ], S(t) ≤ x ≤ F (t)} and b0

(t , z) → (t , x) = L(t , z) = (t , L1 (t , z)) := (t , ze−r(T −t) +

r

(1 − e−r(T −t) )).

1 Notice that L1 (t , S) = S(t) and L1 (t , F ) = F (t). Actually, when the fund value at time t is x, or X (t) = x, then z := L− 1 (t , x) represents the value of the fund that one will have at time T by investing the whole portfolio in the riskless asset from time t onwards. Corresponding to the wealth process X (·), we define the process Z (·) as

Z (t) := [L1 (t , ·)]−1 (X (t)),

t ∈ [0, T ],

which gives b0

Z (t) = X (t)er(T −t) +

r

(1 − er(T −t) ),

t ∈ [0, T ].

Then, by applying the Itô formula, we get dZt = π (t)(µ − r)[Zt +

b0 r

(er(T −t) − 1)]dt + σ π (t)[Zt +

b0 r

(er(T −t) − 1)]dWt .

(3.5)

˜ ad (t), we set Z (·; t , z , π (·)) := [L1 (s, ·)]−1 (X (·; t , x, π (·))), in which z = Moreover, corresponding to any π (·) ∈ Π L1 (t , x). Then, the set of admissible strategies can be represented by the new dynamics as −1

Πad (t) = {π ∈ L2 (Ω × [t , T ]; [0, 10])|π (·) is Ft − prog .meas., S ≤ Z (s; t , z , π (·)) ≤ F , s ∈ [t , T ]}. Analogous to (3.1), for any (t , z) ∈ [0, T ] × [S , F ], we define the following functional on Πad (t) which is represented by the process Z (·; t , z , π (·)) T

∫

J(t , z , π (·)) := E{

κη(s)(F − Z (s))2 ds + η(T )(F − Z (T ))2 }, t

where η(s) = e−ρ s−2r(T −s) . Now, defining the value function V as V (t , z) :=

inf

π (·)∈Πad (t)

J(t , z ; π (·)),

we get the following HJB equation on the rectangular domain [0, T ] × [S , F ] inf {

π ∈Πad

∂V + AV (t , z) + κη(t)(F − z)2 } = 0, ∂t

(3.6)

in which A is the generator of the diffusion process Z (·), A =π (t)(µ − r)[z +

b0 r

(er(T −t) − 1)]

∂ 1 b0 ∂2 + σ 2 π 2 (t , z)[z + (er(T −t) − 1)]2 2 . ∂z 2 r ∂z

(3.7)

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

329

Furthermore, we easily obtain the following boundary conditions V (T , z) = e−ρ T (F − z)2 , V (t , F ) = 0, V (t , S) = (F − S)2 (e−ρ T

z ∈ [S , F ], t ∈ [0, T ], κ e−2rT (2r −ρ )T + (e − e(2r −ρ )t )), 2r − ρ

(3.8) t ∈ [0, T ].

Using the dual equation, Di Giacinto et al. [18] study the existence and uniqueness of regular solution of the above HJB equation. Although the existence of regular solution is a base for the verification theorem and obtaining the optimal feedback control, but the existence of a unique viscosity solution is enough to apply our numerical algorithm. Theorem 3.2.

Equation (3.6) with the boundary conditions (3.8) has a unique viscosity solution.

Proof. It is proved in [18, Proposition 4.8] that the value function V is continuous on the domain [0, T ] × [S , T ]. It can be checked easily that the coefficients of differential operator A and the running cost loss function, (F (t)−z)2 , have continuous partial derivatives with respect to the variables t and z on [0, T ] × [S , T ], for any fixed control variable. Furthermore, the control variable is chosen from a compact interval. So, [19, V. Corollary 3.1] implies the existence of viscosity solution for Eq. (3.6) with the boundary conditions (3.8). Moreover, the uniqueness of viscosity solution is implied by [19, V, Corollary 8.1]. □ 4. Numerical algorithm 1 We discretize the time interval [0, T ] to M = T × 52 subintervals with length ∆t = 52 , where each one corresponds to one week, and the space interval [S , F ] as S = z0 , z1 , . . . , zN +1 = F , in which the length of each subinterval equals to ∆z. To apply the finite difference method, we consider the forward scheme for the time and the space first derivative and the central difference scheme for the space second derivative. So, if we write V (i, j) = V (ti , zj ), for the value function at the node (ti , zj ), 0 ≤ i ≤ M , 0 ≤ j ≤ N + 1, the discretization of (3.6) at any node (ti , zj ), 0 ≤ i ≤ M − 1, 1 ≤ j ≤ N is written as

V (i + 1, j) − V (i, j)

∆t

+ inf {α (i, j)

V (i, j + 1) − V (i, j)

π∈Πad

+ β (i, j)

∆z V (i, j + 1) − 2V (i, j) + V (i, j − 1) (∆z)2

+ κη(ti )(F − zj ) } = 0, 2

(4.1)

in which

α (i, j) = π (ti )(ti , zj )(µ − r)[zj + β (i, j) =

1 2

σ 2 π 2 (ti )(ti , zj )[zj +

b0 r

b0 r

(er(T −ti ) − 1)],

(er(T −ti ) − 1)]2 .

The above discretization can be rewritten as V (i + 1, j) − ∆t inf { π ∈Πad

) α (i, j) β (i, j) +2 V (i, j) ∆t ∆z (∆z)2 ) ( α (i, j) β (i, j) β (i, j) + V (i, j + 1) − V (i, j − 1) − ∆z (∆z)2 (∆z)2

(

1

+

− κη(ti )(F − zj )2 } = 0

(4.2)

Since the coefficients α (i, j) and β (i, j) are positive at any node, the above representation indicates the positive coefficient property for our scheme, see [5, Condition 4.1]. Moreover, as mentioned in the proof of Theorem 3.2, the coefficients of operator A are Lipschitz continuous and the loss function is smooth enough to satisfy the assumptions of [19, V, Theorem 8.1], which guarantee the comparison property for Eq. (3.6). Now, if we prove that our scheme is consistent, monotone and stable in the l∞ norm, then [20, Theorem 2.1] implies the convergence of our numerical approximations to the viscosity solution of Eq. (3.6). The proof of monotonicity and consistency, when the scheme has the positive coefficient property, are quite straight, see [5, Lemmas 5.2 and 5.3]. Here, we give just the proof of stability. Proposition 4.1.

The discretization (4.1) has the l∞ - stability property. In fact, there is a constant C > 0 such that

∥V (i, ·)∥∞ ≤ C (F − S )2 ,

0 ≤ i ≤ M.

(4.3)

330

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

Proof. Inserting the optimal control π ∗ in the coefficients α (i, j) and β (i, j) of (4.1), we have for 0 ≤ i ≤ M − 1 and 1≤j≤N V (i, j) =V (i + 1, j) + ∆t α (i, j)

+ ∆t β (i, j)

V (i, j + 1) − V (i, j)

∆z V (i, j + 1) − 2V (i, j) + V (i, j − 1) (∆z)2

+ ∆t κη(ti )(F − zj ) . 2

(4.4)

So, we have

( ) ∆t ∆t |V (i, j)| 1 + α (i, j) + 2 β (i , j) ≤ ∆z (∆z)2 ( ) ∆t ∆t V (i + 1, j) + ∥V (i, ·)∥∞ α (i, j) + 2 β (i, j) ∆z (∆z)2 + ∆t κη(ti )(F − zj )2 .

(4.5)

If V (i, j1 ) = max1≤j≤N V (i, j) = ∥V (i, ·)∥∞ , then by rewriting the above relation for the node (i, j1 ), we get

) ( ∆t ∆t ∥V (i, ·)∥∞ 1 + α (i, j1 ) + 2 β (i , j ) ≤ 1 ∆z (∆z)2 ( ) ∆t ∆t ∥V (i + 1, ·)∥∞ + ∥V (i, ·)∥∞ α (i, j1 ) + 2 β (i , j ) 1 ∆z (∆z)2 2 + ∆t κη(ti )(F − zj1 ) , which means

∥V (i, ·)∥∞ ≤ ∥V (i + 1, ·)∥∞ + ∆t κη(ti )(F − S)2 .

(4.6)

Moreover, we have the bound e−ρ T (F − S)2 for the value function at the terminal time T . So, by going backward and considering (4.6), we get (4.3). Moreover, the bound is easily achieved on the upper and lower borders of the boundary. □ Remark 4.2. The infimum in Eq. (4.1) indicates that inserting any admissible strategy in the coefficients α (i, j) and β (i, j), equality in (4.4) turns to an inequality which preserves the inequality (4.5). So, the bound (4.3) holds uniformly over the set of admissible strategies. Since the coefficients α (i, j) and β (i, j) in (4.1) depend on the control variable π , the implicit time stepping method gives nonlinear algebraic equations at each time step. To get rid of this difficulty, at each time step, we apply the policy iteration method, in which the value function and the investment strategy are deriving iteratively to converge to the solution of (4.1). The value function is known at the terminal time T . So, we start from the last column t = T − ∆t inside the domain [0, T ] × [S , F ] and go backward in time. On each column t = ti , as a starting point of the value function, we take the vector V (ti + ∆t , ·), the value function at the next time step, and apply the following algorithm. Policy Iteration Algorithm: (I) Policy improvement: At any node (ti , zj ), 1 ≤ j ≤ N, solve the static optimization problem in (4.1) to find new investment strategy π new . (II) Policy evaluation: Given investment strategy π on the column t = ti and concerning the boundary conditions (3.8), solve Eq. (4.1) for the value function V on the nodes (ti , zj ), 1 ≤ j ≤ N. (III) End of iteration: If the following convergence criterion does not hold, return to Step (I) max ⏐V new (ti , zj ) − V old (ti , zj )⏐ ≤ {max ⏐V new (ti , zj )⏐} × 10−6 .

⏐

j

⏐

⏐

⏐

j

It should be mentioned that considering the space steps equal to ∆z = 0.05, the criterion expressed in Step (III) satisfies in less than 10 iterations of the above algorithm, for the values of S , F that are considered in the next section. Moreover, in a Core-i7 CPU, the running times of each iteration for the values of S , F that correspond to high, medium and low level of risk aversion, which are specified in the next section, are about 10, 20 and 50 s, respectively. To prove the convergence of above algorithm, at first we write the discrete representation (4.2) in the matrix form. Our proof is a modification of the proof given in [21, Theorem 1]. Let V i = (V (i, j))0≤j≤N +1 be the vector that represents the value function on t = ti and Di be an (N + 2) × (N + 2) tridiagonal matrix in which the first and last rows are zero and

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

[Di V i ]j =

β (i, j) i V (j − 1) + (∆z)2

(

α (i, j) β (i, j) + V i (j + 1) − ∆z (∆z)2

)

(

α (i, j) β (i, j) +2 V i (j), ∆z (∆z)2

)

331

1 ≤ j ≤ N.

Then Eq. (4.2), for any 0 ≤ i ≤ M − 1, can be rewritten in the matrix form as

[I − ∆tDi ]V i = V i+1 + ∆tE i + Gi − Gi+1 ,

(4.7)

in which E i represents the last term in (4.2) and Gi is a vector with just one nonzero element which represents the boundary value of V at (ti , S). Notice that (I − ∆tDi ) is a tridiagonal matrix in which the diagonal entries are positive, the off-diagonal entries are negative and the row sums are all equal to 1. So, it is an M-matrix which implies that its inverse is a nonnegative matrix, (I − ∆tDi )−1 ≥ 0. Theorem 4.3. The policy iteration algorithm gives a monotone sequence of functions which converge to the unique value function. Proof. For a fixed 0 ≤ i ≤ M − 1, let W 0 = V i+1 and W k , k ≥ 1 be the vector obtained in Step (II) in the kth iteration as an approximation of V i . Then, from (4.7) we get

[I − ∆tDi,k (π k )]W k = V i+1 + ∆tE i,k (π k ) + Gi − Gi+1 , in which for any k ≥ 1

π k = argminπ k ∈Πad {Di,k (π k )W k−1 + E i,k (π k )}. i,k

Note that the matrix D and the vector E By a few manipulations, we get

i,k

(4.8)

are modified in each iteration due to the modification of strategy in Step (I).

[ ] [I − ∆tDi,k+1 ](W k+1 − W k ) = ∆t (Di,k+1 W k + E i,k+1 ) − (Di,k W k + E i,k ) .

(4.9)

The construction (4.8) implies that, given W k , by employing the matrix Di,k+1 and the vector E i,k+1 we get the minimum value for all elements of the vector (Di,k+1 W k + E i,k+1 ). Therefor we have for all 0 ≤ j ≤ N + 1 (Di,k+1 W k + E i,k+1 ) − (Di,k W k + E i,k )

[

] j

≤ 0.

Now, since (I − ∆tDi,k+1 )−1 ≥ 0, we conclude from (4.9) that the sequence of vectors W k is decreasing. Furthermore, from Remark 4.2, the sequence (W k )k≥1 is uniformly bounded and therefore it is convergent. Moreover, the construction process of W k implies that it converges to V i , the unique solution of (4.1), which is an approximation of the value function on t = ti . □ 5. Simulation results To compare our simulation results with the results obtained in the relevant works, we take the same market parameters as in [2] and [1]. So, we consider the interest rate r = 0.03 and the expected return and the volatility of the risky asset µ = 0.08 and σ = 0.15, respectively, which implies a Sharpe ratio equal to β = 0.33. Furthermore, we consider a retiree with age a0 = 60 and the initial wealth x0 = 100 and set the decumulation period equal to T = 15 years, which means a1 = 75. The constant consumption rate is set to be b0 = 6.22, which is assumed to be equal to the annuity purchasable at the retirement time. We give the simulation results for three different values of target and safety level which correspond to three different degrees of risk aversion. It is clear that someone with higher risk aversion chooses lower target level and higher safety level. So, like to the assumptions of [1], we associate to the high risk aversion class, the target level F = 1.5b0 a75 and the safety level S = 32 b0 a75 , in which a75 is the actuarial value of a unitary lifetime annuity issued to an individual aged 75. In other words, in this class the final annuity will be at least two third of b0 and moreover by reaching the target F the member can purchase an annuity which is 50% higher than b0 . To the medium and the low level of risk aversion classes, we associate the target levels F = 1.75b0 a75 , F = 2b0 a75 and the safety levels S = 0.5b0 a75 , S = 0, respectively. To simulate the optimal wealth process by using the dynamics (2.1), a same 10 000 generated stream of pseudo random numbers for different levels of risk aversion and different impact factors of the running cost is applied. By employing the simulated optimal wealth process for three different levels of risk aversion when κ = 1, we present the histograms of the final annuity, Figs. 1–3, and reveal some percentiles of the optimal investment strategy, Figs. 4–6, and the optimal wealth amount, Figs. 7–9, during the decumulation period. Figs. 1–9 indicate that when the risk aversion decreases the distribution of the final annuity, the optimal wealth amount and the optimal investment strategy become more and more spread. Moreover, in this case we get higher average of final annuity and optimal wealth amount. Particularly, in the class of low risk aversion, in most of the times the final annuity is greater than 10 and consequently is much greater than b0 = 6.22. To provide more information about the distribution of the final annuity that are of immediate use for the members of the plan, Table 1 is given. In this table, the mean, standard deviation and some percentiles of the final annuity for different levels of risk aversion and different impact factors of the running cost are reported. You can find in this table

332

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

Fig. 1. Final annuity (high risk aversion).

Fig. 2. Final annuity (medium risk aversion).

the probability of getting a final annuity greater than b0 which provides more insights for comparing the immediate annuitization and postponing the annuity purchase. Particularly, these probabilities indicate that in our model postponing the annuity purchase is highly beneficial for investors with low and medium levels of risk aversion. However, for investors with high level of risk aversion, it does not give valuable final annuity. Comparing the statistics of the distributions of the final annuities for different impact factors of the running cost, κ = 0, 0.5, 1, 2, in which the case κ = 0 corresponds to the results obtained in [1], we see a considerable difference between the results obtained in the case κ = 0 and the results of the other scenarios. Moreover, we see that for higher κ we get greater mean and of course greater deviation for the final annuity. Furthermore, Table 1 shows that when the parameter κ increases, the probability that the wealth process attracts to the safety level decreases. The reason is that for higher values of κ , the paths of the wealth process that are close to the target function are more granted by the value function. 6. Conclusions Dealing with a nonlinear HJB equation on a bounded domain arising from the optimal investment problem for a member of a defined contribution plan, we give a numerical algorithm based on the policy iteration method to find

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

Fig. 3. Final annuity (low risk aversion).

Fig. 4. Investment strategy (high risk aversion).

Fig. 5. Investment strategy (medium risk aversion).

333

334

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

Fig. 6. Investment strategy (low risk aversion).

Fig. 7. Optimal wealth (high risk aversion).

Fig. 8. Optimal wealth (medium risk aversion).

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

335

Fig. 9. Optimal wealth (low risk aversion).

Table 1 Distribution of final annuity.

a

High risk aversion

Medium risk aversion

Low risk aversion

κ = 0, κ = 0.5 κ = 1, κ = 2

κ = 0, κ = 0.5 κ = 1, κ = 2

κ = 0, κ = 0.5 κ = 1, κ = 2

Mean

5.71, 5.73 5.73, 5.73

7.43, 7.53 7.53, 7.54

9.43, 9.50 9.51, 9.52

st. dev.

1.76, 1.93 1.95, 1.96

2.74, 2.97 3.00, 3.03

3.46, 3.72 3.76, 3.78

min

4.14, 4.14 4.14, 4.14

3.11, 3.11 3.11, 3.11

0.00, 0.00 0.00, 0.00

5th perc.

4.146, 4.149 4.15, 4.15

3.111, 3.113 3.113, 3.114

0.01, 0.005 0.005, 0.005

25th perc.

4.151, 4.152 4.152, 4.152

4.78, 3.22 3.12, 3.12

8.38, 8.92 9.00, 9.06

50th perc.

4.76, 4.152 4.152, 4.152

8.39, 8.86 8.90, 8.94

10.79, 11.17 11.20, 11.22

75th perc.

7.38, 7.75 7.81, 7.85

9.81, 10.13 10.16, 10.17

11.78, 11.92 11.94, 11.94

95th perc.

8.78, 8.98 9.00, 9.01

10.58, 10.69 10.70, 10.71

12.24, 12.30 12.31, 12.31

max

9.30, 9.31 9.31, 9.31

10.87, 10.88 10.88, 10.88

12.43, 12.44 12.44, 12.44

P (FAa > b0 ) (%)

52.58, 45.94 44.52, 43.60

76.42, 72.78 72.16, 71.70

88.53, 87.58 87.50, 87.45

P (FA = min. guarantee) (%)

5.45, 2.91 2.47, 2.25

4.34, 2.55 2.35, 2.14

2.23, 1.65 1.61, 1.60

FA = Final annuity.

an approximation of the solution. The numerical results show that the policy iteration method gives a relatively fast algorithm in finding good approximations of solution of this class of HJB equations. Furthermore, the simulation results show that considering the running cost as one term of the loss function, we get annuities with greater means and of course with greater deviations. Furthermore, in this case the probability of touching the safety level decreases considerably. Our aim for the next work is to study the portfolio selection problem post retirement in a market model with jumpdiffusion dynamics. Moreover, as a straight extension to the present model, the consumption rate can be considered as a control variable too.

336

H. Dadashi / Journal of Computational and Applied Mathematics 363 (2020) 325–336

References [1] M. Di Giacinto, S. Federico, F. Gozzi, E. Vigna, Income drawdown option with minimum guarantee, European J. Oper. Res. 234 (3) (2014) 610–624. [2] R. Gerrard, S. Haberman, E. Vigna, Optimal investment choices post retirement in a defined contribution pension scheme, Insurance Math. Econom. 35 (2004) 321–342. [3] R. Gerrard, S. Haberman, E. Vigna, The management of decumulation risk in a defined contribution environment, N. Am. Actuar. J. 10 (2006) 84–110. [4] T. Chellathurai, T. Draviam, Dynamic portfolio selection with fixed and/or proportional transaction costs using non-singular stochastic optimal control theory, J. Econom. Dynam. Control 31 (2007) 2168–2195. [5] P. Forsyth, G. Labahn, Numerical methods for controlled Hamilton–Jacobi-Bellman PDEs in finance, J. Comput. Finance 11 (2) (2007) 1–44. [6] P. Forsyth, L. Wang, Maximal use of central differencing for Hamilton–Jacobi-Bellman PDEs in finance, SIAM J. Numer. Anal. 46 (3) (2008) 1580–1601. [7] M.A. Milevsky, V.R. Young, Annuitization and asset allocation, J. Econom. Dynam. Control 31 (9) (2007) 3138–3177. [8] G. Stabile, Optimal timing of annuity purchases: a combined stochastic control and optimal stopping problem, Int. J. Theor. Appl. Finance 9 (2) (2006) 151–170. [9] D. Blake, A.J.G. Cairns, K. Dowd, Pensionmetrics 2: stochastic pension plan design during the distribution phase, Insurance Math. Econom. 33 (1) (2003) 29–47. [10] P. Emms, S. Haberman, Income drawdown schemes for defined contribution pension plan, J. Risk Insurance 75 (3) (2008) 739–761. [11] P. Albrecht, P. Maurer, Self-annuitization, consumption shortfall in retirement and asset allocation: the annuity benchmark, J. Pension Econ. Finance 1 (3) (2002) 269–288. [12] R. Gerrard, B. Højgaard, E. Vigna, Choosing the optimal annuitization time post-retirement, Quant. Finance 12 (7) (2012) 1143–1159. [13] E. Vigna, On efficiency of mean-variance based portfolio selection in DC pension schemes, Quant. Finance 14 (2) (2014) 237–258. [14] J.F. Boulier, S. Huang, G. Taillard, Optimal management under stochastic interest rates: the case of a protected defined contribution pension fund, Insurance Math. Econom. 28 (2) (2001) 173–189. [15] N.-W. Han, M.-W. Hung, Optimal asset allocation for DC pension plans under inflation, Insurance Math. Econom. 51 (1) (2012) 172–181. [16] D. Hainaut, G. Deelstra, Optimal timing for annuitization, based on jump diffusion fund and stochastic mortality, J. Econom. Dynam. Control 44 (2014) 124–146. [17] I. Karatzas, S. Shreve, Methods of Mathematical Finance, Springer-Verlag, New York, 1998. [18] M. Di Giacinto, S. Federico, F. Gozzi, E. Vigna, Constrained portfolio choices in the decumulation phase of a pension plan, in: Carlo Alberto Notebooks No. 155, 2010. [19] W. Fleming, H. Soner, Controlled Markov Processes and Viscosity Solutions, second ed., Springer-Verlag, New York, 2006. [20] G. Barles, P. Souganidis, Convergence of approximation schemes for fully nonlinear equations, Asymptot. Anal. 4 (1991) 271–283. [21] D. Pooley, P. Forsyth, K. Vetzal, Numerical convergence properties of option pricing PDEs with uncertain volatility, IMA J. Numer. Anal. 23 (2003) 241–267.