Add code at the appropriate location$($s$)($and make all additional necessary changes throughout the program$)$ to implement the $$fixed percentage$$ rule. When necessary, review the Horneff et al$.$ paper to remind yourself what that rule does. The final program is supposed to either simulate the $$fixed benefit$$ or the $$fixed percentage$$ rule, depending on user input provided when calling the evaluateWithdrawalPlan$()$ function.

Once your program produces the desired result, export the plot for your $$fixed percentage$$ simulation as a pdf file $($manually is fine, no need to write corresponding code$).$

parameters $=$ list$($m $=$ c$(0\mathrm{.}1155,0\mathrm{.}0845),$ # Expected annual $($log$-)$return

v $=$ c$(0\mathrm{.}1533,0\mathrm{.}1028),$ # Annual $($log$-)$return standard dev

rho $=0\mathrm{.}33,$ # Correlation between stocks, bonds

stockFraction $=0\mathrm{.}6,$ # Fraction of periodic wealth invested into stocks

deltaT $=1,$ # set to $1,$ because simulation on annual basis

ageInitial $=65,$

ageEnd $=100,$

seed $=123456,$ # Arbitrary "start value" for random number generator

nIterations $=100000,$ # Number of simulation runs

initialWealth $=100000,$

fixedBenefitAmount $=7200)$

generateLogReturns $=$ function$($parameters$)\{$

# Simulation of asset returns

# Extract individual parameters from parameters list $($not technically necessary

# but improves readability of code below$)$

m $=$ parameters$m

v $=$ parameters$v

rho $=$ parameters$rho

deltaT $=$ parameters$deltaT

nIterations $=$ parameters$nIterations

# Generate $($nIterations x $2)$ matrix of N$(0,1)$ shocks

z $=$ rnorm$(2*$nIterations$)$

z $=$ matrix$($z$,$ nrow $=$ nIterations, ncol $=2)$

# Transform N$(0,1)$ shocks into $($nIterations x $2)$ matrix of log return

# realizations for the $2$ assets using the multivariate Geometric Brownian

# Motion approach

returns $=$ matrix$($nrow $=$ nIterations, ncol $=2)$

returns$[,1]=$ m$[1]*$deltaT $+$ v$[1]*$sqrt$($deltaT$)*$z$[,1]$ # Returns for asset $1($stock$)$

returns$[,2]=$ m$[2]*$deltaT $+$ # Returns for asset $2($bond$)$

rho$*$v$[2]*$sqrt$($deltaT$)*$z$[,1]+$

sqrt$((1-$rho$^2))*$v$[2]*$sqrt$($deltaT$)*$z$[,2]$

return$($returns$)$

$\}$

calculateBenefit $=$ function $($withdrawalRule$,$ beginningOfPeriodWealth, t$,$ parameters$)\{$

# Determine benefit based on withdrawal rule and wealth at a particular point in time

# Inputs:

# withdrawalRule: Character string naming the type of withdrawal rule to be analyzed

# $("$fixedBenefit$",$ "fixedPercentage", "oneOverT", "oneOverET"$)$

# beginningOfPeriodWealth: $($nIterations x $1)$ vector holding the wealth at the beginning

# of a period for each of the nIterations simulation runs

# t: Index of current period $($for later use$)$

# In each if structure below, an $($nIterations x $1)$ vector of benefit payments will

# be determined depending on the rules of the particular withdrawal plan design

if $($withdrawalRule $==$ "fixedBenefit"$)\{$

# In this plan design, the benefit is a pre$-$specified constant amount as long

# as sufficient funds are available to pay the benefit, i$.$e$.$ as long as

# $($benefit amount $<$ available wealth$).$ In case funds are no longer sufficient

# to pay the full benefit, the benefit is reduced to the amount of funds remaining.

# pmin$()$ compares each individual entry in the beginningOfPeriodWealth vector with

# the fixedBenefitAmount and finds the smaller of the two values.

benefit $=$ pmin$($beginningOfPeriodWealth$,$ parameters$fixedBenefitAmount$)$

$\}$

if $($withdrawalRule $==$ "fixedPercentage"$)\{$

# t$.$b$.$d$.$

$\}$

if $($withdrawalRule $==$ "oneOverT"$)\{$

# t$.$b$.$d$.$

$\}$

if $($withdrawalRule $==$ "oneOverET"$)\{$

# t$.$b$.$d$.$

$\}$

return $($benefit$)$

$\}$