Optimal Annuitization and Investment
Problems
Anton Chtepine
Master’s Thesis to obtain the degree in Actuarial Science and Mathematical Finance University of Amsterdam
Faculty of Economics and Business Amsterdam School of Economics Author: Anton Chtepine Student nr: 10603034
Email: a.chtepine@gmail.com Date: August 14, 2014
Supervisor: Prof. Dr. M.H. Vellekoop Second reader: ir. J. de Kort
1 Introduction and Motivation 4
2 Objectives and Literature Overview 5
2.1 The Ruin Problem . . . 5
2.2 The Optimal Annuitization Problem . . . 7
2.2.1 Restricted market . . . 8
2.2.2 Unrestricted market . . . 9
2.3 Supplementary material . . . 10
3 The Ruin Problem 12 3.1 Reduction to a Diffusion Equation . . . 12
3.2 Finite-Difference Approximation . . . 13
3.3 Projected SOR method . . . 14
3.4 Solution Using PSOR . . . 15
3.5 The Method of Policy Iteration . . . 15
3.6 Discretizing the Initial Form of the Problem . . . 17
3.7 Solution Using Policy Iteration . . . 18
3.8 Numerical Experiments . . . 18
4 The Optimal Annuitization Problem 22 4.1 Restricted Market . . . 22
4.1.1 The Optimal Time of Annuitization . . . 22
4.1.2 Optimal Consumption and Investment Strategies . . . 23
4.1.3 Finite-Difference Approximation . . . 23
4.1.4 Solution via Policy Iteration in Matrix Form . . . 24
4.1.5 Policy Iteration - PSOR Hybrid Algorithm . . . 26
4.1.6 Numerical Experiments . . . 26
4.2 Unrestricted Market . . . 34
5 Conclusions 36
6 Bibliography 38
Chapter 1
Introduction and Motivation
As more and more individuals are becoming concerned with securing their welfare state post-retirement, it is now an absolute necessity to apply mathematical modeling in order to plan for the future. Such models aim to advise the retirees on, for instance, the optimal investment strategies, minimizing the probability of lifetime ruin (i.e. the probability of outliving one’s wealth) etc. But not all such model have analytic solutions. In fact, almost none of them do. This paper will attempt at developing a general numerical algorithm for solving the lifetime ruin probability minimization problem and the optimal annuitization problem presented by Milevsky, Moore and Young (2006) in [7] and Milevsky, Young (2007) in [9], respectively. The resulting algorithm will be based on the Policy Iteration scheme proposed by Howard (1960) [3]. Also, a second method will be presented and compared with the first one, which is a hybrid of Policy Iteration and the Projected Successive Over-Relaxation (PSOR) scheme.
The first problem deals with finding the optimal strategy of investing in a risky asset, so that the probability of lifetime ruin is minimized. It is assumed that the considered individual consumes at a fixed rate and has no bequest motives. We will describe in detail all the transformations and changes of variables done in [7] and, then, present both the algorithm proposed by Milevsky et al. (PSOR) and the method developed by us. Afterwards, several numerical experiments will be conducted in order to compare the results.
The assumptions of constant consumption rate and no bequest motives no longer hold for the second problem, where we aim at finding both optimal investment and consumption strategies in order to maximize the value function associated with the retiree. We will consider two cases, with and without bequest motives, and determine the optimal annuitization time in the former case, since, as proven by Yaari (1965) [16], in absence of bequest motives, it is optimal to annuitize completely. For this problem we will compare the exact solution obtained in [9] with the two numerical solutions we calculated using the plain Policy Iteration and its hybrid with the PSOR method. Afterwards, several numerical experiments will demonstrate the optimal annuitization times for different parameters along with a comparison of the value functions and investment and consumption strategies.
Our text is organized as follows. The next chapter will explain in detail the above mentioned problems. It will also provide a short summary of articles on similar subjects with comments on the differences between them. The first part of the third chapter will deal with a finite-difference approx-imation and solution of the transformed lifetime ruin probability minimization problem using PSOR. The second part will demonstrate the approximation of the problem’s initial form and explain the Policy Iteration scheme and its application. Numerical experiments will follow afterwards, in which we will compare the two methods. Chapter 4 is dedicated to the optimal annuitization problem. Similarly to the previous problem, we will present the finite-difference approximation, the two algorithms for obtaining the solution and several numerical experiments. In the last chapter the conclusions and recommendations regarding the choice of methods will be presented.
Along with Jin and Yin (2011) [5], this paper is among the first to attempt at finding the solutions of optimal annuitization and investment problems in a straightforward way, without having to resort to complex transformations of the initial form of the problem.
Objectives and Literature Overview
This chapter is divided into three sections. The first part is dedicated to the examination of an article by Milevsky, Moore and Young (2006) [7]. The second section deals with the article by Milesky and Young (2007) [9], which might be regarded as a modification of the first paper. The third section describes several articles on similar topics with some additions or different methods for obtaining solutions to our problems. It contains comments on the differences between all the works considered.
2.1
The Ruin Problem
The goal of [7] by Milevsky et al. is to obtain the probability of lifetime ruin and an optimal investment strategy as a solution to a variational inequality. This inequality is derived using optimal stochastic control and an optimization model presented later in the section.
It is assumed that an individual can invest into two types of assets: • riskless, whose price at time t follows a process for a fixed rate r ≥ 0:
(
dXt= rXtdt, X0 = x > 0
(2.1.1)
• risky, which follows the process at time t: (
dSt= St(µdt + σdBt), S0 = S > 0,
(2.1.2) where µ > 0, σ > 0 and Bt is a standard Brownian motion.
To account for the forces of mortality two functions are presented: λS(t) and λO(t). They are hazard rates for a specific individual and an objective hazard rate to price annuities. The objec-tive probability that an individual currently aged t lives s more years is represented by tps = exp(−Rt+s
t λO(u)du).Then, the function ¯a(t), ¯ a(t) =
Z ∞ 0
e−rstpsds, (2.1.3)
can be interpreted as the present value of an annuity paying $1 per unit of time to a person aged t. The form of ¯a(t) is easily explained by analyzing the integral: for a person aged t the present value of the annuity payment each year s (assuming the compounding is continuous) is equal to e−rs, and the probability that he survives until year s is tps. This way the integral is nothing more than the mathematical expectation of the total annuity payments to an individual.
Let As− and As denote the income rate of an individual at time s before and after any annuity purchase, respectively, and πs− — the amount of wealth invested in risky asset before any annuity
6 Anton Chtepine — Solving Optimal Annuitization and Investment Problems
purchase. The process for the individual’s wealth can, therefore, be written as: dWs= [rWs−+ (µ − r)πs−− Xs−]ds + σπs−dBs+ ¯a(s)dXs, Wt= w ≥ 0, Xt= x ≥ 0, (2.1.4)
where Xs = c − As can be viewed as excess consumption. We assume that an individual consumes a fixed and given amount c per year until he dies or his wealth reaches zero, the latter of which we will call “ruin”. Note, that the problem combines both continuous (πs) and singular (Xs = c − As) controls. To prove that the problem of choosing optimal excess consumption is on of singular control, we first assume the opposite, i.e. suppose that it is optimal to purchase annuities continuously. This means we can write dXs = χsds for some rate χs. Therefore, from 2.1.4 it follows that the resulting HJB equation would contain a term of the form χ¯a(t)ψx(w, x, t), which is maximized by setting χ = 0 or χ = ∞, depending on the sign ψx. Both of these rates are unacceptable, thus, the problem is of singular control.
Let the random time of death be denoted by τdand the time of ruin by τ0. With this notation the optimized probability of lifetime ruin defined on the set {(w, x, t) : 0 ≤ w ≤ x¯a(t), x ≥ 0, t ≥ 0} has the form
ψ(w, x, t) = inf πs,Xs
Pr[τ0 < τd|Wt= w, Xt= x, τd> t, τ0 > t]. (2.1.5) The reason for considering values of wealth under ¯a(t) will be stated later in the section.
Using Itˆo’s lemma for ψ in case it is optimal not to annuitize at all and in case it it optimal to annuitize instantaneously at the point (w, x, t), and combining the results, we come to the associated Hamilton-Jacobi-Bellman variational inequality
max λS(t)ψ − ψt− (rw − x)ψw− min π (µ − r)πψw+ 1 2σ 2π2ψ ww , ¯a(t)ψw+ ψx = 0, (2.1.6) where ψt= ∂ψ ∂t, ψw= ∂ψ ∂w, ψww= ∂2ψ ∂2w. (2.1.7)
By observation of Milevsky and Robinson (2000) in [8] the probability of lifetime ruin does not depend on w and x separately, but only on their ratio z = w/x. Therefore, the dimensions of the minimization problem can be reduced and a proposition is proved, stating that the probability of lifetime ruin is given by ψ(w, x, t) = V (z, t) if z = w/x < ¯a(t) and ψ(w, x, t) = 0 otherwise, where V (z, t) solves the following equation:
λS(t)V = Vt+ (rz − 1)Vz+ min π (µ − r)πVz+ 1 2σ 2π2V zz , (2.1.8)
for z < ¯a(t), with boundary conditions V (0, t) = 1 and V (¯a(t), t) = 0. If w ≥ x¯a(t) an individual can immediately buy an annuity guaranteeing him a yearly income of x = c − A, which in combination with his existing income of A allows him to have income exactly equal to the consumption rate c. Therefore, when z = w/x ≥ ¯a(t) the probability of lifetime ruin is 0.
The minimization over π can be calculated directly by setting the derivative of the minimized function w.r.t. π equal to zero. This way the minimizing value of π and the minπ term are found to be
π∗ = −(µ − r)Vz σ2V zz , (2.1.9) min π (µ − r)πVz+ 1 2σ 2π2V zz = −(µ − r) 2V2 z 2σ2V zz . (2.1.10)
Linearization is conducted in two steps. Firstly, to get rid of the term λS(t)V , define f (z, t) = exp − Z t 0 λS(u)du V (z, t). (2.1.11)
It follows that 2.1.8 becomes ft+ (rz − 1)fz−(µ − r) 2f2 z 2σ2f zz = 0. (2.1.12)
Then, we switch from f (z, t) to its Legendre dual ˜f (y, t) = minz>0[f (z, t) + zy]. Assuming differentia-bility, the derivative of f (z, t) + zy w.r.t. z equals fz(z, t) + y, so the optimal z∗ solves fz(z, t) + y = 0. Therefore, one can write z∗ = I(−y, t) where I(−y, t) is the inverse of f
z(z, t) w.r.t. z. It can be checked that
˜
f (y, t) = f (I(−y, t), t) + yI(−y, t), (2.1.13) ˜
fy(y, t) = I(−y, t). (2.1.14)
One can show that the function f can be retrieved from ˜f by the relationship f (z, t) = max
y>0[ ˜f (y, t) − zy]. (2.1.15) Analogously to the previous paragraph, the critical value y∗ solves ˜f
y(y, t) − z = 0. Therefore, due to 2.1.13-2.1.14, y∗= −fz(z, t) and ˜f (y∗, t) − zy∗ = f (z, t). The necessary derivatives are given by
˜ fyy= −
1 fzz(I(−y, t), t)
, f˜t= ft(I(−y, t), t). (2.1.16)
To get the final form of the problem rewrite 2.1.12 in terms of ˜f : ˜
ft(y, t) − ry ˜fy(y, t) + my2f˜yy(y, t) + y = 0, (2.1.17) in which m = 12 µ−rσ 2
. The boundary conditions are implicit and given by ˜
f (y0(t), t) = e− Rt
0λS(u)du, for ˜fy(y0(t), t) = 0, (2.1.18)
˜
f (yb(t), t) = ¯a(t)yb(t), for ˜fy(yb(t)), t) = ¯a(t). (2.1.19) The algorithm for solving2.1.17 is presented in chapter3.
2.2
The Optimal Annuitization Problem
In this section we will consider a model similar to the one in the previous section, but more complex. Milevsky and Young (2007) in [9] consider a process identical to2.1.4. Their article considers restricted and unrestricted markets, i.e. markets where an individual annuitizes all or nothing and anything anytime. We will divide this section into two subsections dealing with the former and the latter case, respectively.
Unlike the previous problem, here we do not aim at calculating the probability of financial ruin, but instead we calculate the optimal time of annuitization and the utility-maximizing strategies of consumption and investment in risky asset of a retiree facing a stochastic time of death.
In both cases the attention is restricted to utility functions with constant relative risk-aversion (CRRA), i.e. u(c) = ( c1−γ 1−γ, γ > 0, γ 6= 1, ln(c), γ = 1, (2.2.1) for c > 0.
8 Anton Chtepine — Solving Optimal Annuitization and Investment Problems
2.2.1 Restricted market
In case of restricted markets it is proved, for example, by Yaari (1965), that in absence of bequest motives the optimal way to annuitize is to annuitize completely. So, we assume that an individual annuitizes his entire wealth Wτ at a stochastic time τ determined by the agent and consumes Wτ/¯aOx+τ afterwards. With these assumptions the value function associated with our problem has a form
U (w, t) = sup {cs,πs,τ } Ew,t Z τ t e−r(s−t)s−tpSx+tu(cs)ds + Z ∞ τ e−r(s−t)s−tpSx+tu Wτ ¯ aO x+τ ds = sup {cs,πs,τ } Ew,t Z τ t e−r(s−t)s−tpSx+tu(cs)ds + e−r(τ −t)τ −tpSx+tu Wτ ¯ aO x+τ ¯ aSx+τ . (2.2.2) In a way analogous to the previous section, the value function is found as a smooth solution to a variational inequality similar to 2.1.8, but with an extra term maxc[u(c) − cVw]. This with the restraining condition V (w, t) ¯ aS x+t ≥ u w ¯ aO x+t (2.2.3) results in the following variational inequality:
(λSx+t+ r)V ≥ Vt+ rwVw− min π (r − µ)πVw−1 2σ 2π2V ww − minc [cVw− u(c)], (2.2.4) V (w, t) ¯ aS x+t ≥ u w ¯ aO x+t (2.2.5) with equation in at least one of them. The term w/¯aOx+tin the restriction2.2.3 can be regarded as the annuity income, i.e. the amount an individual will consume annually after annuitizing his wealth w, calculated using the actual market pricing hazard rates. Therefore, the whole constraint means that the value of one’s personal and subjectively assessed annuity income should be at least that of the objective annuity income.
With the use of utility consumption the solution of this problem does not require the complex transformations of the previous section. It is shown that the solution is obtained simply by considering a special form of V (w, t), namely V (w, t) = 1/(1−γ)w1−γψγ(t). As we can see, w and t are multiplicable separable, therefore, the optimal time to annuitize is independent of wealth and, thus, deterministic. So, in this case both the value function and the optimal time are given explicitly. We assume that the time an individual annuitizes is T . We can solve this problem for any value of T and then optimize over this value. Let φ = φ(t; T ) be the solution. For t < T the function φ has the form
φ(t; T ) = ¯a S x+T (¯aO x+T)1−γ !1/γ e−((r−δ(1−γ)/γ))(T −t)(T −tpSx+t)1/γ (2.2.6) + Z T t e−((r−δ(1−γ)/γ))(s−t)(s−tpSx+t)1/γds, (2.2.7) with δ = r + (((µ − r)/σ)2)/2γ. For t ≥ T φ(t; T ) = ¯a S x+T (¯aOx+T)1−γ !1/γ (2.2.8)
Let ¯U (w, t; T ) = 1/(1−γ)w1−γφγ(t; T ). The optimal time to annuitize is obtained by differentiating ¯ U (w, t; T ) with respect to T, t < T : ∂ ¯U ∂T ∝ γ 1 − γ ¯ aSx+T ¯ aOx+T !(γ−1)/γ − 1 1 − γ + ¯ aSx+T ¯ aOx+T + ¯aSx+T[δ − (r + λOx+T)]. (2.2.9)
If there exist such time T′ > 0 that the right-hand side above changes sign from positive to negative at T′, then it is optimal to annuitize at T′, and U (w, t) = ¯U (w, t; T ). If the expression is negative for all T ≥ 0, one should annuitize immediately, and U(w, t) = ¯U (w, t; 0).
In the later chapters we will prove that the optimal values of π and c corresponding to the solution ¯
U (w, t; T ) are, in fact, linear in w.
2.2.2 Unrestricted market
In case of unrestricted markets the trick of guessing the solution does not work. Milevsky and Young (2007) in [9] show that in this case the individual’s optimal annuitization strategy is defined by a barrier policy and he will pick moments and amounts to purchase annuities in order to stay on one side of the barrier in wealth-annuity space. An individual is allowed to choose between lump-sum and continuous annuity purchasing, whichever is optimal. Here the associated value function U (w, t) is defined using the strictly increasing, concave utility functions of consumption and bequest, u1(c) and u2(w), respectively: U (w, A, t) = sup {cs,πs,As} Ew,A,t Z ∞ t e−r(s−t)s−tpSx+t{u1(cs) + λSx+su2(Ws−)}ds . (2.2.10) The new variable A denotes the annuity income rate, Ew,A,tis the expectation conditional on Wt− = w and At− = A. The value function U is proven to be jointly concave in w and A by Milevsky and Young in [10].
For this case the variational inequality is similar to2.2.4-2.2.5 but with extra term λSx+tu2(w): (λSx+t+ r)U ≥ Ut+ (rw + A)Uw− min
π (r − µ)πUw− 1 2σ 2π2U ww (2.2.11) − min c≥0[cUw− u(c)] − λ S x+tu2(w), ¯ aOx+tUw ≥ UA. (2.2.12)
Note that the marginal utility of annuity income can be viewed as the marginal utility of benefit, while the adjusted marginal utility of wealth is the utility of cost. Similarly to many results in economics, the constraint means that the utility of cost is at least that of later income.
The authors consider the same utility of consumption 2.2.1 and the utility of bequest u2(w) = ku1(w), k ≥ 0. For utility functions chosen in such way, the value function U is homogenous of degree 1 − γ w.r.t. w and A, i.e. U(αw, αA, t) = α1−γU (w, A, t) for α ≥ 0. Therefore, U(w, A, t) = A1−γV (w/A, t) where V (z, t) = U (z, 1, t). With this notation the variational inequality for V has the form (λSx+t+ r)V ≥ Vt+ (rz + 1)Vz− min ˆ π (r − µ)ˆπVz− 1 2σ 2πˆ2V zz (2.2.13) − min ˆ c≥0[ˆcVz− u1(ˆc)] − kλ S x+tu1(z), (z + ¯aOx+t)Vz ≥ (1 − γ)V, (2.2.14)
in which ˆπ = π/A and ˆc = c/A.
Milevsky and Young prove that for each time t ≥ 0 there exists z0(t) such that (z0(t)+¯aox+t)Vz(z0(t), t) = (1 − γ)V (z0(t), t), and if z = w/A > z0(t) the individual immediately buys an annuity so that V (z, t) = V (z0(t), t) and
z0(t) = w − ∆A¯a O x+t
a + ∆A . (2.2.15)
If z < z0(t) the individual does not annuitize at all, and V solves 2.2.13 as an equation. Thus, for every t ≥ 0 the barrier w = z0(t)A is a ray emanating form the origin and lies in the first quadrant of (w, A) space.
10 Anton Chtepine — Solving Optimal Annuitization and Investment Problems
2.3
Supplementary material
In [1] E. Bayraktar et al. propose to regard volatility σ as a stochastic variable given at time t by the function σt= f (Yt, Zt). They do not consider the possibility and necessity of annuitization. The article only deals with minimizing the probability of lifetime ruin by investing in the assets correctly. This way the probability of lifetime ruin has the form of ψ(w, y, z) = infπPrw,y,z[τ0 < τd], where τ0 is the time of ruin, τd is the time of death and y and z are explained below. The processes for the price of a risky asset, Ytand Zt are given by
dSt= St(µdt + σtdBt(1)), S0= S > 0, dYt= 1ǫ(m − Yt)dt + ν q 2 ǫdB (2) t , Y0 = y, dZt= δg(Zt)dt + √ δh(Zt)dBt(3), Z0 = z. (2.3.1)
where f is bounded, bounded away from zero, positive and smooth, and Y and Z are two diffusions. Y is a fast mean-reverting Gaussian Ornstein-Uhlenbeck process with 0 < ǫ ≪ 1 and Bt(2) — standard Brownian motion. Z is a slow diffusion process behind the risky asset’s price; g and h are smooth; 0 < δ ≪ 1, and Bt(3) is a standard Brownian motion. Let ρij denote the correlation coefficients between Bt(i) and Bt(j). With these correlation coefficients the authors construct three other standard Brownian motions ˜Bt(1), ˜Bt(2), ˜Bt(3), which are independent. They are used to define the controlled process
dXtγ= −(r − λ)Xtγdt − µ − r f (Yt, Zt)
XtγB˜t(1)+ γt(2)d ˜Bt(2)+ γt(3)d ˜Bt(3), X0 = x > 0, (2.3.2) in which γ = (γ(2), γ(3)) is the control. For x > 0 the function ˆψ is defined by
ˆ ψ(x, y, z) = inf τ supγ E x,y,z Z τ 0 e−λtcXtγdt + e−λτmin ((c/r)Xtγ, 1) , (2.3.3)
where c is the individual’s rate of consumption. Bayraktar et al. find ˆψ(x, y, z) as the solution to a free-boundary problem specified in the article and prove that the minimum probability of lifetime ruin is given by the Legendre transform of ˆψ(x, y, z) with respect to x. As the algorithm for solving the free-boundary problem the authors propose to asymptotically expand ˆψ in powers of√ǫ and √δ. For applicability of this method they refer to an article by Jonsson and Sircar (2002) [6]. Expanding ˆψ instead of ψ is motivated by the necessity of solving non-linear differential equations for each term in the expansion of ψ, while for the terms of ˆψ only linear differential equations have to be solved, and it is done explicitly. As a result, the proposed algorithm is very much analogous to the one developed in section 2.1. The difference is, of course, in extending the problem to the case of stochastic volatility and in the fact that they optimize only in π while our two problems require optimization in (π, c) and (π, c, τ ).
A different algorithm for solving the optimal control problem directly, the Markov Chain approxi-mation, has been discussed by both E. Bayraktar et al. [1] and Jin and Yin in [5]. The former authors consider a process completely identical to ours, containing continuous dynamics and discrete events. To approximate the diffusion process they construct a discrete-time, finite-state, controlled Markov chain, which is defined to be locally consistent with2.1.4. Their Markov chain approximation method is motivated by the fact that a direct discretization greatly depends on the properties of the variational inequalities. The proposed algorithm uses mainly probabilistic methods and does not need any analytic properties of the solutions of the system of variational inequalities. As a result the proposed algorithm is mathematically more complex but theoretically can be applied to a wider range of problems.
Another realistic feature is introduced in [2] by Gerrard, Højgaard and Vigna. It is assumed that up until annuitization the pensioner has three degrees of freedom: he can select a strategy for investing his own funds; he can withdraw any amount from his own funds at any time until annuitization (if any); he can select the time for annuitization (if ever). Two functions, U1(b) and U2(kx), are introduced to account for utilities from a payment of size b before annuitization and from the same payment after it. The force of mortality is assumed to be constant and equal to δ. The size of annuity that can
be purchased with x is kx, where k > r (r - interest rate on a riskless asset). As the authors state themselves, the main contribution of their work is the closed-form solution of the optimal annuitization problem when quadratic loss functions are present, as Vigna (2009) [13] proved that optimal portfolios obtained with CARA and CRRA utility functions are inefficient in the mean-variance setting. Gerrard et al. deal with the problem of choosing control variables πt (the fraction of the fund invested in the risky asset at time t), bt (total amount withdrawn from own funds between time t and annuitization) and a stopping time T , so that the following expectation is maximized:
Ex " v Z τ 0 e−(r+δ)t(b0− bt)2dt + w e−(r+δ)τ r + δ (b1− kxτ) 2 # , (2.3.4)
where τ = min(T, T0), T0 is the time of ruin, v and w are positive weights. With the quadratic loss function the authors are able to pose the constraints for the value function of the problem depending on the region where x lies, and construct a solution satisfying the constraints. Then, via the verification theorem proved in the article, they demonstrate that the constructed solution is, indeed, correct. This work ends with numerical experiments addressing the dependence of the solution on the parameters.
Chapter 3
The Ruin Problem
When attempting to solve the problem introduced in section 2.1 numerically by explicitly mini-mizing the minπ term, we inevitably encounter trouble created by taking differentiating w.r.t. π. This operation leads to the derivatives of V in inconvenient powers. Therefore, a stable numerical algorithm is required to solve this problem. In this chapter we will solve it using two different methods.
The first method, the Projected Successive Over-Relaxation algorithm, allows us to solve 2.1.8 after all the transforms and changes of variables, i.e. we will demonstrate step-by-step the numerical method for solving the following partial differential equation:
− ˆft(y, t) + ry ˆfy(y, t) − my2fˆyy(y, t) − y = 0, y ∈ (yb(t), y0(t)), (3.0.1) subject to the constraint
ˆ f (y, t) ≤ ˆu(y, t), y ∈ (yb(t), y0(t)), (3.0.2) where ˆ u(y, t) = min exp − Z t 0 λS(v)dv , ¯a(t)y . (3.0.3)
This equation can be turned into a diffusion equation via several successive changes of variables and then solved as described by Wilmott, Howison and Dewynne (1995) in [15]. We solve the constrained PDE on a domain containing the free boundary (yb(t), y0(t)) and then recover its location after com-puting the solution.
The second method, the Policy Iteration algorithm, can be applied to the initial equation without any transforms: λS(t)V = Vt+ (rz − 1)Vz+ min π (µ − r)πVz+ 1 2σ 2π2V zz . (3.0.4)
Unlike the free-boundary problem solved by the first method, this problem has defined boundary conditions V (0, t) = 1 and V (¯a(t), t) = 0. We will present the correct finite-difference approximation of our equation and the algorithm itself as described by Wang and Forsyth (2006) in [14]. We only consider the values of π ≥ 0, since we do not allow an individual to invest a negative amount.
Also, since if z > ¯a(t) the probability of lifetime ruin is zero, we will restrict our attention to values z ∈ [0, ¯a(t)] and set a certain high enough value as the upper bound for π, as well.
3.1
Reduction to a Diffusion Equation
To get rid of the y and y2 terms in3.0.1 we use the transformation y = eξ. The derivatives of ˆf are, therefore: ˆfy = e−ξfˆξ and ˆfyy= e−2ξ( ˆfξξ− ˆfξ). This results in the following equation:
− ˆft(ξ, t) + (r + m) ˆfξ(ξ, t) − m ˆfξξ(ξ, t) − eξ= 0. (3.1.1) 12
To turn it into a forward equation we set t = −τ/m, so ˆft= −m ˆfτ. Therefore, with k = −r/m, ˆ
fτ(ξ, τ ) = ˆfξξ(ξ, τ ) + (k − 1) ˆfξ(ξ, τ ) + eξ/m. (3.1.2) Let ˆf = Af , where A = eαξ+βτ. Then,
ˆ
fτ = A(βf + fτ), fˆξ = A(αf + fξ), fˆξξ= A(α2f + 2αfξ+ fξξ), (3.1.3) and the resulting equation is:
βf + fτ = α2f + 2αfξ+ fξξ+ (k − 1)(αf + fξ) + eξ/Am. (3.1.4) The constraining condition3.0.2 is, therefore, transformed into
f (ξ, τ ) ≤ u(ξ, τ), (3.1.5)
where u(ξ, τ ) = e−αξ−βτu(eˆ ξ, τ ).
By choosing α = (1 − k)/2 and β = −(k − 1)2/4 we eliminate the f and fξ terms to get the final form of our equation:
fτ = fξξ+ 1 me (k+1)ξ 2 + (k−1)2τ 4 , (3.1.6)
which can be rewritten as
fτ = fξξ+ Q(ξ, τ ). (3.1.7)
This is a diffusion equation with a replenishment (or consumption) term Q(ξ, τ ).
3.2
Finite-Difference Approximation
As stated in [15] explicit and implicit finite-difference schemes have several limitations. Firstly, both of them have O(∆τ ) rates of convergence to the solution of the PDE. Secondly, the explicit scheme is unstable for ∆τ /(∆ξ)2 > 1/2. To overcome these limitations we approximate fτ and fξξ using Crank-Nicolson method. This way we have a finite-difference scheme which converges to the solution at a rate O((∆τ )2) for any ∆τ /(∆ξ)2> 0.
Consider a regular mesh with step sizes ∆τ and ∆ξ, and N+ with N− suitably large for the following inequality to hold:
N−∆ξ ≤ ξn= n∆ξ ≤ N+∆ξ. (3.2.1)
According to the Crank-Nicolson method fτ and fξξ are approximated as follows: fτ(ξ, τ + ∆τ /2) = fk+1 n − fnk ∆τ + O (∆τ ) 2 , (3.2.2) fξξ(ξ, τ + ∆τ /2) = 1 2 fn+1k+1− 2fnk+1+ fn−1k+1 (∆ξ)2 ! +1 2 fk n+1− 2fnk+ fn−1k (∆ξ)2 ! + O (∆ξ)2 , (3.2.3)
where fnk = f (n∆ξ, k∆τ ). Thus, dropping the last terms of 3.2.2 and 3.2.3, 3.1.7 and 3.1.5 are ap-proximated the following way:
fnk+1−1 2α(f k+1 n+1− 2fnk+1+ fnk+11 ) = f k n+ 1 2α(f k n+1− 2fnk+ fn−1k ) + ∆τ Qkn, (3.2.4) fnk≤ ukn, (3.2.5)
14 Anton Chtepine — Solving Optimal Annuitization and Investment Problems where α = ∆τ /∆ξ2, Qk n= Q(n∆ξ, k∆τ ) and ukn= u(n∆ξ, k∆τ ). Define Znk by Znk= fnk+ 1 2α(f k n+1− 2fnk+ fn−1k ) + ∆τ Qkn. (3.2.6) Then,3.2.4 becomes (1 + α)fnk+1−1 2α(f k+1 n+1 + fn−1k+1) = Znk. (3.2.7) Let vectors fk, uk and bk be given by
fk=fNk−+1, · · · , fNk+−1 T , (3.2.8) uk=ukN−+1, · · · , ukN+−1 T , (3.2.9) bk= bk N−+1 .. . .. . bkN+−1 = Zk N−+1 .. . .. . ZNk+−1 +1 2α fNk+1− 0 .. . 0 fNk+1+ , (3.2.10)
with bk only containing information on time step k and boundary conditions for f on time step k + 1 which are stipulated in the next section. Using the introduced notation we can rewrite our problem in matrix form: Cfk+1= bk, fk+1 ≤ uk+1, (Cfk+1− bk)(fk+1− uk+1) = 0. (3.2.11) where C = 1 + α −α/2 0 · · · 0 −α/2 1 + α −α/2 ... 0 −α/2 . .. . .. 0 .. . . .. 1 + α −α/2 0 · · · 0 −α/2 1 + α . (3.2.12)
3.3
Projected SOR method
The SOR method is a refinement of a well-known Gauss-Seidel method, which is itself a refined Jacobi method. The idea behind time-stepping from fkto fk+1using the above mentioned algorithms is to approach fk+1with a specially constructed series {fk+1,l}, l = 0, 1, 2, . . . until the approximations fk+1,l stop changing significantly. When this happens fk+1,lis considered the solution fk+1.
The iterations for the SOR algorithm adapted for Crank-Nicolson finite-difference approximation have the following form:
yk+1,l+1n = 1 1 + α(b k n+ 1 2α(f k+1,l+1 n−1 + f k+1,l n+1 )), (3.3.1) fnk+1,l+1= fnk+1,l+ ω(ynk+1,l+1− fnk+1,l), (3.3.2) where 1 < ω < 2 is the over-relaxation parameter. As with all iterative methods, the initial approx-imation fk+1,0 plays a very important role. A good choice for fk+1,0 is usually fk, the value on the previous time step. Note, that the ynk+1,l+1 is constructed using information from both iteration l and l + 1 due to fn−1k+1,l+1 being already known at the moment of calculating the value at node n.
To ensure that3.2.5 is satisfied we modify the second equation:
The initial guess is also changed to fk+1,0 = max(fk, uk+1). The process 3.3.1, 3.3.3 is called the Projected SOR, and it is the easiest way to solve the problem 3.2.11 with the condition 3.2.5. The specified steps are iterated until the condition
||fk+1,l+1− fk+1,l|| ≤ ǫ (3.3.4) is satisfied, where ǫ is the required tolerance. After the tolerance is reached we put fk+1 = fk+1,l+1. For the numerical experiments we assume || . . . || to be Euclidean distance, i.e. for x = (x1, . . . , xn) and y = (y1, . . . , yn)
||x − y|| =p(x1− y1)2+ · · · + (xn− yn)2. (3.3.5) The constraint is enforced at the very time the iterate fnk+1,l+1 is computed, therefore, it directly affects the calculation of fn+1k+1,l+1, fn+2k+1,l+1 etc. Thus, there is an internal consistency in this method: we do not simply solve 3.2.11 and then apply the constraint 3.2.5 to the elements of the solution. Such actions are invalid in case of implicit finite-difference formulations of the problem because all the components fnk+1 of fk+1 are implicitly dependent on each other, so changing one without affecting the others is impossible.
As mentioned in [4], the convergence of the method depends on the choice of relaxation factor ω. Let ρ be the spectral radius of the Jacobi iteration matrix G = D−1(C − D), where D is the diagonal of C. Then the classical choice for ω, given the eigenvalues of G are real, is
ω = 2 1 +p1 − ρ2, ρ ≈ maxi 1 Ci,i X j6=i |Ci,j|. (3.3.6)
3.4
Solution Using PSOR
To obtain the solution ˆf to our problem the changes of variables made in section 3.1 have to be reversed. As stated in 3.1, ˆf (ξ, τ ) = A(ξ, τ )f (ξ, τ ) where A(ξ, τ ) = eαξ+βτ with α = (1 − k)/2 and β = −(k − 1)2/4. Then we switch from (ξ, τ ) back to (y, t) via the transforms y = eξ, τ = −mt.
The boundary points yb(t) and y0(t) are the points where inequality in3.0.2 changes to equality. It is proved by Moore and Young (2006) in [11] that yb(t) = 0, i.e. the unknown left boundary is zero, which is useful for numerical solution. Therefore, we only have to retrieve the right boundary y0(t) using2.1.18,2.1.19.
Our numerical solution ˆf solves the free-boundary problem 2.1.17, 2.1.18 and 2.1.19. Thus, we retrieve the function f (z, t) from ˆf (y, t) via the transform 2.1.15:
f (z, t) = max
y>0[ ˆf (y, t) − zy]. (3.4.1) Having obtained f (z, t) we find V (z, t) = eR0tλS(u)duf (z, t), and using the proposition from [7]
we have the probability of lifetime ruin ψ(w, x, t) = V (z, t), if z = w/x < ¯a(t) and ψ(w, x, t) = 0 otherwise.
3.5
The Method of Policy Iteration
Now that a method to solve our problem using the PSOR algorithm has been obtained, we will attempt at solving it with an alternative as well: the Policy Iteration method. Before applying it, we will demonstrate the idea behind the algorithm as presented by Wang and Forsyth (2006) in [14]. Consider the following problem:
Vτ = max
16 Anton Chtepine — Solving Optimal Annuitization and Investment Problems
After discretization and rearrangement of the terms we can rewrite it in the following form for 1 < i < N : Viτ +1− Vτ i ∆τ = maxqiτ +1 Aτ +1 i (qτ +1i )Vi−1τ +1− Biτ +1(qτ +1i )Viτ +1+ Ciτ +1(qiτ +1)Vi+1τ +1+ Diτ +1(qiτ +1) . (3.5.2)
For further work it is convenient to rewrite the above in terms of vectors and matrices. Let Fτ +1 be the following tridiagonal matrix:
Fτ +1 = −B2τ +1 C2τ +1 0 · · · 0 Aτ +13 −B3τ +1 C3τ +1 ... 0 Aτ +14 . .. . .. 0 .. . . .. −BN −2τ +1 CN −2τ +1 0 · · · 0 Aτ +1N −1 −BN −1τ +1 . (3.5.3)
Denote by Dτ +1(qτ +1) a diagonal matrix with elements Dτ +1ii = Diτ +1(qτ +1i ). Let I be the identity matrix. Then 3.5.2 can be modified into
Vτ +1− Vτ
∆τ = maxqτ +1[F
τ +1(qτ +1)Vτ +1+ Dτ +1(qτ +1)]. (3.5.4)
Now, if we find the maximizing vector q∗τ +1 and denote by Fτ +1 and Dτ +1 the optimized values Fτ +1(qτ +1∗ ) and Dτ +1(q∗τ +1), respectively, the equation turns into
(I − ∆τFτ +1)Vτ +1 = Vτ + ∆τ Dτ +1, (3.5.5) where (qτ +1∗ )i = arg maxqτ +1 i [F τ +1(qτ +1)Vτ +1+ Dτ +1(qτ +1)] i. (3.5.6)
The policy iteration algorithm of reaching time step τ +1 from τ can be sketched with the following steps:
1. Define a new iterative process ˆVk,
2. Guess the initial value, for example, ˆV0= Vτ, 3. For k = 0, 1, 2, . . . until convergence solve
(I − ∆τFk) ˆVk+1= Vτ+ ∆τ Dk, (3.5.7) given boundary values V1τ, VNτ and V1τ +1, VNτ +1, with Fk= Fk(qk), Dk= Dk(qk),
4. Obtain the solution at time τ + 1: Vτ +1 = ˆVk+1.
By convergence we mean that the elements of ˆVkseize to change significantly. In other words, we stop the process when
max i
| ˆVik+1− ˆVk i |
max(scale, | ˆVik+1|) < ǫ, (3.5.8) where ǫ is the desired tolerance. The parameter scale is used to insure that unrealistic levels of accuracy are not required, but in most cases scale = 1.
Wang and Forsyth state that in order for the iteration scheme 3.5.7 to converge, the matrix [I − ∆τFk] has to be an M-matrix, i.e. a matrix whose non-diagonal elements are less than or equal to zero and whose eigenvalues’ real parts are positive. Further in the text we will show that this requirement is satisfied in our problem.
3.6
Discretizing the Initial Form of the Problem
Like before, we will find the solution to our problem going backwards in time with the variable change τ = T − t, τ ∈ [0, T ]. This way Vτ = −Vt.
Wang and Forsyth emphasize the importance of selecting the differencing scheme appropriately, so the positive coefficient condition (PCC) is satisfied. What it means is that if we have a term such as maxπ in3.0.4, and it is approximated like
max π (µ − r)πVz+ 1 2σ 2π2V zz i = max π αiV t+1 i−1 + βiVi+1t+1− (αi+ βi+ νi)Vit+1 , (3.6.1) where N− < i < N+, the coefficients α
i, βi, νi have to satisfy
αi≥ 0, βi ≥ 0, νi≥ 0. (3.6.2)
Keeping this in mind, we can view our problem 3.0.4 in discrete form for N−< i < N+: Vik+1− Vik ∆τ + λ S T −τVik+1= rzi− 1 ∆z (V k+1 i − Vi−1k+1) + min πk+1i " (µ − r)πk+1i ∆z (V k+1 i − Vi−1k+1) + σ2(πk+1i )2 2(∆z)2 (V k+1 i+1 − 2Vik+1+ Vi−1k+1) # . (3.6.3) Let πk+1 = (πk+1 N−+1, ..., π k+1
N+−1) be the vector of optimizing values of πk+1 for the above equations. If
we assume the values πk+1i to be known the terms in the approximated equation can be rearranged into Aik+1Vi−1k+1+ Bik+1Vik+1+ Cik+1Vi+1k+1 = 1 ∆τV k i , (3.6.4) in which Ak+1i = (µ − r)π k+1 i + rzi− 1 ∆z − σ2(πk+1i )2 2(∆z)2 , Bik+1 = 1 ∆τ − (µ − r)πk+1i + rzi− 1 ∆z + σ2(πk+1 i )2 (∆z)2 + λ S T −τ, Cik+1 = −σ 2(πk+1 i )2 2(∆z)2 . (3.6.5)
Like in the previous section, 3.6.4 can be rewritten in matrix form: Mk+1Vk+1= 1 ∆τV k − Gk+1, (3.6.6) where Gk+1 = (Ak+1N−+1V k+1 N− , 0, . . . , 0, C k+1
N+−1VNk+1+ )T and the matrix Mk+1 is of the form
Mk+1= BNk+1−+1 C k+1 N−+1 0 · · · 0 Ak+1N−+2 B k+1 N−+2 C τ +1 N−+2 ... 0 Ak+1N−+3 . .. . .. 0 .. . . .. BNk+1+−2 CNk+1+−2 0 · · · 0 Ak+1N+−1 BNk+1+−1 . (3.6.7)
It is easy to verify, that M is a Z-matrix: Ak+1i and Cik+1 are the only off-diagonal elements and they are negative. The numerical experiments show that the real parts of M ’s eigenvalues are always positive, therefore, M is likely to be an M-matrix. M is also diagonally dominant, since |Ak+1i | + |Cik+1| < |Bk+1i |. This fact guarantees the algorithm’s stability.
18 Anton Chtepine — Solving Optimal Annuitization and Investment Problems
3.7
Solution Using Policy Iteration
Now we can apply the previously explained Policy Iteration algorithm to our problem. The most important detail of solving our problem using Policy Iteration is formulating the initial condition. We assume that if an individual has reached a significant age Ξ, his probability of ruin equals zero independent of his wealth. Thus, for the experiments we choose Ξ = 130, and τ changes from 0 to Ξ − T , where T is the attained age, so t starts at Ξ and brings us to T . As specified by the method, to get from time k to k + 1 we choose ˆV0= Vk and solve
MnVˆn+1 = 1 ∆τV
k+ Gn, (3.7.1)
for n = 0, 1, 2, . . . until convergence using the obtained minimizing values πnon each iteration. When the values ˆVn seize to change significantly, we set the solution at time k + 1 equal to ˆVn+1.
In the next section numerical experiments will be presented and, as it turns out, the results of using both methods are identical. Solving the problem with the help of Policy Iteration proved to be computationally faster. Also, it did not require any transforms and changes of variables like in PSOR, but required slightly more programming, which some might consider an advantage and some a drawback.
3.8
Numerical Experiments
In this section we will present the results of several numerical experiments. All of the tests have been performed using both algorithms with the same parameters. These examples are aimed at discovering the effect of small changes of parameters on the problem’s solution.
For the experiments we will use the parameters of the base scenario from Milevsky (2006) [7], if not stated otherwise:
• Gompertz hazard rate λO
t = 1bexp t− ¯bm, where ¯m is a modal value and b is a scale parameter. We choose ¯m = 90 and b = 9. The objective and subjective hazard rates are related by λSt = λOt + η, where η quantifies the individual’s mortality relative to the pricing mortality. We begin with η = 0.
• T = 50; the investor is 50 years old.
• r = 0.02; the risk-free rate is 2% over inflation.
• µ = 0.06; the drift on the risky asset is 6% over inflation. • σ = 0.2; the volatility is 20%.
• c = 1; the individual consumes 1 unit of wealth per unit of time.
All the calculations have been conducted in Matlab. As the wealth of an individual lies between 0 and ¯at, by choosing the number of nodes N = 350 and the size of a time step ∆t = 0.01 we can keep the ratio (∆w)∆t2 ≈ 1.
Note that in case of constant hazard rates Milevsky et al. in [7] obtained an analytical solution for the problem. This solution was lated supported by numerical experiments using the Projected SOR algorithm, therefore, we will compare the results obtained with the Policy Iteration methods with the PSOR results, which are more accurate.
Example 3.8.1. In figures3.1and 3.2we see the ruin probabilities ψ(w, 1, t) of individuals aged 30, 50 and 70 for the base scenario obtained using PSOR and Policy Iteration algorithms. Note, that in both cases the ruin probability declines with age and wealth, i.e. of the two investors of the same age, the older one is less likely to ruin; the less funds an individual has, the more likely he is to ruin. It is an intuitively pleasing result.
The resulting probabilities of lifetime ruin obtained by both methods are alike, although, as we see here and later will see from the numerical experiments in the next section, the Policy Iteration
0 5 10 15 20 25 30 35 0 0.5 1 Probability of ruin Wealth T=30 T=50 T=70
Figure 3.1: Projected SOR
0 5 10 15 20 25 30 35 0 0.5 1 Wealth Probability of ruin T=30 T=50 T=70
Figure 3.2: Policy Iteration
algorithm gives a slightly different curve. The most likely reason for this is the initial condition in case of Policy Iteration. Our assumption in this respect may have affected the solution, while the first algorithm did not involve such an approximation.
Example 3.8.2. This experiment deals with the impact of individual-specific mortality on the ruin probability of an individual age 65 (to avoid negative force of mortality). We use the Gompertz hazard rate as the pricing mortality λOt to define the individual-specific mortality by λSt = λOt + η. The parameter η takes 3 values for this experiment: -0.005, 0 and 0.005. This result does not require a visual representation because we do not observe any significant changes in the probability of ruin. This can be explained by a sort of “adaptation” to the mortality, i.e. change of investment strategy to compensate for the change in subjective mortality. However, in the next, similar experiment the difference in ruin probabilities is more pronounced.
Example 3.8.3. This is basically the same experiment as the previous one, but instead of Gompertz mortality rates, we take λOt = 0.04. We consider an individual whose subjective mortality is given by λS
t = λOt + η for η = -0.015, 0 and 0.015. In this case the difference in ruin probabilities is more noticeable: the higher subjective mortality, the lower ruin probability is, which is intuitively correct, since an individual will live less. This is consistent with the first experiment. The experiment is illustrated by figures3.3and 3.4.
Example 3.8.4. For the next experiment we restrict our attention to the individual aged 50. Our goal now is to see how the probability of ruin depends on the risky asset’s volatility. Specifically, we consider σ = 0.1, 0.2, and 0.5.
Both figures 3.5 and 3.6 show an increase in ruin probability with an increase in volatility for a fixed value of wealth. This result is similar to the one presented by Bayraktar in [1]. Again, notice the slightly larger ruin probabilities in case of the Projected SOR.
Example 3.8.5. Note that the term (µ−r)/σ2 in our problem3.0.4can be regarded as the equity risk premium, i.e. the excess return of a risky asset over the risk-free rate. The next numerical experiment researches the effect of varying this relation on the ruin probability. Let ER = (µ − r)/σ2. Figures3.7
20 Anton Chtepine — Solving Optimal Annuitization and Investment Problems 0 5 10 15 0 0.5 1 Probability of ruin Wealth Eta=−0.015 Eta=0 Eta=0.015
Figure 3.3: Projected SOR
0 5 10 15 0 0.5 1 Wealth Probability of ruin Eta=−0.015 Eta=0 Eta=0.015
Figure 3.4: Policy Iteration
0 5 10 15 20 25 0 0.5 1 Probability of ruin Wealth Sigma=0.1 Sigma=0.2 Sigma=0.5
Figure 3.5: Projected SOR
0 5 10 15 20 25 0 0.5 1 Wealth Probability of ruin Sigma=0.1 Sigma=0.2 Sigma=0.5
Figure 3.6: Policy Iteration
and 3.8illustrate the ruin probability for the initial set of parameters, where ER = 0.16, and for the values ER = 0.40 and ER = 0.64. For instance, such values of ER can be obtained with the values of µ = 0.06, µ = 0.12 and µ = 0.18.
We notice that as ER increases, the ruin probability decreases. The result is intuitively pleasing because of the following: ER can only be increased when either µ is increased, or one or both of r and σ is lowered. If µ is increased, the return on the risky asset gain in value, so an individual is less likely to ruin; if r is decreased, the retiree’s purchasing power will be maintained on a higher level, so
0 5 10 15 20 25 0 0.5 1 Probability of ruin Wealth ER=0.16 ER=0.40 ER=0.64
Figure 3.7: Projected SOR
0 5 10 15 20 25 0 0.5 1 Wealth Probability of ruin ER=0.16 ER=0.40 ER=0.64
Figure 3.8: Policy Iteration
he is again less likely to ruin; and, finally, the effects of varying volatility are presented in one of the previous experiments and are consistent with the current example.
In the figures presented in this section one can observe a slight difference in the solutions obtained by using two different methods to solve one problem. The Projected SOR algorithm may be regarded as a method somewhat “tailored” to the task, while the Policy Iteration was initially intended as an algorithm for solving a wider range of problems. Moreover, the PSOR algorithm leads to accurate results in case of constant hazard rates, as it was compared with the exact results in this case by Milevsky et al. in [7]. Since the Policy Iteration algorithm depends on initial and boundary conditions to a larger extent than PSOR, it is only logical that a minor difference is present, because of the assumption that the probability of ruin at a sufficiently high age is zero.
Chapter 4
The Optimal Annuitization Problem
This chapter deals with the second problem stated in chapter2both in restricted and unrestricted markets. For both markets we restrict our attention to CRRA utility functions like 2.2.1. In case of a restricted market, our goal is to obtain a solution to the following variational inequality:
(λSx+t+ r)V ≥ Vt+ rwVw+ max π (µ − r)πVw+ 1 2σ 2π2V ww + max c [u(c) − cVw], (4.0.1) V ¯ aS x+t ≥ u w ¯ aO x+t (4.0.2) with equation in at least one of them.
If the market is unrestricted, as explained in chapter2, we can limit our attention to the function V (z, t) = U (z, 1, t), z = w/A, instead of the original value function U (w, A, t). Therefore, in the second case we must solve the following variational inequality (with equation in at least one of them):
(λSx+t+ r)V ≥ Vt+ (rz + 1)Vz+ max ˆ π (µ − r)ˆπVz+ 1 2σ 2ˆπ2V zz + max ˆ c≥0[u1(ˆc) − ˆcVz] − kλ S x+tu1(z), (4.0.3) (z + ¯aOx+t)Vz ≥ (1 − γ)V, (4.0.4)
in which ˆπ = π/A and ˆc = c/A.
Both π and c have to be appropriately bounded. We assume the lower bound for both π and c to be zero. The upper bounds are selected to be positive high enough values cmax and πmax.
In all calculations λO
x+t= 1be(x+t− ¯m)/b (Gompertz hazard rate) and λSx+t is specified separately.
4.1
Restricted Market
4.1.1 The Optimal Time of Annuitization
To solve the variational inequality 4.0.1, 4.0.2 we assume that there exists a time T when an individual annuitizes his wealth. This way, for any time t ∈ (0, T )4.0.1 holds as an equality and4.0.2 holds as an inequality. At t = T 4.0.2 turns into an equality. Therefore, it is of great importance to find that time T , when it becomes optimal for the individual to annuitize.
Consider the point t = T where4.0.2 is an equality: V = ¯aSx+tu w ¯ aO x+t = a¯ S x+t (¯aO x+t)1−γ w1−γ 1 − γ. (4.1.1)
To find T we substitute the above into4.0.1 and solve it for t. To do that, the maxc and maxπ terms need to be computed. This can be performed explicitly by taking the derivatives w.r.t. c and π, so the maximizing c and π are found to be
c∗= Vw−1/γ = w ¯ ax+t , π∗ = −Vw(µ − r) Vwwσ2 = w(µ − r) σ2γ , (4.1.2) 22
and max π (µ − r)πVw+ 1 2σ 2π2V ww = −V 2 w(µ − r)2 2Vwwσ2 , max c [u(c) − cVw] = − γ 1 − γ(Vw) (γ−1) γ . (4.1.3)
If we assume relationship 4.0.2 to be an equality which is differentiable, the derivatives of 4.1.1 can calculated to be Vt= V (¯aS x+t) ′ t ¯ aS x+t − (1 − γ)(¯a O x+t) ′ t ¯ aO x+t ! , Vw = V1 − γ w , Vww= −V γ(1 − γ) w2 , (4.1.4)
and, after substituting4.1.3 and 4.1.4 into4.0.1, the equation for finding T has the form (¯aSx+t)′t ¯ aS x+t − (1 − γ)(¯a O x+t) ′ t ¯ aO x+t + γ ¯ aO x+t ¯ aSx+t ¯ aO x+t −1/γ + (1 − γ)(µ − r) 2 2γσ2 − rγ − λ S x+t= 0. (4.1.5) Solving it analytically doesn’t seem possible, but a numerical solution is fairly easy. Like in [7], for the numerical experiments we choose the so-called “proportional hazard” transformation ¯aSx+t= (1 + f )¯aO
x+t, where f ranges from −1 (immortal) to infinity (death at any moment now). It can be shown that (¯aOx+t)′t= λOx+t ¯ aOx+t− Z ∞ 0
eu(−r+1/b)+bλOx+t(1−eu/b)du
, (4.1.6) (¯aSx+t)′t= λSx+t ¯ aSx+t− Z ∞ 0
eu(−r+1/b)+bλSx+t(1−eu/b)du
. (4.1.7)
Due to bλS
x+t> 0, bλOx+t> 0 and the term −eu/b, the integrals in the expressions above converge.
4.1.2 Optimal Consumption and Investment Strategies
Consider the exact solution to our problem, given by ¯U (w, t; T ) = 1/(1 − γ)w1−γφγ(t; T ) in the previous chapter, with φ(t; T ) represented by2.2.6. In this subsection we will prove that the optimal consumption and investment strategies c and π are linear in w. First of all, the first and second derivatives of U are
Uw = w−γφγ(t; T ) = 1 − γ
w U, Uww= −
γ(1 − γ)
w2 U. (4.1.8)
Now, using4.1.2, we get
c∗= w
φ(t; T ), π
∗ = wµ − r
γσ2 . (4.1.9)
Therefore, c and π are linear in w, indeed.
4.1.3 Finite-Difference Approximation
Similarly to the previous chapter we consider a regular mesh with steps ∆w and ∆t such that 0 ≤ i∆t ≤ T and wmin ≤ i∆w ≤ wmax, where T is the optimal annuitization time and wmin and wmax are the minimum and maximum values of an individual’s wealth corresponding to nodes 1 and N , respectively.
As it usually is in finance, we solve our problem backwards in time from T to 0 via change of variables to τ = T − t, so Vτ = −Vt. Therefore, we can rewrite 4.0.1 in the form
Vτ = −(λSx+t+ r)V + rwVw+ max π (µ − r)πVw+ 1 2σ 2π2V ww + max c [u(c) − cVw]. (4.1.10)
24 Anton Chtepine — Solving Optimal Annuitization and Investment Problems
Like in the previous chapter, we have to select the differencing scheme appropriately, so the positive coefficient condition (PCC) is satisfied. If we have terms such as maxπ and maxc in4.1.10, and they are approximated by max π (µ − r)πVw+ 1 2σ 2π2V ww i = max π αiV t+1 i−1 + βiVi+1t+1− (αi+ βi+ νi)Vit+1 , max c [u(c) − cVw]i = maxc [ωiV t+1 i−1 + ζiVi+1t+1− (ωi+ ζi+ φi)Vit+1], (4.1.11) where 1 < i < N , the coefficients αi, βi, νi, ωi, ζi, φi have to satisfy
αi≥ 0, βi ≥ 0, νi≥ 0, ωi ≥ 0, ζi≥ 0, φi ≥ 0. (4.1.12) As we have shown in4.1.3, explicitly maximizing the terms in4.0.1leads to a differential equation with derivatives in inconvenient powers. Instead of doing so, we will show that for solving the problem numerically it is much more convenient to discretize first and optimize afterwards. Thus, as a result of maximizing 4.1.11, we have a discrete optimal controls ci and πi for 1 < i < N . Note, that the coefficients αi, βi, νi, ωi, ζi, φi are, therefore, all functions of the corresponding control:
αi= αi(πi), βi= βi(πi), νi = νi(πi), ωi = ωi(ci), ζi = ζi(ci), φi = φi(ci). (4.1.13) We can now write down our problem and mark the appropriate differencing scheme:
Vτf wd = −(λSx+t+ r)V + rwVwf wd+ max c [u(c) − cV bwd w ] + maxπ (µ − r)πVwf wd+ 1 2σ 2π2V ww , (4.1.14) where superscript f wd stands for forward differencing. This choice satisfies PCC when parameters µ > 0 and r > 0 are selected in a way that µ > r. The necessary differencing schemes are given by
Vτf wd = V τ +1 i − Viτ ∆τ , (4.1.15) Vwf wd = V τ +1 i+1 − Viτ +1 ∆w , (4.1.16) Vwbwd = V τ +1 i − Vi−1τ +1 ∆w , (4.1.17) Vww= Vi+1τ +1− 2Viτ +1+ Vi−1τ +1 (∆w)2 . (4.1.18)
Taking this and4.1.14 into account, we obtain the discrete form of our problem: Viτ +1− Vτ i ∆τ = − (λ S x+t+ r)Vτ +1+ rwi ∆w(V τ +1 i+1 − Viτ +1) + max ci h ci ∆w(V τ +1 i−1 − Viτ +1) + u(ci) i + max πi σ2π2i 2(∆w)2(V τ +1 i+1 − 2Viτ +1+ Vi−1τ +1) +(µ − r)πi ∆w (V τ +1 i+1 − Viτ +1) , (4.1.19) Viτ +1 ≥¯aSx+τ +1u wi ¯ aO x+τ +1 ! = ¯a S x+τ +1 (¯aS x+τ +1)1−γ w1−γi 1 − γ. (4.1.20)
4.1.4 Solution via Policy Iteration in Matrix Form
Now we have the tool for solving 4.1.19. Let cτ +1 = (cτ +1
2 , ..., cτ +1N −1) and πτ +1 = (π2τ +1, ..., πτ +1N −1) be the maximizing vectors. We can now omit the max operators, and after rearranging the terms in our equation we get
Aiτ +1Vi−1τ +1+ Biτ +1Viτ +1+ Ciτ +1Vi+1τ +1 = 1 ∆τV
τ
where Aτ +1i = −c τ +1 i ∆w − σ2(πiτ +1)2 2(∆w)2 , (4.1.22) Biτ +1= 1 ∆τ + λ S x+t+ r + rwi ∆w + cτ +1i ∆w + σ2(πiτ +1)2 (∆w)2 + (µ − r)πτ +1i ∆w , (4.1.23) Ciτ +1= −(µ − r)π τ +1 i ∆w − σ2(πτ +1i )2 2(∆w)2 − rwi ∆w. (4.1.24)
Similarly to [14] it is convenient to rewrite 4.1.21 in matrix form. Denote by Mτ +1 the matrix
Mτ +1 = B2τ +1 C2τ +1 0 · · · 0 Aτ +13 B3τ +1 C3τ +1 ... 0 Aτ +14 . .. . .. 0 .. . . .. BN −2τ +1 CN −2τ +1 0 · · · 0 Aτ +1N −1 Bτ +1N −1 . (4.1.25)
Then 4.1.21 turns into
Mτ +1Vτ +1= 1 ∆τV τ+ Gτ +1, (4.1.26) in which Gτ +1 = u(cτ +12 ) u(cτ +13 ) .. . u(cτ +1N −2) u(cτ +1N −1) − Aτ +12 V1τ +1 0 .. . 0 CN −1τ +1VNτ +1 . (4.1.27)
Due to unspecified boundary values VNτ and V1τ +1 in all numerical experiments we will presume that on the boundaries the second-order derivatives of V are equal to zero, i.e. V1 = 2V2− V3 and VN = 2VN −1− VN −2. There have been problems with this kind of boundary conditions but this will be discussed in more detail further in the text.
Now that we have constructed the matrix M , it is trivial to prove that all non-diagonal elements are negative, while all the elements on the main diagonal are positive. Moreover, we can see that M is diagonally dominant, which reassures us that the algorithm is computationally stable. For M to be an M-matrix its eigenvalues’ real parts have to be positive. While this paper does not formally prove this fact, the numerical experiments show that it is, in fact, true.
To proceed from time τ to τ + 1 in accordance with the policy iteration method we guess the initial value ˆV0 to be equal to the solution at time τ , i.e. ˆV0 = Vτ. Then, until convergence, we solve 4.1.26:
MkVˆk+1= 1 ∆τV
τ+ Gk, (4.1.28)
using the optimal values of ck and πk obtained on each iteration of the process. This can easily be solved by, for instance, the tridiagonal algorithm, or simply by inverting M . The iterations stop when the desired accuracy is reached and the final ˆVk+1 is set to be the solution at time τ + 1.
There are two options for enforcing the constraint4.0.2. The inequality itself means that at each node i we take the maximum of the obtained solution Viτ +1 and the constraint at that time, i.e. for each i: Viτ +1 = max Viτ +1, ¯aSx+T −τ −1u wi ¯ aO x+T −τ −1 !! . (4.1.29)
26 Anton Chtepine — Solving Optimal Annuitization and Investment Problems
This operation can be performed either after the policy iteration algorithm has delivered us to the solution at time τ + 1 or at every iteration of the algorithm itself. The numerical experiments will show that these options are absolutely equal, and no differences have been observed.
The explained algorithm for solving4.0.1,4.0.2may give rise to concerns about its lack of something Wilmott in [15] calls “internal consistency”. It means that when enforcing the constraint a portion of the solutions might get simply cut off and substituted by the constraint itself. Solving our problem with a hybrid algorithm of Policy Iteration and the Projected SOR may prove to be a better way to obtain a result.
4.1.5 Policy Iteration - PSOR Hybrid Algorithm
To ensure that “internal consistency” is present we propose to combine the explained policy itera-tion algorithm and the PSOR method stated in the previous chapter. The idea is to get from time τ to τ + 1 using an iterative process like before, but instead of solving4.1.28 in matrix form and enforcing the constraint on the solution afterwards, successively construct the values ˆVik+1 using PSOR method and enforce the constraint on each value. This way the constraint on one element immediately affects the calculation of the rest.
Keeping in mind the form of the matrix Mk(4.1.25) and vector Gk(4.1.27), the proposed algorithm can be divided into the following steps:
1. Define a new iterative process ˆVk,
2. Guess the initial value, for example, ˆV0= Vτ, 3. For k = 0, 1, 2, . . . until convergence solve
(Mk) ˆVk+1= 1 ∆τV
τ+ Gk. (4.1.30)
using PSOR algorithm: for every 1 < i < N do yk+1i = ( 1 ∆τV τ i + Uik− AkiVˆi−1k+1− CikVˆi+1k+1)/Bik), Vik+1 = max Vik+ ω(yik+1− Vik), ¯aSx+T −τu wi ¯ aO x+T −τ !! , (4.1.31)
where Uik= (u(ck2), . . . , u(ckN −1))T, and 1 < ω < 2 is the over-relaxation factor, 4. Obtain the solution at time τ + 1: Vτ +1 = ˆVk+1.
The proposed algorithm brings us to the same results as the previous one, although, as the nu-merical experiments will show, this method leads to strong fluctuations of c(w) and π(w), while plain policy iteration keeps these functions relatively smooth. PSOR is also computationally slower due to successive construction of values for each node, whereas matrix operations can be parallelized, therefore, slightly boosting performance.
4.1.6 Numerical Experiments
This subsection presents several numerical examples regarding the optimal time of annuitization and the construction of the value function of the “all-or-nothing” problem, i.e. the restricted market case. We will use the same parameters as in [7]: Gompertz force of mortality λt = 1bexp t−mb with the pair of parameters (m, b) for males and females taking values (88.18, 10.5) and (92.63, 8.78), re-spectively; the risky stock is assumed to have drift µ = 0.12 and volatility σ = 0.2; the risk-free rate is r = 0.06.
Example 4.1.1. In this experiment we will determine the optimal age for annuitization using4.1.5 for three levels of risk-aversion: γ = 1 (logarithmic utility), γ = 2 and γ = 5. Within the brackets next to our results, are the results obtained by Milevsky et al. (2006) [7]. As we will see, the results are exactly the same.
Age Optimal age
γ = 1 γ = 2 γ = 5
Male Female Male Female Male Female
60 80.3 (80.3) 84.5 (84.5) 73.0 (73.0) 78.4 (78.4) 63.4 (63.4) 70.4 (70.4) 65 80.3 (80.3) 84.5 (84.5) 73.0 (73.0) 78.4 (78.4) Now (Now) 70.4 (70.4) 70 80.3 (80.3) 84.5 (84.5) 73.0 (73.0) 78.4 (78.4) Now (Now) 70.4 (70.4) 75 80.3 (80.3) 84.5 (84.5) Now (Now) 78.4 (78.4) Now (Now) Now (Now) 80 80.3 (80.3) 84.5 (84.5) Now (Now) Now (Now) Now (Now) Now (Now) 85 Now (Now) Now (Now) Now (Now) Now (Now) Now (Now) Now (Now) Note, that due to lower mortality rates females annuitize later than males. Also, with the increase in risk-aversion both genders tend to annuitize sooner, but even with such high γ = 5, females prefer not to annuitize before 70 and males before 63. Nevertheless, the result are logical and intuitively pleasing.
Example 4.1.2. Now we will research the dependence of the optimal annuitization age on the sub-jective, individual-specific hazard rate. Let λSx+t= (1 + f )λOx+t, where f ≥ −1. The larger f is, the less “healthy” the retiree is. For this experiment we will consider a male aged 60 with CRRA of γ = 2.
f -1 -0.8 -0.6 -0.4 -0.2 0 0.5 1.0 1.5 2.0
Our results 78.28 74.584 73.709 73.286 73.03 73.078 73.305 74.036 75.211 76.963 Results from [7] 78.28 74.58 73.71 73.29 73.03 73.08 73.31 74.04 75.21 76.96
Table 4.1: Optimal annuitization ages
−1 −0.5 0 0.5 1 1.5 2 2.5 3 70 75 80 85 90 Optimal age f
Figure 4.1: Optimal annuitization ages
As we can see, the numbers are nearly identical. Also, it is clear from figure4.1, that deviation of the subjective hazard rate from the pricing hazard rate increases the optimal annuitization age, no matter in which direction the deviation occurs. See [7] for the proof of this fact in case of ¯aS
x < 2¯aOx. In the following examples we will perform the analysis of the value function’s sensitivity to changes of parameters. In every experiment we consider a male aged 60 with CRRA of γ = 2, if not stated otherwise. In this case the optimal age for annuitization is 73. Note, that for γ > 1 the utility function contains division by zero at w = 0. Therefore, in order to avoid “infinities”, we set a small enough number ǫ2= 0.00001 as the lower boundary for an individual’s wealth. The upper boundary is assumed to be equal to 50 units of wealth. For the numerical solution we use a regular mesh with 400 nodes, so ∆w ≈ 0.125. The time step is ∆t = 0.02. This make the ratio ∆w∆t2 ≈ 1. As mentioned before, we choose the boundary conditions at w = ǫ2 and w = 50 such, that the second-order derivatives of the value function V are zero, i.e. Vt
1 = 2V2t− V3t and VNt = 2VN −1t − VN −2t for all times t.
Note, that from the previous sections and Chapter 2 we know the exact solution to problem 4.0.1,4.0.2 U (w, t; T ) = 1/(1 − γ)w¯ 1−γφγ(t; T ) with φ(t; T ) given by2.2.6. We also know the optimal investment strategies given by4.1.9.
28 Anton Chtepine — Solving Optimal Annuitization and Investment Problems 5 10 15 20 25 30 35 40 45 50 −150 −100 −50 0 Value function Wealth Matrix PI PSOR PI Exact
Figure 4.2: Value functions
0 10 20 30 40 50 0 50 100 150 Optimal strategies Wealth Matrix PI Pi PSOR PI Pi Exact Pi
Figure 4.3: Optimal investments in risky asset
0 10 20 30 40 50 0 2 4 6 8 10 Optimal strategies Wealth Matrix PI C PSOR PI C Exact C
Figure 4.4: Optimal investments in riskless asset
Example 4.1.3. We will begin by demonstrating the differences between the exact value function and investment strategies and the ones obtained by the matrix form of the policy Iteration (Matrix PI) and PSOR and Policy Iteration hybrid algorithm (PSOR PI).
As we can see in figure4.2, there is a slight difference between the two numerical solutions and the exact solution. Due to the fact that our numerical results are almost identical, their deviation from the exact solution may be explained by the boundary conditions which, apparently, do not fully suit the problem.
The results in figures 4.3 and 4.4 show sharp fluctuations of the investment strategies obtained using PSOR PI hybrid algorithm, which is in contradiction with the fact that both π and c are smooth in w. However, the strategies in case of Matrix PI appear to be smooth everywhere but the boundaries. It is worth mentioning, that Matrix PI starts off with the correct π and c on the first iteration, but the boundary conditions force the maximizing π to take very large values. This can be explained by the second-order derivative being zero while the first-order derivative is positive, therefore, the larger π is, the larger maxπ term is. This observation testifies in favor of the boundary conditions being slightly inadequate. However, despite the significant differences between the two numerically obtained investment strategies, the resulting value functions are nearly identical.
In the following experiments we will restrict our attention only to the value function and investment strategies obtained using Matrix PI, because of their similarity to the ones in case of PSOR PI hybrid algorithm.
Example 4.1.4. As mentioned in one of the examples from the previous chapter, the term ER = (µ − r)/σ2 can be interpreted as the market value of equity risk. This numerical experiment will illustrate the dependency of the value function and the optimal investment strategies on ER. The value of ER with all the parameters defined like in the base scenario is 1.5. We will consider values of ER = 1 and ER = 2. For example, this can be achieved by setting µ = 0.10 and µ = 0.14, respectively. The results are presented by figures4.5-4.9.
5 10 15 20 25 30 35 40 45 50 −150 −100 −50 0 Value function Wealth ER=1.0 ER=1.5 ER=2.0
Figure 4.5: Value functions (numerically obtained)
Figure 4.5 demonstrates the numerically obtained value functions for each value of ER. Notice, that the value function increases as the market value of equity risk increases. The growth of ER is a result of either the excess return on the risky asset being larger than the risk-free rate, or the volatility decrease. Both cases justify the results. While the differences between the numerically obtained value functions are slightly noticeable, the exact results are much less pronounced but similar. The exact result, therefore, do not require visualization.
As the excess return on the risky asset increases, the logical thing to do is to invest more in it. This strategy is supported by figure4.6, where the exact amounts invested in the risky asset are illustrated. Notice, that the numerical results (figure4.7) deviate from the exact strategies, which can be explained by the boundary conditions similarly to the previous example. As we see, larger ER leads to a less pronounced bend in the graph of π, so, based on the numerical results, the wealth-axis can be split into two segments wa and wb, where for a pair of values of ER the amount π is larger for the first value of ER in wa, and the other way around in wb.
30 Anton Chtepine — Solving Optimal Annuitization and Investment Problems 0 10 20 30 40 50 0 50 100 150 Optimal strategies Wealth ER=1.0 ER=1.5 ER=2.0
Figure 4.6: Optimal investments in risky asset (exact)
0 10 20 30 40 50 0 50 100 150 Optimal strategies Wealth ER=1.0 ER=1.5 ER=2.0
Figure 4.7: Optimal investments in risky asset (numerically obtained)
0 10 20 30 40 50 0 2 4 6 Optimal strategies Wealth ER=1.0 ER=1.5 ER=2.0
Figure 4.8: Optimal investments in riskless asset (exact)
Figures4.8 and4.9demonstrate the dependency of the investment in a riskless asset on the value of ER. Both, exact and numerical results, decrease with lowering ER. This tendency might be the result of larger returns on the risky asset in case of a larger ER which, in turn, leads to the retiree being able to consume more. As before, there is a bend in the numerical solution, possibly due to a flaw in the boundary conditions.
