#!/usr/bin/php
array(
0.006761, 0.000460, 0.000286, 0.000218, 0.000176, 0.000164, 0.000151, 0.000140, 0.000124, 0.000105,
0.000091, 0.000094, 0.000132, 0.000209, 0.000314, 0.000426, 0.000529, 0.000627, 0.000715, 0.000796,
0.000881, 0.000963, 0.001017, 0.001034, 0.001023, 0.001003, 0.000990, 0.000983, 0.000988, 0.001005,
0.001030, 0.001060, 0.001099, 0.001146, 0.001201, 0.001264, 0.001340, 0.001434, 0.001548, 0.001685,
0.001836, 0.002000, 0.002188, 0.002400, 0.002629, 0.002864, 0.003107, 0.003369, 0.003661, 0.003984,
0.004337, 0.004709, 0.005091, 0.005474, 0.005863, 0.006275, 0.006726, 0.007220, 0.007773, 0.008389,
0.009081, 0.009839, 0.010657, 0.011534, 0.012491, 0.013600, 0.014722, 0.015959, 0.017288, 0.018755,
0.020424, 0.022385, 0.024679, 0.027320, 0.030299, 0.033636, 0.037216, 0.041160, 0.045503, 0.050281,
0.055531, 0.061293, 0.067611, 0.074528, 0.082091, 0.090346, 0.099341, 0.109125, 0.119744, 0.131244,
0.143668, 0.157056, 0.171442, 0.186853, 0.203309, 0.220822, 0.239389, 0.258999, 0.279625, 0.301225,
0.301225, 0.301225, 0.301225, 0.301225, 0.301225, 0.301225, 0.301225, 0.301225, 0.301225, 0.301225, # 100-109 conservatively extrapolated
1, # Terminal value needed by simulate web
),
'male' => array(
0.007390, 0.000490, 0.000316, 0.000242, 0.000201, 0.000182, 0.000170, 0.000156, 0.000134, 0.000107,
0.000085, 0.000089, 0.000143, 0.000256, 0.000411, 0.000573, 0.000725, 0.000873, 0.001014, 0.001149,
0.001292, 0.001427, 0.001512, 0.001529, 0.001497, 0.001448, 0.001409, 0.001382, 0.001376, 0.001390,
0.001412, 0.001437, 0.001474, 0.001516, 0.001570, 0.001634, 0.001716, 0.001821, 0.001956, 0.002120,
0.002303, 0.002505, 0.002735, 0.002992, 0.003270, 0.003556, 0.003855, 0.004187, 0.004570, 0.005001,
0.005474, 0.005969, 0.006473, 0.006971, 0.007469, 0.007995, 0.008567, 0.009179, 0.009843, 0.010571,
0.011378, 0.012264, 0.013227, 0.014275, 0.015434, 0.016771, 0.018156, 0.019682, 0.021327, 0.023144,
0.025204, 0.027616, 0.030417, 0.033598, 0.037153, 0.041097, 0.045315, 0.049944, 0.055019, 0.060576,
0.066655, 0.073296, 0.080542, 0.088435, 0.097021, 0.106343, 0.116446, 0.127371, 0.139160, 0.151850,
0.165475, 0.180063, 0.195635, 0.212205, 0.229779, 0.248348, 0.267897, 0.288394, 0.309795, 0.332043,
0.332043, 0.332043, 0.332043, 0.332043, 0.332043, 0.332043, 0.332043, 0.332043, 0.332043, 0.332043, # 100-109 conservatively extrapolated
1, # Terminal value needed by simulate web
),
'female' => array(
0.006103, 0.000430, 0.000255, 0.000193, 0.000149, 0.000145, 0.000132, 0.000122, 0.000112, 0.000103,
0.000096, 0.000100, 0.000120, 0.000160, 0.000212, 0.000271, 0.000325, 0.000369, 0.000400, 0.000422,
0.000443, 0.000467, 0.000488, 0.000504, 0.000518, 0.000532, 0.000548, 0.000565, 0.000583, 0.000605,
0.000634, 0.000670, 0.000714, 0.000767, 0.000824, 0.000887, 0.000959, 0.001040, 0.001137, 0.001248,
0.001367, 0.001495, 0.001644, 0.001812, 0.001994, 0.002182, 0.002373, 0.002569, 0.002775, 0.002995,
0.003236, 0.003494, 0.003763, 0.004041, 0.004330, 0.004639, 0.004981, 0.005372, 0.005826, 0.006347,
0.006942, 0.007595, 0.008293, 0.009029, 0.009826, 0.010753, 0.011692, 0.012722, 0.013830, 0.015062,
0.016484, 0.018170, 0.020151, 0.022445, 0.025056, 0.028016, 0.031215, 0.034767, 0.038707, 0.043073,
0.047907, 0.053254, 0.059160, 0.065676, 0.072854, 0.080749, 0.089416, 0.098914, 0.109300, 0.120630,
0.132959, 0.146339, 0.160816, 0.176428, 0.193208, 0.211174, 0.230333, 0.250679, 0.272186, 0.294812,
0.294812, 0.294812, 0.294812, 0.294812, 0.294812, 0.294812, 0.294812, 0.294812, 0.294812, 0.294812, # 100-109 conservatively extrapolated
1, # Terminal value needed by simulate web
),
'die_exp' => array( # For testing.
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1, # Terminal value needed by simulate web
),
'male_2003' => array( # Old life table data for comparison purposes.
0.007475, 0.000508, 0.000326, 0.000250, 0.000208, 0.000191, 0.000182, 0.000171, 0.000152, 0.000125, # 0-9
0.000105, 0.000111, 0.000162, 0.000274, 0.000431, 0.000608, 0.000777, 0.000935, 0.001064, 0.001166,
0.001266, 0.001360, 0.001419, 0.001435, 0.001419, 0.001390, 0.001365, 0.001344, 0.001336, 0.001341,
0.001352, 0.001371, 0.001408, 0.001469, 0.001553, 0.001653, 0.001770, 0.001911, 0.002075, 0.002254,
0.002438, 0.002632, 0.002853, 0.003113, 0.003412, 0.003735, 0.004071, 0.004428, 0.004806, 0.005206,
0.005648, 0.006121, 0.006594, 0.007045, 0.007488, 0.007946, 0.008459, 0.009064, 0.009810, 0.010706,
0.011763, 0.012934, 0.014159, 0.015362, 0.016558, 0.017847, 0.019331, 0.020992, 0.022858, 0.024921,
0.027065, 0.029363, 0.032031, 0.035178, 0.038734, 0.042414, 0.046171, 0.050325, 0.055085, 0.060498,
0.066557, 0.072986, 0.079682, 0.086593, 0.094013, 0.102498, 0.111640, 0.121472, 0.132023, 0.143319,
0.155383, 0.168232, 0.181880, 0.196334, 0.211592, 0.227645, 0.244476, 0.262057, 0.280351, 0.299312,
0.32, 0.34, 0.36, 0.38, 0.40, 0.42, 0.44, 0.46, 0.48, 0.50, # 100-109 extrapolated
1, # Terminal value needed by simulate web
),
);
# Sanity check.
(count($death['combined']) == count($death['male'])) or die("Mismatched table lengths.");
(count($death['combined']) == count($death['female'])) or die("Mismatched table lengths.");
(count($death['combined']) == count($death['die_exp'])) or die("Mismatched table lengths.");
(count($death['combined']) == count($death['male_2003'])) or die("Mismatched table lengths.");
# End of data.
function plus_1_geomean($a) {
$prod = 1;
foreach ($a as $i) {
$prod *= 1 + $i;
}
return pow($prod, 1 / count($a));
}
function standard_deviation($a, $sample = TRUE) {
$mean = array_sum($a) / count($a);
$variance = 0;
foreach ($a as $i) {
$variance += pow($i - $mean, 2);
}
$variance /= ($sample ? count($a) - 1 : count($a));
return sqrt($variance);
}
# Returns processing.
function load_shiller_data($shiller_file) {
global $stock, $bond;
$stock = array();
$bond = array();
$recent_d = array();
$lines = file($shiller_file);
$processed = 0;
foreach ($lines as $line) {
$fields = explode(',', $line);
$year_month = $fields[0];
$nominal_p = $fields[1];
$nominal_d = $fields[2];
$_ = $fields[3];
$cpi = $fields[4];
$_ = $fields[5];
$nominal_interest_rate = $fields[6];
$year = intval($year_month);
$month = round(($year_month - $year) * 100);
if (($year == 0) || ($month == 0)) {
continue;
}
if (($nominal_d == '') || ($nominal_p == '')) {
continue;
}
$nominal_p = floatval($nominal_p);
$nominal_d = floatval($nominal_d);
$cpi = floatval($cpi);
$nominal_interest_rate = floatval($nominal_interest_rate) / 100;
$p = $nominal_p / $cpi;
$d = $nominal_d / $cpi;
#echo "$year $month $nominal_p $nominal_d $p $d $cpi $nominal_interest_rate\n";
array_push($recent_d, $d / 12);
if ($processed >= 12) {
array_shift($recent_d);
}
if ($month == 12) {
if ($processed >= 12) {
array_push($stock, ($p + array_sum($recent_d)) / $old_p - 1);
}
$old_p = $p;
}
if ($month == 6) {
if ($processed >= 12) {
$interest_rate = (1 + $nominal_interest_rate) / ($cpi / $old_cpi) - 1;
array_push($bond, $interest_rate);
}
$old_cpi = $cpi;
}
$processed++;
}
$min_len = min(count($stock), count($bond));
$stock = array_slice($stock, 0, $min_len);
$bond = array_slice($bond, 0, $min_len);
}
function load_sbbi_data($sbbi_file) {
global $stock, $bond;
$stock = array();
$bond = array();
$lines = file($sbbi_file);
foreach ($lines as $line) {
if (preg_match('/^\s*#/', $line)) {
continue;
}
$fields = explode(',', $line);
$year = $fields[0];
$cpi_delta = $fields[1];
$_ = $fields[2];
$_ = $fields[3];
$nominal_b_delta = $fields[4];
$nominal_e_delta = $fields[5];
$year = intval($year);
$nominal_e_delta = floatval($nominal_e_delta);
$nominal_b_delta = floatval($nominal_b_delta);
$cpi_delta = floatval($cpi_delta);
$e = (1 + $nominal_e_delta / 100) / (1 + $cpi_delta / 100) - 1;
$b = (1 + $nominal_b_delta / 100) / (1 + $cpi_delta / 100) - 1;
array_push($stock, $e);
array_push($bond, $b);
}
}
function massage_returns($ret_equity, $ret_bonds, $diversify_sds, $initial_year, $sim_start_year, $sim_end_year) {
global $stock, $bond;
$start = $sim_start_year - $initial_year;
$count = (is_null($sim_end_year) ? count($stock) - $start : $sim_end_year - $sim_start_year + 1);
$returns = array_slice($stock, $start, $count);
if (is_null($ret_bonds)) {
$bond_returns = array_slice($bond, $start, $count);
} else {
$bond_returns = array();
for ($i = 0; $i < $count; $i++) {
array_push($bond_returns, $ret_bonds);
}
}
# 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 = (is_null($ret_equity) ? pow($samp_prod, 1 / count($diversify)) - 1 : 1);
$adjust = (1 + $ret_equity) / (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 = (is_null($ret_equity) ? 1 : (1 + $ret_equity) / (1 + $samp_geomean));
for ($i = 0; $i <= count($returns) - 1; $i++) {
$equity_returns[$i] = (1 + $returns[$i]) * $adjust - 1;
}
return array('equity_returns' => $equity_returns, 'bond_returns' => $bond_returns, 'diversify' => $diversify);
}
# Simulation engine.
function simulate_web($total, $ef, $limit, $savings, $income, $ret_equity, $ret_bonds, $pe_adjust, $sim_start_year, $seq, $deplete, $num_experiments,
$age1, $sex1, $age2, $sex2, $years, $years_to_retire, $go_to_bonds, $income_dead_fraction, $diversification, $correlation, $diversify_sds) {
global $death, $shiller_initial;
$res = massage_returns($ret_equity, $ret_bonds, $diversify_sds, $shiller_initial, $sim_start_year, NULL);
$equity_returns = $res['equity_returns'];
$bond_returns = $res['bond_returns'];
$diversify = $res['diversify'];
$death1 = $death[$sex1];
$death2 = $death[$sex2];
# 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); $y++) {
$alive1 = 1;
$alive2 = 1;
$death_prob[$y] = array();
for ($e = $y; min($age1, $age2) + $e < count($death1); $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); $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); $y++) {
$alive1 = 1;
$death_prob[$y] = array();
for ($e = $y; $age1 + $e < count($death1); $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); $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;
}
}
}
# Simulate.
$diversify_factor = (1 - $correlation) * $diversification;
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;
$result_array = array();
$fail = 0;
$p_array = array();
$charitable_array = array();
$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($equity_returns) - 1);
$index2 = rand(0, count($diversify) - 1);
$gone_to_bonds = 0;
for ($y = 0; ($y < $years) && !($dead1 && $dead2); $y++) {
$e *= (1 - $diversify_factor) * (1 + $equity_returns[$index1]) + $diversify_factor * (1 + $diversify[$index2]);
$c *= 1 + $bond_returns[$index1];
if ($seq == 'seq') {
$index1 = ($index1 + 1) % count($equity_returns);
$index2 = ($index2 + 1) % count($diversify);
} else {
$index1 = rand(0, count($equity_returns) - 1);
$index2 = rand(0, count($diversify) - 1);
}
$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);
array_push($charitable_array, $charitable);
if ($p < 0) {
$fail++;
}
}
array_multisort($p_array, SORT_ASC, $charitable_array, SORT_ASC);
return array('success' => ($experiments - $fail) / $experiments, 'p_array' => $p_array, 'charitable_array' => $charitable_array);
}
#function cli_simulate($age, $wr, $ef) {
#
# global $years;
#
# $total = 1;
#
# $limit = 0;
# $savings = 0;
# $income = $wr;
#
# $ret_equity = 0.05;
# $ret_bonds = 0.02;
# $pe_adjust = 1;
# $sim_start_year = 1872;
# $seq = 'seq';
# $deplete = 'rb';
# $num_experiments = 1000;
#
# $age1 = $age;
# $sex1 = 'male';
# $age2 = -1;
# $sex2 = 'female';
# # $years = ...
# $years_to_retire = 0;
# $go_to_bonds = -1;
# $income_dead_fraction = 0.7;
# $diversification = 0;
# $correlation = 0.8;
# $diversify_sds = 1.3;
#
# return simulate_web($total, $ef, $limit, $savings, $income, $ret_equity, $ret_bonds, $pe_adjust, $sim_start_year, $seq, $deplete, $num_experiments,
# $age1, $sex1, $age2, $sex2, $years, $years_to_retire, $go_to_bonds, $income_dead_fraction, $diversification, $correlation, $diversify_sds);
#}
#function wr_ef($age, $pctl_success) {
#
# global $equity_range;
#
# $max_wr = -1;
# $max_ef = -1;
#
# foreach ($equity_range as $ef) {
#
# # Binary search.
# $wr_low = 0;
# $wr_high = 0.10;
# $wr_delta = 0.00001;
# $success_delta = 0.005;
# $res = cli_simulate($age, $wr_low, $ef);
# assert($res['success'] > $pctl_success);
# $res = cli_simulate($age, $wr_high, $ef);
# assert($res['success'] < $pctl_success);
# while (TRUE) {
# $wr_mid = ($wr_low + $wr_high) / 2;
# if (($wr_high - $wr_low) / 2 <= $wr_delta) {
# echo "Binary search failed - skipping: $age $pctl_success $ef\n";
# continue 2;
# }
# $res = cli_simulate($age, $wr_mid, $ef);
# if (abs($res['success'] - $pctl_success) < $success_delta) {
# break;
# }
# if ($res['success'] > $pctl_success) {
# $wr_low = $wr_mid;
# } else {
# $wr_high = $wr_mid;
# }
# }
#
# if ($wr_mid > $max_wr) {
# $max_wr = $wr_mid;
# $max_ef = $ef;
# }
# }
#
# assert($max_swr != -1);
#
# return array('wr' => $max_wr, 'ef' => $max_ef);
#}
function pre_compute_life_expectancy($death1) {
$le = array();
for ($s = 0; $s < count($death1); $s++) {
$le[$s] = 0;
$alive = 1;
for ($i = 0; $y = $s + $i, $y < count($death1); $i++) {
$prev_alive = $alive;
$alive *= (1 - $death1[$y]);
$avg_alive = ($alive + $prev_alive) / 2;
$le[$s] += $avg_alive;
}
}
return $le;
}
#function search_ef($equity_range, $savings, &$future, $future_base, $bucket_min_wf, $min_wf, $wf,
# &$le, &$returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
# $total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps) {
#
# if (is_null($ef_search_steps)) {
# return search_ef_all($equity_range, $savings, $future, $future_base, $bucket_min_wf, $min_wf, $wf,
# $le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
# $total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
# } else {
# return search_ef_fast($equity_range, $savings, $future, $future_base, $bucket_min_wf, $min_wf, $wf,
# $le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
# $total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
# }
#}
#
#function search_ef_all($equity_range, $savings, &$future, $future_base, $bucket_min_wf, $min_wf, $wf,
# &$le, &$returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
# $total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps) {
#
# $max_wf = NULL;
# $max_ef = NULL;
# $all_success = TRUE;
# $max_success = NULL;
# $success_prob_array = array();
# #$target_success = pow($total_pctl_success, $le[$age] / $le[$original_age]); This doesn't work as expected since it limits the options considered.
# $target_success = $total_pctl_success;
# $prev_success = $target_success;
# foreach (array_reverse($equity_range) as $ef) {
# $age_array = array($age);
# $wr = 1 * $wf;
# $single_wr_range = array($wr);
# $single_equity_range = array($ef);
# $res = simulate(1, $equity_range, $single_equity_range, $savings, $single_wr_range, $future, $future_base, $min_wf,
# $le, $returns, $deplete, $original_age, $age_array, $sex, $years, $years_to_retire, $wrap,
# $total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
# $success = (($success_mode == 'stw') ? $res['stw_success'] : $res['success']);
# $ef_array = $success[$age][strval($wr)];
# if (array_key_exists(strval($ef), $ef_array)) {
# $curr_success = $ef_array[strval($ef)];
# $success_prob_array[strval($ef)] = $curr_success;
# if ($curr_success < 1) {
# $all_success = FALSE;
# }
# # Select the ef with the highest success probability. We favor the more aggressive of two equal ef values for the non-success based expected upside.
# if (($curr_success > $prev_success) || (is_null($max_success) && ($curr_success == $prev_success))) {
# $max_wf = $wr / 1;
# $max_ef = $ef;
# $max_success = $curr_success;
# $prev_success = $curr_success;
# } else {
# # It might be expected that we can early exit the loop if we find success probabilities are descending.
# # In practice we find this is not the case, probabilities can rise again, by as much as 1.5% over a 0.25 equity interval.
# # It is unclear if this behavour is real and the result of our use of real world data or or the result of some unknown bug.
# # Trying to set a small delta tolerance we find it only speeds things up if the tolerance is at the point where
# # it discards useful records and effects the end result.
# #if (!is_null($max_success) && ($curr_success < $prev_success - 0.001)) {
# # break; # Assume success probability doesn't go down then up as a function of ef (up then down as a function of ef reversed).
# #}
# }
# } else {
# $all_success = FALSE;
# }
# }
# assert(!$all_success || ($curr_success == 1));
#
# return array('bucket' => $bucket_min_wf, 'wf' => $max_wf, 'ef' => $max_ef, 'all_success' => $all_success, 'success_prob' => $max_success,
# 'method' => 'calculated all for bucket', 'success' => $success_prob_array);
# # Record bucket mwr, method, and success for analysis and debugging purposes.
#}
function search_ef($equity_range, $savings, &$future, $future_base, $bucket_min_wf, $min_wf, $wf,
&$le, &$returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps) {
# Discover ef_index, the key for ef in the equity_range array.
if (is_null($ef_search_steps)) {
$ef_search_steps = count($equity_range) - 1;
$ef_index = count($equity_range) - 1;
} else {
# Search for hint in older bucket.
$hint = lookup_bucket_soft($future, $future_base, NULL, $bucket_min_wf, $age + 1); # Float rounding errors in bucket_min_wf don't matter, this is only a hint.
$scan_limit = 20;
# Search for hint in buckets above and below.
if (is_null($hint)) {
for ($offset = 1; $offset < $scan_limit; $offset++) {
$hint = lookup_bucket_soft($future, $future_base, NULL, $bucket_min_wf * pow($future_base, $offset), $age);
if (!is_null($hint)) {
break;
}
$hint = lookup_bucket_soft($future, $future_base, NULL, $bucket_min_wf / pow($future_base, $offset), $age);
if (!is_null($hint)) {
break;
}
}
}
if (is_null($hint)) {
$ef_index = count($equity_range) - 1;
} else {
foreach ($equity_range as $ef_index => $ef) {
if ($ef == $hint['ef']) {
break;
}
}
assert($ef == $hint['ef']);
}
}
$age_array = array($age);
$wr = 1 * $wf;
$single_wr_range = array($wr);
$max_wf = NULL;
$max_ef = NULL;
$max_success = NULL;
$all_success = TRUE;
$success_prob_array = array();
$simulated = array();
$index_start = max($ef_index - $ef_search_steps, 0);
$index_end = min($ef_index + $ef_search_steps, count($equity_range) - 1);
for ($index = $index_start; $index <= $index_end; $index++) {
if (array_key_exists($index, $simulated)) {
continue;
}
$ef = $equity_range[$index];
$single_equity_range = array($ef);
$res = simulate(1, $equity_range, $single_equity_range, $savings, $single_wr_range, $future, $future_base, $min_wf,
$le, $returns, $deplete, $original_age, $age_array, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
$simulated[$index] = TRUE;
$success = (($success_mode == 'stw') ? $res['stw_success'] : $res['success']);
$ef_array = $success[$age][strval($wr)];
if (array_key_exists(strval($ef), $ef_array)) {
$curr_success = $ef_array[strval($ef)];
$success_prob_array[strval($ef)] = $curr_success;
if ($curr_success < 1) {
$all_success = FALSE;
}
if (is_null($max_success) ? ($curr_success >= $total_pctl_success) : ($curr_success >= $max_success)) {
$ef_index = $index;
$max_wf = $wr / 1;
$max_ef = $ef;
$max_success = $curr_success;
$index = max($ef_index - $ef_search_steps, 0) - 1; # -1 because value will be incremented by loop.
$index_end = min($ef_index + $ef_search_steps, count($equity_range) - 1);
}
} else {
$all_success = FALSE;
}
}
if ($all_success) {
# Verify all_success by checking the end points.
$end_points = array(0, count($equity_range) - 1);
foreach ($end_points as $index) {
if (!array_key_exists($index, $simulated)) {
$ef = $equity_range[$index];
$single_equity_range = array($ef);
$res = simulate(1, $equity_range, $single_equity_range, $savings, $single_wr_range, $future, $future_base, $min_wf,
$le, $returns, $deplete, $original_age, $age_array, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
$simulated[$index] = TRUE;
$success = (($success_mode == 'stw') ? $res['stw_success'] : $res['success']);
$ef_array = $success[$age][strval($wr)];
if (array_key_exists(strval($ef), $ef_array)) {
$curr_success = $ef_array[strval($ef)];
$success_prob_array[strval($ef)] = $curr_success;
if ($curr_success < 1) {
$all_success = FALSE;
break;
}
if ($index == count($equity_range) - 1) {
$max_wf = $wr / 1;
$max_ef = $ef;
$max_success = $curr_success;
}
} else {
$all_success = FALSE;
}
}
}
}
ksort($success_prob_array);
$max_er_success = (is_null($max_success) ? NULL : pow($max_success, $le[$original_age] / $le[$age]));
return array('bucket' => $bucket_min_wf, 'wf' => $max_wf, 'ef' => $max_ef, 'all_success' => $all_success,
'success_prob' => $max_success, 'er_success' => $max_er_success, 'method' => 'calculated for bucket', 'success' => $success_prob_array);
}
function buckets_equal($a, $b) {
return ((is_null($a['wf']) == is_null($b['wf'])) &&
($a['ef'] == $b['ef']) &&
($a['all_success'] == $b['all_success']));
}
function new_bucket($equity_range, $savings, &$future, $future_base, $bucket,
&$le, &$returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps) {
global $allow_interpolate;
$bucket_min_wf = exp($bucket * log($future_base));
$smallest = NULL;
$biggest = NULL;
if (count($future[$age]) > 0) {
$future_keys = array_keys($future[$age]);
$biggest = max($future_keys);
$smallest = min($future_keys);
if (($bucket > $biggest) && is_null($future[$age][$biggest]['wf'])) {
return array('bucket' => $bucket_min_wf, 'wf' => NULL, 'ef' => NULL, 'all_success' => FALSE, 'success_prob' => 0, 'er_success' => 0,
'method' => 'above biggest');
} else if (($bucket < $smallest) && $future[$age][$smallest]['all_success']) {
$min_wf = $bucket_min_wf;
return array('bucket' => $bucket_min_wf, 'wf' => $min_wf, 'ef' => $equity_range[count($equity_range) - 1], 'all_success' => TRUE,
'success_prob' => 1, 'er_success' => 1, 'method' => 'below smallest');
}
}
# See if buckets some distance above and below exist and have the same ef value, if so use that.
# This assumes it is safe to do this.
$below = NULL;
$above = NULL;
if (count($future[$age]) > 0) {
for ($try_below = $bucket - 1; $try_below >= $smallest; $try_below--) {
if (array_key_exists($try_below, $future[$age])) {
$below = $try_below;
break;
}
}
for ($try_above = $bucket + 1; $try_above <= $biggest; $try_above++) {
if (array_key_exists($try_above, $future[$age])) {
$above = $try_above;
break;
}
}
}
if (!is_null($below) && !is_null($above)) {
if (buckets_equal($future[$age][$below], $future[$age][$above])) {
# We interpolate outside the region of interest even when allow_interpolate is false. Doing so produces a 4 fold performance win.
if ($future[$age][$below]['success_prob'] == 0) {
$success_prob = 0;
} else if ($future[$age][$above]['all_success']) {
$success_prob = 1;
} else {
$success_prob = NULL;
}
if ($allow_interpolate || !is_null($success_prob)) {
$min_wf = (is_null($future[$age][$below]['wf']) ? NULL : $bucket_min_wf);
return array('bucket' => $bucket_min_wf, 'wf' => $min_wf, 'ef' => $future[$age][$below]['ef'],
'all_success' => $future[$age][$below]['all_success'], 'success_prob' => $success_prob, 'er_success' => $success_prob,
'method' => 'interpolated constant ef');
}
}
}
# Compute recommendations for our bucket.
$wf = $bucket_min_wf;
$min_wf = $wf;
$bucket_res = search_ef($equity_range, $savings, $future, $future_base, $bucket_min_wf, $min_wf, $wf,
$le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
# Pre-populate cache down to next change in ef.
if (is_null($below)) {
# Exponentially search down.
$prev_int_i = 0;
$res = $bucket_res;
$start = (is_null($smallest) ? $bucket + 1 : $smallest);
for ($i = 2 * ($start - $bucket); !$res['all_success']; $i *= 1.4) { # Adjust to improve performance.
$int_i = intval($i);
if ($int_i == $prev_int_i) {
continue;
}
$prev_int_i = $int_i;
$bottom = $start - $int_i;
$wf = pow($future_base, $bottom);
$min_wf = $wf;
$bucket_min_wf = $wf;
$res = search_ef($equity_range, $savings, $future, $future_base, $bucket_min_wf, $min_wf, $wf,
$le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
$future[$age][$bottom] = $res;
}
} else {
if (buckets_equal($future[$age][$below], $bucket_res)) {
# No need to pre-populate value below matches current value. It is responsible for pre-populating down.
# Empty.
} else {
# Binary search up.
$low = $below;
$high = $bucket;
while ($low < $high - 1) {
$mid = round(($high + $low) / 2);
$wf = pow($future_base, $mid);
$min_wf = $wf;
$bucket_min_wf = $wf;
$res = search_ef($equity_range, $savings, $future, $future_base, $bucket_min_wf, $min_wf, $wf,
$le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
$future[$age][$mid] = $res;
if (buckets_equal($res, $bucket_res)) {
$high = $mid;
} else {
$low = $mid;
}
}
}
}
return $bucket_res;
}
# future_age[$age] is indexed by bucket since the optimal strategy may depend on the current min_wf value.
function lookup_bucket($equity_range, $savings, &$future, $future_base, $bucket, $min_wf,
&$le, &$returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps) {
assert(is_null($bucket) xor is_null($min_wf));
# Speedup assuming 100% returns aren't possible.
if (is_null($min_wf)) {
$min_wf = exp($bucket * log($future_base));
}
if ($min_wf >= 2) {
return array('wf' => NULL, 'ef' => NULL, 'all_success' => FALSE, 'method' => 'impossible wr');
}
if (!array_key_exists($age, $future)) {
$future[$age] = array();
}
$future_age = &$future[$age];
if (is_null($bucket)) {
$bucket = intval(ceil(log($min_wf) / log($future_base)));
}
if (!array_key_exists($bucket, $future_age)) {
$future_age[$bucket] = new_bucket($equity_range, $savings, $future, $future_base, $bucket,
$le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
}
return $future_age[$bucket];
}
function lookup_bucket_soft(&$future, $future_base, $bucket, $min_wf, $age) {
assert(is_null($bucket) xor is_null($min_wf));
if (is_null($bucket)) {
$bucket = intval(ceil(log($min_wf) / log($future_base)));
}
if (!array_key_exists($age, $future)) {
return NULL;
}
$future_age = $future[$age];
if (!array_key_exists($bucket, $future_age)) {
return NULL;
}
return $future_age[$bucket];
}
# Helper for simulate().
function compute_success($death1, &$le, $all_solvent, &$solvent_prob, $years, $age) {
if ($all_solvent) {
return array('success_prob' => 1, 'insolvent_years' => 0, 'total_years' => $le[$age]);
}
$success_prob = 0;
$alive = 1;
$solvent = 1;
$insolvent_years = 0;
$total_years = 0;
$loop_completed = TRUE;
$limit = min($years, count($death1) - $age);
for ($i = 0; $i < $limit; $i++) {
$solvent = (array_key_exists($i, $solvent_prob) ? $solvent_prob[$i] : 0);
if ($solvent == 0) {
$loop_completed = FALSE;
$insolvent_years += $alive * $le[$age + $i];
$total_years += $alive * $le[$age + $i];
break;
}
$dying = $alive * $death1[$age + $i];
$prev_alive = $alive;
$alive -= $dying;
$avg_alive = ($alive + $prev_alive) / 2;
$success_prob += $dying * $solvent;
$insolvent_years += $avg_alive * (1 - $solvent);
$total_years += $avg_alive;
}
if ($loop_completed) {
$success_prob += $alive * $solvent;
}
return array('success_prob' => $success_prob, 'insolvent_years' => $insolvent_years, 'total_years' => $total_years);
}
# Compute $success[$age][$wr][$ef] precisely.
function simulate($total, $equity_range, $doit_equity_range, $savings, $wr_range, &$future, $future_base, $min_wf,
&$le, &$returns, $deplete, $original_age, $age_range, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps) {
global $death;
assert(is_null($future) || (count($age_range) == 1));
$equity_returns = $returns['equity_returns'];
$bond_returns = $returns['bond_returns'];
$death1 = $death[$sex];
$actual_years = (is_null($years) ? count($death1) : $years);
$first_age = $age_range[0];
$max_years = min($actual_years, count($death1) - $first_age);
# Skip wr values below min_wf.
if (!is_null($future)) {
$new_wr_range = array();
foreach ($wr_range as $wr) {
if ($wr / $total >= $min_wf) {
array_push($new_wr_range, $wr);
}
}
$wr_range = $new_wr_range;
}
foreach ($age_range as $age) {
$success[$age] = array();
foreach ($wr_range as $wr) {
$success[$age][strval($wr)] = array();
$stw_success[$age][strval($wr)] = array();
}
}
$len_equity_returns = count($equity_returns); # 25% speedup.
if ($wrap) {
$num_equity_returns = $len_equity_returns;
} else {
$num_equity_returns = $len_equity_returns - $max_years + 1;
assert($num_equity_returns >= 0);
}
foreach ($doit_equity_range as $initial_ef) {
$age_last_safe = 1000000;
foreach (array_reverse($wr_range) as $wr) {
# Compute $solvent_prob: map from offset year to prob being solvent in that year across all return sequences.
$equity = $initial_ef * $total;
$cash = $total - $equity;
$min_wr = $min_wf * $total;
$solvent_prob = array();
$all_solvent = TRUE;
for ($s = 0; $s < $num_equity_returns; $s++) {
$e = $equity;
$c = $cash;
$wf = $wr / $total;
$scaled_min_wf = $min_wf;
$ef = $initial_ef;
$prev_p = $total;
$all_success = FALSE;
for ($y = 0, $index = $s; $y < $max_years; $y++, $index = ($index + 1) % $len_equity_returns) {
if ($all_success) {
if (!array_key_exists($y, $solvent_prob)) {
$solvent_prob[$y] = 0;
}
$solvent_prob[$y] += 1;
continue;
}
$e *= 1 + $equity_returns[$index];
$c *= 1 + $bond_returns[$index];
$p = $e + $c;
#echo "lo: $age $y wr=$wr ef=$ef prev_p=$prev_p p=$p $s\n";
if ($y >= $years_to_retire) {
if (is_null($future)) {
$p -= $wr;
} else {
if (is_null($wf)) {
$p = NULL;
} else {
$actual_wr = max($wf * $prev_p, $min_wr);
$p -= $actual_wr;
if ($p > 0) {
$scaled_min_wf = $min_wf * $total / $p;
} else {
$scaled_min_wf = NULL;
}
#echo "dd: $age $y wf=$wf min_wf=$min_wf actual_wr=$actual_wr $s\n";
#echo "sm: $age $y p=$p smw=$scaled_min_wf $s\n";
if (($y + 1 < $max_years) && !is_null($scaled_min_wf)) {
$new_age = $age_range[0] + $y + 1;
$future_years = $actual_years - $y - 1;
$res = lookup_bucket($equity_range, $savings, $future, $future_base, NULL, $scaled_min_wf,
$le, $returns, $deplete, $original_age, $new_age, $sex, $future_years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
$wf = $res['wf'];
$ef = $res['ef'];
assert(is_null($wf) == is_null($ef));
if (is_null($wf)) {
# Just because the bucket was unsuccessful, doesn't mean we will be. Our wf is slightly smaller than the bucket wf, so we might succeed.
$wf = $min_wf;
$ef = $equity_range[count($equity_range) - 1]; # Use most aggresive asset allocation. Higher chance of success in the long run.
# Should we use $equity_range[0] if we are close to death? Or just increase the bucket resolution.
} else {
$wf *= $p / $total;
#$wf = $min_wf; # The bucket withdrawal recommendation is at the bucket ceiling.
# Could experiment with taking a tad less like this. No noticeable difference in results.
}
$old_all_success = $all_success;
$all_success = $res['all_success'];
} else {
$wf = NULL;
$ef = NULL;
$all_success = FALSE;
}
#echo "if: $age $y wf=$wf ef=$ef $s\n";
}
}
} else {
$p += $savings;
}
if (!is_null($p) && !is_null($ef)) {
if ($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) {
$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 (is_null($p)) {
$all_solvent = FALSE;
break;
}
if (!array_key_exists($y, $solvent_prob)) {
$solvent_prob[$y] = 0;
}
if ($p >= 0) {
$solvent_prob[$y] += 1;
} else {
# Interpolate solvency in year of bankruptcy to produce smother transitions as $ef changes, and ultimately a smoother SWR EF versus age graph.
$solvent_prob[$y] += $prev_p / ($prev_p - $p);
$all_solvent = FALSE;
break;
}
$prev_p = $p;
}
}
# We scale back rather than adding the fractional probabilities to begin with to prevent rounding errors.
foreach ($solvent_prob as $y => $_) {
$solvent_prob[$y] /= $num_equity_returns;
}
$skip_set_age_last_safe = FALSE;
foreach (array_reverse($age_range) as $age) {
if ($age >= $age_last_safe) {
$success_prob = 1;
$stw_success_prob = 1;
$total_years = $le[$age];
$insolvent_years = 0;
} else {
$res = compute_success($death1, $le, $all_solvent, $solvent_prob, $actual_years, $age);
$success_prob = $res['success_prob'];
$insolvent_years = $res['insolvent_years'];
$total_years = $res['total_years'];
if ($total_years == 0) {
$stw_success_prob = 1;
} else {
$stw_success_prob = 1 - $insolvent_years / $total_years;
}
if (!$skip_set_age_last_safe) {
if (($success_prob == 1) && ($stw_success_prob == 0)) {
$age_last_safe = $age;
} else {
$skip_set_age_last_safe = TRUE;
}
}
}
if (($success_prob < $min_pctl_success) && ($stw_success_prob < $min_pctl_stw_success)) {
# Can only safely early exit if $death1 is monotone increasing.
break;
}
$success[$age][strval($wr)][strval($initial_ef)] = $success_prob;
$stw_success[$age][strval($wr)][strval($initial_ef)] = $stw_success_prob;
}
}
}
return array('success' => $success, 'stw_success' => $stw_success);
}
function precise_wr_ef($success, $age_range, $wr_range, $equity_range, $pctl_success_range) {
$safe = array();
foreach ($age_range as $age) {
$safe[$age] = array();
$pctl_success_array = $pctl_success_range;
$pctl_success = array_shift($pctl_success_array);
$first_iteration = TRUE;
foreach (array_reverse($wr_range) as $wr) {
$ef_array = $success[$age][strval($wr)];
$max_success = -1;
foreach ($equity_range as $ef) {
if (array_key_exists(strval($ef), $ef_array)) {
$curr_success = $ef_array[strval($ef)];
if ($curr_success > $max_success) {
$max_success = $curr_success;
$max_wr = $wr;
$max_ef = $ef;
}
}
}
while ($max_success >= $pctl_success) {
if ($first_iteration) {
$safe[$age][strval($pctl_success)] = NULL;
} else {
$safe[$age][strval($pctl_success)] = array('wr' => $max_wr, 'ef' => $max_ef);
}
$pctl_success = array_shift($pctl_success_array);
if (is_null($pctl_success)) {
break 2;
}
}
$first_iteration = FALSE;
}
while (!is_null($pctl_success)) {
$safe[$age][strval($pctl_success)] = NULL;
$pctl_success = array_shift($pctl_success_array);
}
}
return $safe;
}
# Time varying withdrawal rates and equity fractions.
function precise_tv($equity_range, $savings, $future_base, $gnuplot_bottom_bucket, $gnuplot_top_bucket,
$le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps) {
global $death;
assert($years_to_retire == 0);
$future = array();
# Main code.
for ($min_wf = pow($future_base, intval(log(1.2) / log($future_base))); $min_wf <= 1.2; $min_wf /= $future_base) {
# Search down from bucket closest to 120%. Must be above the actual mwr. Setting it too high increases the run time slightly.
$record = lookup_bucket($equity_range, $savings, $future, $future_base, NULL, $min_wf,
$le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
if ($record['all_success']) {
break;
}
}
# Fill in any missing data for plotting.
foreach ($future as $age => $_) {
for ($bucket = $gnuplot_bottom_bucket; $bucket <= $gnuplot_top_bucket; $bucket++) {
$record = lookup_bucket($equity_range, $savings, $future, $future_base, $bucket, NULL,
$le, $returns, $deplete, $original_age, $age, $sex, $years, $years_to_retire, $wrap,
$total_pctl_success, $min_pctl_success, $min_pctl_stw_success, $success_mode, $ef_search_steps);
}
}
return $future;
}
function dump_data($prefix, $success, $safe, $pctl_success_range) {
global $age_range, $equity_range;
# Dump success.
$f = fopen("$prefix-success.data", "w");
foreach ($success as $age => $wr_array) {
foreach ($wr_array as $wr => $ef_array) {
foreach ($ef_array as $ef => $success_prob) {
fwrite($f, "$age $wr $ef $success_prob\n");
}
}
}
fclose($f);
# Graph safe withdrawal rates and equity fraction.
foreach ($pctl_success_range as $pctl_success) {
$f = fopen("$prefix-swr-$pctl_success.data", "w");
foreach ($age_range as $age) {
$res = $safe[$age][strval($pctl_success)];
if (is_null($res)) {
$wr = '-';
$ef = '-';
$success_prob = '-';
} else {
$wr = $res['wr'];
$ef = $res['ef'];
$success_prob = $success[$age][strval($wr)][strval($ef)];
}
fwrite($f, "$age $wr $ef $success_prob\n");
}
fclose($f);
}
## Graph success probability as a function of age.
#foreach (array(0.03, 0.04, 0.05) as $wr) {
# foreach (array(0, 0.2, 0.4, 0.6, 0.8, 1) as $ef) {
# $f = fopen("$prefix-success-$wr-$ef.data", "w");
# foreach ($success as $age => $wr_array) {
# if (!array_key_exists(strval($wr), $wr_array)) {
# continue;
# }
# $ef_array = $wr_array[strval($wr)];
# if (array_key_exists(strval($ef), $ef_array)) {
# $success_prob = $ef_array[strval($ef)];
# fwrite($f, "$age $wr $ef $success_prob\n");
# }
# }
# fclose($f);
# }
#}
# Scatter plot equity fraction as a function of age for near each max percentile success.
foreach (array(0.001, 0.002, 0.003) as $wr_delta) {
foreach ($pctl_success_range as $pctl_success) {
$f = fopen("$prefix-ef-$wr_delta-$pctl_success.data", "w");
foreach ($success as $age => $wr_array) {
$res = $safe[$age][strval($pctl_success)];
if (is_null($res)) {
continue;
}
$wr = $res['wr'];
$ef = $res['ef'];
$delta_dec_wr = round($wr - $wr_delta, 6);
if (array_key_exists(strval($delta_dec_wr), $wr_array)) {
foreach ($wr_array[strval($delta_dec_wr)] as $ef => $success_prob) {
if ($success_prob >= $pctl_success) {
fwrite($f, "$age $delta_dec_wr $ef $success_prob\n");
}
}
}
}
fclose($f);
}
}
## Graph the sensitivity of success to changes in withdrawal rate near the safe withdrawal rate.
#foreach (array(0.001, 0.002, 0.003) as $wr_delta) {
# foreach ($pctl_success_range as $pctl_success) {
# $f = fopen("$prefix-success-wr-$wr_delta-$pctl_success.data", "w");
# foreach ($success as $age => $wr_array) {
# $res = $safe[$age][strval($pctl_success)];
# if (is_null($res)) {
# continue;
# }
# $wr = $res['wr'];
# $ef = $res['ef'];
# $delta_dec_wr = round($wr - $wr_delta, 6);
# if (array_key_exists(strval($delta_dec_wr), $wr_array)) {
# if (array_key_exists(strval($ef), $wr_array[strval($delta_dec_wr)])) {
# $success_prob = $wr_array[strval($delta_dec_wr)][strval($ef)];
# fwrite($f, "$age $delta_dec_wr $ef $success_prob\n");
# }
# }
# }
# fclose($f);
# }
#}
## Graph sensitivity of success to changes in equity fraction near the safe withdrawal rate.
#foreach ($pctl_success_range as $pctl_success) {
# $f = fopen("$prefix-sef-$pctl_success.data", "w");
# foreach ($age_range as $age) {
# $res = $safe[$age][strval($pctl_success)];
# if (is_null($res)) {
# continue;
# }
# $wr = $res['wr'];
# $ef = $res['ef'];
# $ef_array = $success[$age][strval($wr)];
# $init_success = $ef_array[strval($ef)];
# $delta1_inc_ef = min(1, round($ef + 0.1, 3));
# $delta1_dec_ef = max(0, round($ef - 0.1, 3));
# if (array_key_exists(strval($delta1_inc_ef), $ef_array)) {
# $delta1_inc_success = $ef_array[strval($delta1_inc_ef)];
# } else {
# $delta1_inc_success = NULL;
# }
# if (array_key_exists(strval($delta1_dec_ef), $ef_array)) {
# $delta1_dec_success = $ef_array[strval($delta1_dec_ef)];
# } else {
# $delta1_dec_success = NULL;
# }
# $delta1_success = ((is_null($delta1_inc_success) || is_null($delta1_dec_success)) ? '-' : min($delta1_inc_success, $delta1_dec_success));
# $delta2_inc_ef = min(1, round($ef + 0.2, 3));
# $delta2_dec_ef = max(0, round($ef - 0.2, 3));
# if (array_key_exists(strval($delta2_inc_ef), $ef_array)) {
# $delta2_inc_success = $ef_array[strval($delta2_inc_ef)];
# } else {
# $delta2_inc_success = NULL;
# }
# if (array_key_exists(strval($delta2_dec_ef), $ef_array)) {
# $delta2_dec_success = $ef_array[strval($delta2_dec_ef)];
# } else {
# $delta2_dec_success = NULL;
# }
# $delta2_success = ((is_null($delta2_inc_success) || is_null($delta2_dec_success)) ? '-' : min($delta2_inc_success, $delta2_dec_success));
# $deltap_success = (count($ef_array) == count($equity_range) ? min($ef_array) : '-');
# fwrite($f, "$age $pctl_success $init_success $delta1_success $delta2_success $deltap_success\n");
# }
# fclose($f);
#}
}
# Helper for dump_tv().
function dump_solvency($prefix, $future, $success_prob_range, $success_prob_mode) {
$bucket_array = array();
foreach ($future as $future_age) {
foreach($future_age as $bucket => $_) {
$bucket_array[$bucket] = TRUE;
}
}
$bucket_array = array_keys($bucket_array);
$min_bucket = min($bucket_array);
$max_bucket = max($bucket_array);
# Success/solvency probability percentile lines versus age and wr.
$f = fopen("$prefix-$success_prob_mode.data", "w");
foreach ($future as $age => $future_age) {
fwrite($f, $age);
$prob_range = $success_prob_range;
array_push($prob_range, 2); # Any value greater than 1.
for ($bucket = $max_bucket; $bucket >= $min_bucket; $bucket--) {
$bucket_res = (array_key_exists($bucket, $future_age) ? $future_age[$bucket] : NULL);
if (is_null($bucket_res)) {
continue;
}
$candidate_prob = (($success_prob_mode == 'fp') ? $bucket_res['success_prob'] : $bucket_res['er_success']);
while ($candidate_prob >= $prob_range[0]) {
$min_wr = $bucket_res['bucket'];
fwrite($f, " $min_wr");
array_shift($prob_range);
}
}
array_pop($prob_range);
while (!is_null(array_shift($prob_range))) {
fwrite($f, ' -');
}
fwrite($f, "\n");
}
}
# Dump time varying data files.
function dump_tv($prefix, $future, $gnuplot_bottom_bucket, $gnuplot_top_bucket, $success_prob_range, $equity_range) {
$bucket_array = array();
foreach ($future as $future_age) {
foreach($future_age as $bucket => $_) {
$bucket_array[$bucket] = TRUE;
}
}
$bucket_array = array_keys($bucket_array);
$min_bucket = min($bucket_array);
$max_bucket = max($bucket_array);
# Success/solvency probability fixed probability percentile lines versus age and wr.
dump_solvency($prefix, $future, $success_prob_range, 'fp');
# Success/solvency probability equal risk percentile lines versus age and wr.
dump_solvency($prefix, $future, $success_prob_range, 'er');
# Old style heatmap.
#
## Asset allocation heatmap.
##
#$f = fopen("$prefix-aa.data", "w");
#foreach ($future as $age => $future_age) {
# for ($bucket = $min_bucket; $bucket <= $max_bucket; $bucket++) {
# $bucket_res = (array_key_exists($bucket, $future_age) ? $future_age[$bucket] : NULL);
# $min_wr = $bucket_res['bucket'];
# fwrite($f, $age);
# $prev_res_ef = 0;
# $equity_range_modified = $equity_range;
# array_push($equity_range_modified, 2); # Any value greater than 1.
# foreach ($equity_range_modified as $ef) {
# if (!is_null($bucket_res) && !is_null($bucket_res['wf']) && ($prev_res_ef <= $bucket_res['ef']) && ($bucket_res['ef'] < $ef)) {
# fwrite($f, " $min_wr");
# } else {
# fwrite($f, ' -');
# }
# $prev_res_ef = $ef;
# }
# fwrite($f, "\n");
# }
#}
#fclose($f);
# Asset allocation heatmap.
#
# Gnuplot doesn't support heatmaps with a logarithmic scale, so we have to de-log the data.
#
$f = fopen("$prefix-aa.data", "w");
foreach ($future as $age => $future_age) {
$bucket_res = $future_age[$gnuplot_bottom_bucket];
$prev_display_wr = round($bucket_res['bucket'], 4);
$ef = $bucket_res['ef'];
for ($bucket = $gnuplot_bottom_bucket + 1; $bucket <= $gnuplot_top_bucket; $bucket++) {
$bucket_res = $future_age[$bucket];
$curr_wr = $bucket_res['bucket'];
assert($prev_display_wr < $curr_wr); # Ensure display resolution not set too low; would loose data.
for ($display_wr = $prev_display_wr; $display_wr < $curr_wr; $display_wr += 0.0001) {
$display_wr = round($display_wr, 4);
fwrite($f, "$age $display_wr $ef\n");
}
$ef = $bucket_res['ef'];
$prev_display_wr = $display_wr;
}
fwrite($f, "\n");
}
fclose($f);
}
# Human readable log of time varying data.
function log_tv($future) {
# Cleanup up future by pruning off any extraneous records.
foreach ($future as $age => $future_age) {
$future_keys = array_keys($future_age);
$min_bucket = min($future_keys);
$max_bucket = max($future_keys);
$min_success = NULL;
$max_success = NULL;
foreach ($future_age as $bucket => $record) {
if (!$record['all_success'] && (is_null($min_success) || ($bucket < $min_success))) {
$min_success = $bucket;
}
}
foreach ($future_age as $bucket => $record) {
if (!is_null($record['wf']) && (is_null($min_success) || ($bucket > $max_success))) {
$max_success = $bucket;
}
}
if (is_null($min_success)) {
$min_success = $min_bucket + 1;
}
if (is_null($max_success)) {
$max_success = $max_bucket - 1;
}
foreach ($future_age as $bucket => $record) {
if ((!is_null($max_success) && ($bucket < $min_success - 1)) || ($bucket > $max_success + 1)) {
unset($future[$age][$bucket]);
}
}
}
# Sort.
foreach ($future as $age => $future_age) {
ksort($future[$age]);
}
# Quick at a glance report of equal risk 95th percentile level.
foreach ($future as $age => $future_age) {
foreach ($future_age as $bucket => $record) {
if ($record['er_success'] >= 0.95) {
$wf = $record['wf'];
$ef = $record['ef'];
$er_success = $record['er_success'];
}
}
echo "mwr: $age $wf $ef $er_success\n";
}
print_r($future);
}
function run_tvaa($dump_prefix, $tv_equity_range, $savings, $future_base, $gnuplot_bottom_bucket, $gnuplot_top_bucket,
$le, $returns, $deplete, $original_age, $tv_age, $sex, $years, $years_to_retire, $wrap,
$tv_pctl_success, $tv_min_pctl_success, $tv_min_pctl_stw_success, $success_mode, $success_prob_range, $tv_dump_equity_range, $ef_search_steps) {
$time = time();
$future = precise_tv($tv_equity_range, $savings, $future_base, $gnuplot_bottom_bucket, $gnuplot_top_bucket,
$le, $returns, $deplete, $original_age, $tv_age, $sex, $years, $years_to_retire, $wrap,
$tv_pctl_success, $tv_min_pctl_success, $tv_min_pctl_stw_success, $success_mode, $ef_search_steps);
dump_tv($dump_prefix, $future, $gnuplot_bottom_bucket, $gnuplot_top_bucket, $success_prob_range, $tv_dump_equity_range);
log_tv($future);
$elapsed = time() - $time;
echo "tv done - $elapsed\n";
}
function run_simulate($dump_prefix, $total, $equity_range, $savings, $wr_range,
$le, $returns, $deplete, $age_range, $sex, $years, $years_to_retire, $wrap,
$pctl_success_range, $pctl_stw_success_range, $min_pctl_success, $min_pctl_stw_success) {
$time = time();
$null_variable = NULL; # So we can pass null by reference.
$res = simulate($total, NULL, $equity_range, $savings, $wr_range, $null_variable, NULL, NULL,
$le, $returns, $deplete, NULL, $age_range, $sex, $years, $years_to_retire, $wrap,
NULL, $min_pctl_success, $min_pctl_stw_success, NULL, NULL);
$success = $res['success'];
$stw_success = $res['stw_success']; # Solvency time weighted success.
$elapsed = time() - $time;
echo "successes done - $elapsed\n";
$time = time();
$safe = precise_wr_ef($success, $age_range, $wr_range, $equity_range, $pctl_success_range);
$stw_safe = precise_wr_ef($stw_success, $age_range, $wr_range, $equity_range, $pctl_stw_success_range);
$elapsed = time() - $time;
echo "safes done - $elapsed\n";
$time = time();
dump_data("$dump_prefix-ntw", $success, $safe, $pctl_success_range); # Non-time weighted.
dump_data("$dump_prefix-stw", $stw_success, $stw_safe, $pctl_stw_success_range); # Solvency time weighted.
$elapsed = time() - $time;
echo "dumping done - $elapsed\n";
}
# Main code.
srand(0);
if (PHP_SAPI == "cli") {
$dump_prefix = 'tvaa';
$data_source = 'shiller'; # 'sbbi' or 'shiller'.
$age_start = 0; # Normally 0.
$age_stop = 100; # Normally 100.
$age_step = 1;
$age_range = array();
for ($i = $age_start; $i <= $age_stop; $i += $age_step) {
array_push($age_range, $i);
}
$wr_start = 0.02; # Normally 0.02.
$wr_stop = 0.06; # Normally 0.19.
$wr_step = 0.01; # Fast 0.001, precise 0.0001.
$wr_small_delta = 0.000001;
$wr_range = array();
for ($i = $wr_start; $i <= $wr_stop + $wr_small_delta; $i += $wr_step) {
array_push($wr_range, round($i, 6));
}
$equity_steps = 50; # Fast 50, precise 100.
$equity_range = array();
for ($i = 0; $i <= $equity_steps; $i++) {
array_push($equity_range, round($i / $equity_steps, 3));
}
$pctl_success_range = array(0.95);
$pctl_stw_success_range = array(0.99);
$min_pctl_success = 0.9;
$min_pctl_stw_success = 0.95;
# Cutoff below which we don't generate success/stw_success results for that age or younger.
# Use caution when setting these value. They can yield a major performance win.
# But it can result in success/stw_success values greater or equal this value not getting generated if death array is not monotone increasing.
# Set to zero for safety.
# Time varying simulation specific parameters.
$allow_interpolate = FALSE; # Enable if safe to assume ef constant between two identical ef values as a function of wr.
# In practice ef may have a local minimum or maximum, so the results from doing this will not be precise. But provides a 60% speedup.
$future_base = pow(1.02, 1); # Future buckets spaced in min_wf space are a factor of this much appart.
# Higher resolution increases run time, but provides more opportunities for and a higher success probability.
# Switching resolutions might not result in continuous improvements unless the bucket grids are factors of one another.
# And even then possibly not for reasons not understood.
$gnuplot_bottom_bucket = intval(floor(log(0.02) / log($future_base))); # Fill in and report data all the way down to the gnuplot floor.
$gnuplot_top_bucket = intval(ceil(log(0.10) / log($future_base))); # Fill in and report data all the way up to the gnuplot ceiling.
$tv_equity_steps = 4; # Use 4 steps to mirror 5 band Trinity study. 20 for expanded asset allocation choices.
$tv_equity_range = array();
for ($i = 0; $i <= $tv_equity_steps; $i++) {
array_push($tv_equity_range, round($i / $tv_equity_steps, 3));
}
$tv_dump_equity_range = array(0.125, 0.375, 0.625, 0.875); # Unused. Was used with old-style heatmaps. Inter-equity resolution points with which to display.
$success_prob_range = array(0.9, 0.95, 0.98, 0.99, 1); # Used when dumping data to produce percentile success graphs.
$original_age = 50; # Determines risk fraction base on uniform risk allocation of risk used for equal risk percentile lines.
$tv_age = 50; # Analyze from this age on.
$tv_pctl_success = 0; # Probabilities below this we don't consider.
# Should consider probabilites below the end probability we are interested in because there may only be a small chance of improbable paths being taken.
$tv_min_pctl_success = 0; # Has no impact on execution time; setting to non-zero would prevent all results from being considered.
$tv_min_pctl_stw_success = 0;
$success_mode = 'stw'; # 'stw' for solvency time weighted, or 'ntw' for non-time weighted.
$ef_search_steps = NULL; # Number of tv_equity_range steps either side of equity range value that must be less to determine when we have found the maximum.
# NULL to search all.
# Due to local maxima, will produce sub-optimal results unless search all equity range values.
$total = 1;
$savings = 0;
$ret_equity = NULL; # NULL for no adjustment, float for target geomean return.
$ret_bonds = NULL; # NULL for nominals, float for TIPS.
$sim_start_year = 1872;
$sim_end_year = 2010; # NULL for until end of data.
$deplete = 'rb';
$sex = 'combined'; # 'combined' for avareage male/female death probabilities, 'die_exp' for living forever, or whatever else that death map may define.
$years = NULL; # Years to run. Set to NULL to use the full death array length.
$years_to_retire = 0;
$wrap = TRUE; # Allow time periods to wrap.
if ($data_source == 'sbbi') {
load_sbbi_data($sbbi_file);
$initial_year = $sbbi_initial;
} else {
load_shiller_data($shiller_file);
$initial_year = $shiller_initial;
}
$le = pre_compute_life_expectancy($death[$sex]);
$returns = massage_returns($ret_equity, $ret_bonds, 0, $initial_year, $sim_start_year, $sim_end_year);
#$gm = plus_1_geomean($stock) - 1;
#$sd = standard_deviation($stock);
#echo "stock $gm +/- $sd\n";
#
#$gm = plus_1_geomean($bond) - 1;
#$sd = standard_deviation($bond);
#echo "bond $gm +/- $sd\n";
if (TRUE) {
# Time Varying Asset Allocator.
run_tvaa($dump_prefix, $tv_equity_range, $savings, $future_base, $gnuplot_bottom_bucket, $gnuplot_top_bucket,
$le, $returns, $deplete, $original_age, $tv_age, $sex, $years, $years_to_retire, $wrap,
$tv_pctl_success, $tv_min_pctl_success, $tv_min_pctl_stw_success, $success_mode, $success_prob_range, $tv_dump_equity_range, $ef_search_steps);
} else {
# Batch form HARS simulator.
run_simulate($dump_prefix, $total, $equity_range, $savings, $wr_range,
$le, $returns, $deplete, $age_range, $sex, $years, $years_to_retire, $wrap,
$pctl_success_range, $pctl_stw_success_range, $min_pctl_success, $min_pctl_stw_success);
}
} else {
echo "\n";
echo "HARS - Historical Actuarial Retirement Simulator\n";
echo "";
echo "HARS - Historical Actuarial Retirement Simulator
\n";
echo "\n";
if (isset($_GET['total'])) {
# Run the simulator.
load_shiller_data($shiller_file);
# Get simulation parameters.
$total = $_GET['total'] + 0;
$equity_steps = 10; # Simulate different bond/stock ratios $total/$equity_steps appart
$limit = $_GET['limit'] + 0;
$savings = $_GET['savings'] + 0;
$income = $_GET['income'] + 0;
$ret_equity = ($_GET['ret_equity'] + 0) / 100;
if ($_GET['ret_bonds'] == '') {
$ret_bonds = NULL;
} else {
$ret_bonds = ($_GET['ret_bonds'] + 0) / 100;
}
$sim_start_year = intval($_GET['sim_start_year']);
(($sim_start_year >= $shiller_initial) && ($sim_start_year <= $shiller_initial + count($stock) - 1)) or die("Invalid simulation start year.");
# Minus 2 rather than minus 1 so we have at least two samples and can compute the standard deviation of returns later on.
$seq = (($_GET['seq'] == 'seq') ? 'seq' : 'rnd');
$deplete = $_GET['deplete'];
preg_match('/^(ef|bf|rb)$/', $deplete) or die("Invalid depletion mode.");
$num_experiments = $_GET['num_experiments'] + 0;
($num_experiments == 1000) || ($num_experiments == 3000) or die("Invalid number of experiments.");
$sex1 = $_GET['sex1'];
preg_match('/^(combined|male|female)$/', $sex1) or die("Invalid first person gender");
$years1 = count($death[$sex1]);
$age1 = intval($_GET['age1']); # Set to -1 to never die.
(($age1 >= -1) && ($age1 < $years1)) or die("Invalid first person age.");
$sex2 = $_GET['sex2'];
preg_match('/^(combined|male|female)$/', $sex2) or die("Invalid second person gender");
$years2 = count($death[$sex2]);
$age2 = intval($_GET['age2']); # Set to -1 if effectively already dead.
(($age2 >= -1) && ($age2 < $years2)) or die("Invalid second person age.");
$years = min($years1, $years2);
$years_to_retire = intval($_GET['years_to_retire']);
($years_to_retire >= 0) or die("Invalid years until retirement.");
$go_to_bonds = $_GET['go_to_bonds'] + 0;
($go_to_bonds == -1) || (($go_to_bonds >= 0) && ($go_to_bonds <= 100)) or die("Invalid go to bonds percentage");
$income_dead_fraction = $_GET['income_dead_fraction'] / 100;
($income_dead_fraction >= 0) or die("Invalid income dead proportion.");
$diversification = $_GET['diversification'] / 100;
(($diversification >= 0) && ($diversification <= 1)) or die("Invalid diversification.");
$correlation = $_GET['correlation'] / 100;
(($correlation >= -1) && ($correlation <= 1)) or die("Invalid correlation.");
$diversify_sds = $_GET['diversify_sds'] + 0;
($diversify_sds >= 0) or die("Invalid diversification standard deviation scaling factor.");
echo "\n";
# Output input parameters.
echo "
\n";
printf("initial assets %4dk max assets %4dk return on equity %3.1f%%", $total / 1000, $limit / 1000, $ret_equity * 100);
if (is_null($ret_bonds)) {
printf(" nominal bonds\n");
} else {
printf(" return on bonds %3.1f%%\n", $ret_bonds * 100);
}
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($stock) - 1, $num_experiments);
# Construct array of equity values to use
$equity_range = array();
for ($i = 0; $i <= $equity_steps; $i++) {
array_push($equity_range, $i / $equity_steps);
}
if (($age1 != -1) && ($deplete == 'rb')) {
array_push($equity_range, 'age-in-bonds');
}
# # Report statistics
#
# $samp_prod = 1;
# for ($i = 0; $i <= count($equity_returns) - 1; $i++) {
# $samp_prod *= 1 + $equity_returns[$i];
# }
# $samp_geomean = pow($samp_prod, 1 / count($equity_returns)) - 1;
#
# $samp_tot = 0;
# for ($i = 0; $i <= count($equity_returns) - 1; $i++) {
# $samp_tot += $equity_returns[$i];
# }
# $samp_mean = $samp_tot / count($equity_returns);
#
# $samp_ss = 0;
# for ($i = 0; $i <= count($equity_returns) - 1; $i++) {
# $samp_ss += pow($equity_returns[$i] - $samp_mean, 2);
# }
# $samp_sd = pow(1 / (count($equity_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);
#
# echo "\n";
echo "
\n";
echo "
\n";
echo "
| \n";
echo " |
|---|
| Equity | Bonds | Success probability \n";
echo " | Mean gifted assets | Mean assets at final death \n";
# Simulate portfolios
foreach ($equity_range as $ef) {
$res = simulate_web($total, $ef, $limit, $savings, $income, $ret_equity, $ret_bonds, $pe_adjust, $sim_start_year, $seq, $deplete, $num_experiments,
$age1, $sex1, $age2, $sex2, $years, $years_to_retire, $go_to_bonds, $income_dead_fraction, $diversification, $correlation, $diversify_sds);
$success = $res['success'];
$p_array = $res['p_array'];
$charitable_array = $res['charitable_array'];
# Compute and report stats
# $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)];
$tot = array_sum($p_array);
$charitable_tot = array_sum($charitable_array);
$mean = $tot / $num_experiments;
$charitable_mean = $charitable_tot / $num_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($success * 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 "\n";
echo "\n";
}
?>