HARS - Historical Actuarial Retirement Simulator

\n"; echo "Author (bug reports welcome) Source code\n"; echo "

\n"; echo "This is free software. Use it at your own risk. No guarantees of correctness are made.\n"; echo "

\n"; echo "

\n"; printf("initial assets %4dk max assets %4dk return on equity %3.1f%% return on bonds %3.1f%%\n", $total / 1000, $limit / 1000, $ret_equity, $ret_bonds); echo "
\n"; printf("years until retire %2d annual savings until retire %3dk annual living expenses once retire %3dk\n", $years_to_retire, $savings / 1000, $income / 1000); echo "
\n"; printf("equity diversification %d%% correlation %d%% standard deviation scaling factor %3.1f\n", $diversification * 100, $correlation * 100, $diversify_sds); echo "
\n"; printf("go to bonds %5.1f%% income dead fraction %3d%% age %d/%d gender %s/%s\n", $go_to_bonds, $income_dead_fraction * 100, $age1, $age2, $sex1, $sex2); echo "
\n"; printf("simulation order mode '%3s' asset depletion mode '%2s'\n", $seq, $deplete); echo "
\n"; printf("simulation start year %d end year %d experiments per result %d\n", $sim_start_year, $shiller_initial + count($shiller_price) / 12 - 1, $num_experiments); # Construct array of equity values to use $equity_fraction = array(); for ($i = 0; $i <= $equity_steps; $i++) { array_push($equity_fraction, $i / $equity_steps); } if (($age1 != -1) && ($deplete == 'rb')) { array_push($equity_fraction, 'age-in-bonds'); } $diversify_factor = (1 - $correlation) * $diversification; $death1 = (($sex1 == 'male') ? $male : $female); $death2 = (($sex2 == 'male') ? $male : $female); # Compute shiller annual returns. $shiller_real_price = array(); $shiller_real_dividend = array(); for ($i = 0; $i <= count($shiller_price) -1 ; $i++) { array_push($shiller_real_price, $shiller_price[$i] * $shiller_cpi[count($shiller_cpi) - 1] / $shiller_cpi[$i]); array_push($shiller_real_dividend, $shiller_dividend[$i] * $shiller_cpi[count($shiller_cpi) - 1] / $shiller_cpi[$i]); } $index = 1; $shiller = array(); $cpi_level = array(); for ($i = 1; $i <= count($shiller_real_price) - 1; $i++) { $index = $index / $shiller_real_price[$i - 1] * ($shiller_real_price[$i] + $shiller_real_dividend[$i] / 12); if ($i % 12 == 11) { # Just simulate every 12 monthly value for speed. array_push($shiller, $index); array_push($cpi_level, $shiller_cpi[$i]); } } # The terminal value of the index differs from the terminal value computed using a spreadsheet by 0.5%. # This represents 0.0003% per iteration. This difference is assumed to be due to rounding differences. # The difference in successive index values, which is what we care about, is far less. # Compute Shiller returns $shiller_offset = $sim_start_year - ($shiller_initial + 1); for ($i = 0; $i < (count($shiller) - ($shiller_offset + 1)); $i++) { $returns[$i] = ($shiller[$shiller_offset + $i + 1] - $shiller[$shiller_offset + $i]) / $shiller[$shiller_offset + $i]; } # Compute diversify to have a scaled standard deviation $diversify = array(); for ($i = 0; $i <= count($returns) - 1; $i++) { $diversify[$i] = $returns[$i] * $diversify_sds; } # Normalize diversify to have a geometric mean of $ret_equity and a correspondingly scaled standard deviation $samp_prod = 1; for ($i = 0; $i <= count($diversify) - 1; $i++) { $samp_prod *= 1 + $diversify[$i]; } $samp_geomean = pow($samp_prod, 1 / count($diversify)) - 1; $adjust = (1 + $ret_equity / 100) / (1 + $samp_geomean); for ($i = 0; $i <= count($diversify) - 1; $i++) { $diversify[$i] = (1 + $diversify[$i]) * $adjust - 1; } # XXX 2009-09 - Don't understand following 5 lines - comment them out. Had been increasing volatility and bringing down success probabilities. # Compute returns to have a scaled standard deviation #$returns_sds = (1 - $diversification) + $diversify_sds * $diversification; #for ($i = 0; $i <= count($returns) - 1; $i++) { # $returns[$i] = $returns[$i] * $returns_sds; #} # Normalize returns to have a geometric mean of $ret_equity and a correspondingly scaled standard deviation $samp_prod = 1; for ($i = 0; $i <= count($returns) - 1; $i++) { $samp_prod *= 1 + $returns[$i]; } $samp_geomean = pow($samp_prod, 1 / count($returns)) - 1; $adjust = (1 + $ret_equity / 100) / (1 + $samp_geomean); for ($i = 0; $i <= count($returns) - 1; $i++) { $returns[$i] = (1 + $returns[$i]) * $adjust - 1; } # Report statistics $samp_prod = 1; for ($i = 0; $i <= count($returns) - 1; $i++) { $samp_prod *= 1 + $returns[$i]; } $samp_geomean = pow($samp_prod, 1 / count($returns)) - 1; $samp_tot = 0; for ($i = 0; $i <= count($returns) - 1; $i++) { $samp_tot += $returns[$i]; } $samp_mean = $samp_tot / count($returns); $samp_ss = 0; for ($i = 0; $i <= count($returns) - 1; $i++) { $samp_ss += pow($returns[$i] - $samp_mean, 2); } $samp_sd = pow(1 / (count($returns) - 1) * $samp_ss, 0.5); echo "

\n"; echo "\n"; printf("equity return = geometric mean %4.2f%% arithmetic mean %4.2f%% +/- %4.2f%%\n", $samp_geomean * 100, $samp_mean * 100, $samp_sd * 100); $samp_prod = 1; for ($i = 0; $i <= count($diversify) - 1; $i++) { $samp_prod *= 1 + $diversify[$i]; } $samp_geomean = pow($samp_prod, 1 / count($diversify)) - 1; $samp_tot = 0; for ($i = 0; $i <= count($diversify) - 1; $i++) { $samp_tot += $diversify[$i]; } $samp_mean = $samp_tot / count($diversify); $samp_ss = 0; for ($i = 0; $i <= count($diversify) - 1; $i++) { $samp_ss += pow($diversify[$i] - $samp_mean, 2); } $samp_sd = pow(1 / (count($diversify) - 1) * $samp_ss, 0.5); echo "
\n"; printf("diversify return = geometric mean %4.2f%% arithmetic mean %4.2f%% +/- %4.2f%%\n", $samp_geomean * 100, $samp_mean * 100, $samp_sd * 100); # Compute percentile life expectancies in order to be able to implement go to bonds strategy if ($go_to_bonds != -1) { if (($age1 != -1) && ($age2 != -1)) { $death_prob = array(); # $death_prob[$y][$e] is the probability both members of a couple are dead at the end of year $e if they were both alive at the start of year $y assert(count($death1) == count($death2)); for ($y = 0; min($age1, $age2) + $y <= count($death1) - 1; $y++) { $alive1 = 1; $alive2 = 1; $death_prob[$y] = array(); for ($e = $y; min($age1, $age2) + $e <= count($death1) - 1; $e++) { $alive1 *= ($age1 + $e >= count($death1)) ? 0 : 1 - $death1[$age1 + $e]; $alive2 *= ($age2 + $e >= count($death2)) ? 0 : 1 - $death2[$age2 + $e]; $death_prob[$y][$e] = (1 - $alive1) * (1 - $alive2); } } $life_epectancy_both = array(); # $life_expectancy_both[$y] is the $go_to_bonds percentile life expectancy for a couple given they were both alive at the start of year $y # Value should be rounded up for ($y = 0; $y <= count($death_prob) - 1; $y++) { $le = 0; for ($e = count($death1) - 1 - min($age1, $age2); $e >= $y; $e--) { if ($death_prob[$y][$e] <= $go_to_bonds / 100) { $le = $e - $y + 1; break; } } $life_expectancy_both[$y] = $le + 1; } echo "
\n"; printf("%5.1fth percentile additional life expectancy=%d\n", $go_to_bonds, $life_expectancy_both[0]); } if ($age1 != -1) { $death_prob = array(); # $death_prob[$y][$e] is the probability the first person is dead at the end of year $e if they were alive at the start of year $y for ($y = 0; $age1 + $y <= count($death1) - 1; $y++) { $alive1 = 1; $death_prob[$y] = array(); for ($e = $y; $age1 + $e <= count($death1) - 1; $e++) { $alive1 *= 1 - $death1[$age1 + $e]; $death_prob[$y][$e] = 1 - $alive1; } } $life_epectancy_1 = array(); # $life_expectancy_1[$y] is the $go_to_bonds percentile life expectancy for the first person given they were alive at the start of year $y # Value should be rounded up for ($y = 0; $y <= count($death_prob) - 1; $y++) { $le = 0; for ($e = count($death1) - 1 - $age1; $e >= $y; $e--) { if ($death_prob[$y][$e] <= $go_to_bonds / 100) { $le = $e - $y + 1; break; } } $life_expectancy_1[$y] = $le + 1; } if ($age2 == -1) { echo "
\n"; printf("%5.1fth percentile additional life expectancy=%d\n", $go_to_bonds, $life_expectancy_1[0]); } } if ($age2 != -1) { $death_prob = array(); # $death_prob[$y][$e] is the probability the second person is dead at the end of year $e if they were alive at the start of year $y for ($y = 0; $age2 + $y <= count($death2) - 1; $y++) { $alive2 = 1; $death_prob[$y] = array(); for ($e = $y; $age2 + $e <= count($death2) - 1; $e++) { $alive2 *= 1 - $death2[$age2 + $e]; $death_prob[$y][$e] = 1 - $alive2; } } $life_epectancy_2 = array(); # $life_expectancy_2[$y] is the $go_to_bonds percentile life expectancy for the second person given they were alive at the start of year $y # Value should be rounded up for ($y = 0; $y <= count($death_prob) - 1; $y++) { $le = 0; for ($e = count($death2) - 1 - $age2; $e >= $y; $e--) { if ($death_prob[$y][$e] <= $go_to_bonds / 100) { $le = $e - $y + 1; break; } } $life_expectancy_2[$y] = $le + 1; } } } echo "\n"; echo "
\n"; echo "
\n"; echo "
\n"; echo "
Equity Bonds Success probability \n"; echo " Mean gifted assets Mean assets at final death \n"; srand(0); # Simulate portfolios foreach ($equity_fraction as $ef) { if (is_string($ef) && ($ef == 'age-in-bonds')) { $equity = (1 - min(100, ($age2 == -1) ? $age1 : ($age1 + $age2) / 2) / 100) * $total; } else { $equity = $ef * $total; } $cash = $total - $equity; $p_array = array(); $fail = 0; $tot = 0; $charitable_tot = 0; $experiments = $num_experiments; for ($experiment = 0; $experiment < $experiments; $experiment++) { $e = $equity * $pe_adjust; $c = $cash; $in = $income; $charitable = 0; $dead1 = 0; $dead2 = ($age2 == -1); $index1 = rand(0, count($returns) - 1); $index2 = rand(0, count($diversify) - 1); $gone_to_bonds = 0; for ($y = 0; ($y < $years) && !($dead1 && $dead2); $y++) { if ($seq == 'seq') { $e *= (1 - $diversify_factor) * (1 + $returns[$index1]) + $diversify_factor * (1 + $diversify[$index2]); $index1 = ($index1 + 1) % count($returns); $index2 = ($index2 + 1) % count($diversify); } else { $e *= (1 - $diversify_factor) * (1 + $returns[rand(0, count($returns) - 1)]) + $diversify_factor * (1 + $diversify[rand(0, count($diversify) - 1)]); } $c *= 1 + $ret_bonds / 100; $p = $e + $c; if ($y >= $years_to_retire) { if ($dead1 || ($dead2 && ($age2 != -1))) { $p -= $in * $income_dead_fraction; } else { $p -= $in; } } else { $p += $savings; } if (($limit > 0) && ($p > $limit)) { $charitable += $p - $limit; $p = $limit; } if (($go_to_bonds != -1) && ! $gone_to_bonds) { if ($dead1) { $le = $life_expectancy_2[$y]; if ($age1 != -1) { $le *= $income_dead_fraction; } } elseif ($dead2) { $le = $life_expectancy_1[$y]; if ($age2 != -1) { $le *= $income_dead_fraction; } } else { $le = $life_expectancy_both[$y]; } $gone_to_bonds = ($p >= $le * $in); } if ($gone_to_bonds) { $c = $p; } elseif (is_string($ef) && ($ef == 'age-in-bonds')) { $c = min(100, ($age2 == -1) ? $age1 : ($age1 + $age2) / 2) / 100 * $p; } elseif ($deplete == 'rb') { $c = (1 - $ef) * $p; } else { if ($y < $years_to_retire) { # use savings contribution to try and maintain bond fraction up to entire contribution amount $c = min(max($c, (1 - $ef) * $p), $c + $savings); } if ($deplete == 'bf') { if ($y >= $years_to_retire) { if ($dead1 || ($dead2 && ($age2 != -1))) { $c -= $in * $income_dead_fraction; } else { $c -= $in; } } } else { # ($deplete == 'ef'); $c = $c; # leave alone } } if ($c < 0) { $c = 0; } $e = $p - $c; if ($e < 0) { $c = $p; $e = 0; } if (($age1 != -1) && !$dead1) { $dead1 = (rand(0, 9999) / 10000 < $death1[$age1 + $y]); } if (!$dead2) { $dead2 = (rand(0, 9999) / 10000 < $death2[$age2 + $y]); } } array_push($p_array, $p); if ($p < 0) { $fail++; } $tot += $p; $charitable_tot += $charitable; } # Compute and report stats # sort($p_array); # $pctl_0 = $p_array[0]; # $pctl_10 = $p_array[intval((count($p_array) + 1) / 10)]; # $pctl_25 = $p_array[intval((count($p_array) + 1) / 4)]; # $pctl_50 = $p_array[intval((count($p_array) + 1) / 2)]; $mean = $tot / $experiments; $charitable_mean = $charitable_tot / $experiments; # if ($pctl_0 < 0) { # $pctl_0_display = "-$" . round(- $pctl_0 / 1000) . "k"; # } else { # $pctl_0_display = "$" . round($pctl_0 / 1000) . "k"; # } if ($mean < 0) { $mean_display = "-$" . round(- $mean / 1000) . "k"; } else { $mean_display = "$" . round($mean / 1000) . "k"; } if (is_string($ef) && ($ef == 'age-in-bonds')) { echo "
average age in bonds\n"; } else { printf("
%d%% %d%%\n", round($ef * 100), round((1 - $ef) * 100)); } printf(" %d%% $%dk %s\n", round(($experiments - $fail) / $experiments * 100), round($charitable_mean / 1000), $mean_display); } echo "
\n"; echo "

\n"; echo "This retirement simulator simulates an initial set of assets with a savings accumulation phase followed by an income withdrawal retirement phase. The results associated with holding different ratios of stocks/bonds are reported.\n"; echo "

\n"; echo "The reported results show the probability of dying without running out of money (success). For a couple the death of one partner reduces the amount of income required to support the other partner, but simulation continues until both partners are dead.\n"; echo "

\n"; echo "To use this tool effectively, you should perform multiple runs of the simulator changing the input parameters in order to assess the sensitivity of your results to small changes in the inputs.\n"; echo "

\n"; echo "Equity returns are based on the historical returns reported by Shiller, but with the mean return adjusted to match a specified simulation parameter. This is because it seems unlikely that the return values of the past will be repeated, but they are still informative in setting the volatility of future returns.\n"; echo "

\n"; echo "Bond returns are assumed to be a fixed amount per year after adjusting for inflation, and as such most closely model TIPS.\n"; echo "

\n"; echo "An actuarial model of longevity is employed, in which the probability of death each year is taken from the U.S. Life Tables. Separate tables are used for men and women.\n"; echo "

\n"; echo "All input values and reported values are adjusted for inflation.\n"; echo "

\n"; echo "

\n"; echo "
Equity	Bonds	Success probability \n"; echo "	Mean gifted assets	Mean assets at final death \n"; srand(0); # Simulate portfolios foreach ($equity_fraction as $ef) { if (is_string($ef) && ($ef == 'age-in-bonds')) { $equity = (1 - min(100, ($age2 == -1) ? $age1 : ($age1 + $age2) / 2) / 100) * $total; } else { $equity = $ef * $total; } $cash = $total - $equity; $p_array = array(); $fail = 0; $tot = 0; $charitable_tot = 0; $experiments = $num_experiments; for ($experiment = 0; $experiment < $experiments; $experiment++) { $e = $equity * $pe_adjust; $c = $cash; $in = $income; $charitable = 0; $dead1 = 0; $dead2 = ($age2 == -1); $index1 = rand(0, count($returns) - 1); $index2 = rand(0, count($diversify) - 1); $gone_to_bonds = 0; for ($y = 0; ($y < $years) && !($dead1 && $dead2); $y++) { if ($seq == 'seq') { $e = (1 - $diversify_factor) (1 + $returns[$index1]) + $diversify_factor * (1 + $diversify[$index2]); $index1 = ($index1 + 1) % count($returns); $index2 = ($index2 + 1) % count($diversify); } else { $e = (1 - $diversify_factor) (1 + $returns[rand(0, count($returns) - 1)]) + $diversify_factor * (1 + $diversify[rand(0, count($diversify) - 1)]); } $c = 1 + $ret_bonds / 100; $p = $e + $c; if ($y >= $years_to_retire) { if ($dead1 \|\| ($dead2 && ($age2 != -1))) { $p -= $in $income_dead_fraction; } else { $p -= $in; } } else { $p += $savings; } if (($limit > 0) && ($p > $limit)) { $charitable += $p - $limit; $p = $limit; } if (($go_to_bonds != -1) && ! $gone_to_bonds) { if ($dead1) { $le = $life_expectancy_2[$y]; if ($age1 != -1) { $le = $income_dead_fraction; } } elseif ($dead2) { $le = $life_expectancy_1[$y]; if ($age2 != -1) { $le = $income_dead_fraction; } } else { $le = $life_expectancy_both[$y]; } $gone_to_bonds = ($p >= $le * $in); } if ($gone_to_bonds) { $c = $p; } elseif (is_string($ef) && ($ef == 'age-in-bonds')) { $c = min(100, ($age2 == -1) ? $age1 : ($age1 + $age2) / 2) / 100 * $p; } elseif ($deplete == 'rb') { $c = (1 - $ef) * $p; } else { if ($y < $years_to_retire) { # use savings contribution to try and maintain bond fraction up to entire contribution amount $c = min(max($c, (1 - $ef) * $p), $c + $savings); } if ($deplete == 'bf') { if ($y >= $years_to_retire) { if ($dead1 \|\| ($dead2 && ($age2 != -1))) { $c -= $in * $income_dead_fraction; } else { $c -= $in; } } } else { # ($deplete == 'ef'); $c = $c; # leave alone } } if ($c < 0) { $c = 0; } $e = $p - $c; if ($e < 0) { $c = $p; $e = 0; } if (($age1 != -1) && !$dead1) { $dead1 = (rand(0, 9999) / 10000 < $death1[$age1 + $y]); } if (!$dead2) { $dead2 = (rand(0, 9999) / 10000 < $death2[$age2 + $y]); } } array_push($p_array, $p); if ($p < 0) { $fail++; } $tot += $p; $charitable_tot += $charitable; } # Compute and report stats # sort($p_array); # $pctl_0 = $p_array[0]; # $pctl_10 = $p_array[intval((count($p_array) + 1) / 10)]; # $pctl_25 = $p_array[intval((count($p_array) + 1) / 4)]; # $pctl_50 = $p_array[intval((count($p_array) + 1) / 2)]; $mean = $tot / $experiments; $charitable_mean = $charitable_tot / $experiments; # if ($pctl_0 < 0) { # $pctl_0_display = "-$" . round(- $pctl_0 / 1000) . "k"; # } else { # $pctl_0_display = "$" . round($pctl_0 / 1000) . "k"; # } if ($mean < 0) { $mean_display = "-$" . round(- $mean / 1000) . "k"; } else { $mean_display = "$" . round($mean / 1000) . "k"; } if (is_string($ef) && ($ef == 'age-in-bonds')) { echo "
average age in bonds\n"; } else { printf("
%d%%	%d%%\n", round($ef * 100), round((1 - $ef) * 100)); } printf("	%d%%	$%dk	%s\n", round(($experiments - $fail) / $experiments * 100), round($charitable_mean / 1000), $mean_display); } echo "