Mathematical Foundations: Fire Planner Engine

This article derives and proves every formula in the Fire Planner projection engine. For a quicker overview of what each feature does and its default values, see How Projections Work.

The engine runs a month-by-month simulation. Each section below corresponds to one piece of that simulation, with the math shown from first principles.

1. Compound Growth Derivation

The future value formula is the backbone of the projection engine. Here’s where it comes from.

Setup

Let P be the initial principal, C the fixed monthly contribution, and r the monthly growth rate. The engine applies growth first, then adds the contribution (ordinary annuity convention).

After month 1:

B1 = P(1 + r) + C

After month 2:

B2 = B1(1 + r) + C = P(1 + r)2 + C(1 + r) + C

After month n:

Bn = P(1 + r)n + C[(1 + r)n-1 + (1 + r)n-2 + … + (1 + r) + 1]

Geometric Series Collapse

The bracketed term is a geometric series with ratio (1 + r) and n terms. The sum of a geometric series a + ar + ar² + … + ar^n-1 is a(rⁿ - 1) / (r - 1). Here a = 1 and the ratio is (1 + r):

1 + (1 + r) + (1 + r)2 + … + (1 + r)n-1 = ((1 + r)n - 1) / r

Substituting back:

FV = P(1 + r)n + C × ((1 + r)n - 1) / r

The first term is compound growth on the initial balance. The second is the future value of the contribution stream. The engine doesn’t use this closed-form directly. It runs the month-by-month loop because contributions, growth rates, and expenses change over time. But the closed form validates the loop’s output for constant-input scenarios, as long as income funds the contributions (unfunded contributions are paid from savings, so they don’t add to the total).

End-of-Month Convention

Contributions are added after growth is applied. This is the ordinary annuity convention: money deposited at month-end earns no interest in the deposit month. The alternative (annuity due) would apply growth to the contribution in the same month, producing slightly higher balances. The ordinary convention is standard for modeling paycheck-to-investment flows.

2. Monthly Rate Conversion

The engine converts annual rates to monthly by simple division:

rmonthly = rannual / 12

This is the nominal rate convention. A 7% annual rate becomes 0.5833% monthly. Over one year:

(1 + 0.07/12)12 = 1.0723

The effective annual rate is 7.23%, not 7.00%. The alternative, computing (1 + rannual)1/12 - 1, would produce exactly 7.00% effective annual. The difference over 30 years at 7%:

• Nominal convention: $100,000 → $811,651
• Exact 12th root: $100,000 → $761,226

The ~6.6% difference is noticeable but both conventions are used across financial planning tools. The key is consistency: every rate in the system (asset growth, liability interest, income growth, expense inflation) uses the same nominal convention. When you enter 7%, the engine treats it as 7% nominal compounding monthly.

3. Income & Expense Growth

Growth Model

Both income and expenses use the same formula:

Valuet = Value0 × (1 + g/100)t

Where t is yearsElapsed, calculated as:

yearsElapsed = yearOffset + (month - startMonth) / 12

This produces fractional years, giving smooth month-to-month growth rather than step-function jumps at year boundaries. A salary of $80,000/year with 3% growth in month 6 of year 5:

yearsElapsed = 5 + 6/12 = 5.5
salary = 80,000 × 1.035.5 = $93,790

Per-Category Inflation

Each expense category has its own inflation rate. This matters enormously over long horizons because of the compounding divergence:

Starting from $500/month, over 25 years:

• Healthcare at 5%: 500 × 1.0525 = $1,693
• Housing at 3%: 500 × 1.0325 = $1,047
• Fun at 2.5%: 500 × 1.02525 = $928

Healthcare triples. Housing doubles. Fun barely doubles. A flat 3% inflation rate for all categories would understate healthcare costs by 38% and overstate fun by 13% at year 25. The per-category model captures these dynamics.

4. Liability Amortization

Monthly Update

Each month, interest accrues first, then the payment is applied:

Balancenew = Balance × (1 + rannual/100/12) - Payment

The payment is capped at the remaining balance:

Paymenteffective = min(Payment, Balance + Interest)

This prevents negative balances. When the final payment arrives, it’s only the remaining amount, not the full scheduled payment.

The Negative Amortization Case

When Payment < Balance × r/12, the monthly payment doesn’t cover interest. The balance grows every month:

Example: $100,000 loan at 6%, $200/month payment.

• Monthly interest: $100,000 × 0.005 = $500
• Payment: $200
• Net change: +$300/month → balance grows

The engine handles this correctly (the balance increases indefinitely) and the warnings system flags it so you can adjust.

5. Retirement Withdrawal Math

Net Expense Calculation

The amount to withdraw each month:

NetExpenses = MonthlyExpenses + DebtPayments - TrackedIncome

If net expenses are negative (income exceeds outflows even in retirement), no withdrawal is needed and the surplus is deposited.

Proportional Withdrawal: Proof of Correctness

The proportional strategy withdraws from all eligible accounts simultaneously, weighted by balance:

Withdrawali = min(Balancei, NetExpenses × Balancei / TotalBalance)

Claim: When TotalBalance ≥ NetExpenses, the sum of all withdrawals equals NetExpenses exactly, and no individual withdrawal exceeds its account balance.

Proof: When TotalBalance ≥ NetExpenses, the proportion for account i is Balancei / TotalBalance. The withdrawal is NetExpenses × Balancei / TotalBalance. Since NetExpenses ≤ TotalBalance, we have:

Withdrawali = NetExpenses × Balancei / TotalBalance ≤ Balancei

So the min() cap never activates. Summing across all accounts:

∑ Withdrawali = NetExpenses × ∑(Balancei / TotalBalance) = NetExpenses × 1 = NetExpenses

The total withdrawal is exactly what’s needed. □

When TotalBalance < NetExpenses, each account is drained to zero and there’s a shortfall. The engine records this as a partial depletion.

Taxable First: Sequential Depletion

Assets are sorted into priority tiers:

1. Taxable brokerage, high-yield savings, cash: no penalty at any age
2. Retirement accounts (401k, IRA) and HSA: tax-deferred
3. Roth (Roth IRA, Roth 401k): tax-free growth
4. Everything else (real estate equity, crypto, other)

The engine drains each account completely before moving to the next. This ordering is tax-optimal under common assumptions: taxable accounts have the lowest tax benefit (already taxed), so spend them first. Roth accounts have the highest benefit (tax-free forever), so preserve them longest.

The 59.5 Rule

When enabled, retirement accounts (types: retirement, roth, hsa) are excluded from the eligible withdrawal pool before age 59.5. The age is calculated from the birth date:

currentAge = (currentDate - birthDate) / 365.25

If all non-retirement accounts are depleted and the retiree still needs money, retirement accounts are accessed anyway, but flagged as early withdrawals. In the real world, this incurs a 10% penalty plus income tax. The engine doesn’t model the penalty amount but the flag tells you it would apply.

6. Monte Carlo: Lognormal Returns

This section derives the most mathematically dense formula in the engine: drift-corrected lognormal return generation.

Why Lognormal?

Stock returns are multiplicative, not additive. If your portfolio drops 50%, you need a 100% gain to recover, not 50%. This asymmetry means returns can’t be normally distributed (normal distributions are symmetric and allow returns below -100%). The solution: model the logarithm of the growth factor as normally distributed.

If log(1 + R) is normally distributed, then 1 + R is lognormally distributed, which guarantees R > -1 (you can never lose more than 100%).

The Drift Correction Derivation

This is the critical piece. We need monthly returns such that their compound product over 12 months has the correct expected value.

Goal: E[∏i=112(1 + Ri)] = 1 + rannual

Step 1. Let log(1 + R) ~ N(μ, σ²). Then 1 + R is lognormal, and its expectation is:

E[1 + R] = exp(μ + σ²/2)

This is a standard result for lognormal random variables. The σ²/2 term is the lognormal bias. It pushes the mean above the median.

Step 2. For 12 independent monthly returns, the expected product is:

E[∏(1 + Ri)] = [E[1 + R]]12 = [exp(μ + σ²/2)]12 = exp(12μ + 6σ²)

Step 3. Set equal to the target annual return:

exp(12μ + 6σ²) = 1 + rannual

Take the natural log:

12μ + 6σ² = ln(1 + rannual)

Solve for μ:

μ = ln(1 + rannual)/12 - σ²/2

The Correction’s Impact

Without the -σ²/2 correction, the simulation is biased upward. At 15% annual volatility (the S&P 500 historical average):

σmonthly = 0.15 / √12 = 0.0433
σ²/2 = 0.0433² / 2 = 0.000938

Annualized, this bias is approximately:

exp(12 × 0.000938) - 1 = 1.13%

A simulation configured for 7% growth would actually produce ~8.1% average returns without the correction. Over 30 years, that turns $100,000 into $1,006,266 instead of $811,651, a 24% overstatement. The correction eliminates this bias.

Monthly Volatility

Annual volatility converts to monthly via:

σmonthly = σannual / √12

This assumes returns are independent and identically distributed (iid) across months. Under iid, variances are additive: σ²annual = 12 × σ²monthly, which gives the √12 scaling. Real markets exhibit some serial correlation (momentum/mean-reversion), but the iid assumption is standard for retirement planning horizons.

Normal Random Number Generation

The engine uses the Marsaglia polar method (a variant of Box-Muller) to convert uniform random numbers to standard normal draws:

1. Generate u, v uniformly in [-1, 1]
2. Compute s = u² + v²
3. If s ≥ 1 or s = 0, reject and retry
4. Compute mul = √(-2 ln(s) / s)
5. Return z₁ = u × mul and cache z₂ = v × mul for the next call

This produces two independent N(0,1) draws per iteration. The polar method avoids the trigonometric functions in the original Box-Muller, trading them for a rejection step (which succeeds ~78.5% of the time, since the rejection region is a unit circle inscribed in a 2×2 square).

Mulberry32 PRNG

The engine uses Mulberry32, a seedable 32-bit pseudorandom number generator. Seedability is important: given the same seed, the simulation produces identical results. This enables reproducible testing and lets users compare scenarios without random noise confounding the comparison.

Every plan runs on the same seed (common random numbers): simulation i sees the same market in month k whatever the plan. So the same plan gets the same count on every device and every load, two plans with the same money get the same count, and a plan that costs more never scores better by luck. Simulations come in pairs, the second with the first's draws negated (antithetic variates), so each month's draws average to zero and the one seed is neither a lucky nor an unlucky market.

Putting It Together

For each Monte Carlo run and each month:

1. Draw z ~ N(0,1) using Marsaglia polar method
2. Compute monthly return: R = exp(μ + σm × z) - 1
3. Apply to market-rate assets: Balance × (1 + R)
4. Custom-rate assets keep their fixed monthly rate

After all runs complete, the engine sorts final portfolio values at each year and extracts percentiles using linear interpolation:

Valuep = Values[i] × (1 - f) + Values[i+1] × f

Where i = floor((p/100) × (N-1)) and f = (p/100) × (N-1) - i. This produces P10, P25, P50 (median), P75, and P90 bands.

The headline count is the number of runs where assets never hit $0 before your end year: “money lasts in 462 of 500 simulations”, a count of real runs and never a percentage.

7. Historical Backtesting

Data Source

The backtester uses Robert Shiller’s S&P 500 monthly dataset, spanning 1871 to present (~1,850 monthly observations). Each record contains the year, month, and nominal monthly return (price change + reinvested dividends).

Why Nominal Returns

This is a critical design decision. The projection engine already applies per-category inflation rates to expenses (healthcare at 5%, housing at 3%, etc.). If the backtester used CPI-adjusted (real) returns, inflation would be subtracted twice:

1. Once by the CPI adjustment to returns (lowering portfolio growth)
2. Once by expense inflation (increasing outflows)

Using nominal returns avoids this double-counting. The portfolio grows at nominal rates, expenses inflate at category-specific rates, and the gap between them correctly represents real purchasing power dynamics.

Period Selection

The engine tests every January-starting period with enough historical data to cover the full projection length. For a 50-year projection (600 months), the backtester needs at least 600 monthly records from each start point. Starting from 1872 through roughly the mid-1970s, this gives ~100+ overlapping test periods.

January-only starts are preferred for interpretability (“what if you started in 1929?”). If no January starts are available (shouldn’t happen with Shiller data), the engine falls back to every available month.

Asset-Level Routing

Only “market-rate” (risky) assets receive historical returns. The engine checks two conditions:

1. The asset uses the market rate (not a custom override)
2. The asset is a “risky” type (not cash or savings)

Assets with custom growth rates (e.g., a rental property at 4%) keep their fixed rate regardless of what the S&P 500 did in 1973. This models the reality that not all assets are correlated with equities.

Results

Same percentile extraction as Monte Carlo (P10/P25/P50/P75/P90 at each year). The backtester additionally identifies:

• Worst-case starting year: the period that produced the lowest final net worth
• Best-case starting year: the period that produced the highest final net worth
• "Money lasts in N of M start years": the count of start years where assets never hit $0 before your end year

8. Survival Heatmap

The heatmap uses a simplified Trinity Study model, deliberately separate from the full projection engine.

The Model

For each combination of withdrawal rate and historical starting year:

1. Start with a portfolio of 1 in 100% stocks: the S&P 500 with dividends reinvested (survival is scale-free, so the table shows rates, never amounts)
2. For each year, compute the annual return from 12 monthly Shiller records
3. Grow the portfolio: portfolio × yearlyMultiplier
4. At year end, after that return, subtract the inflation-adjusted withdrawal: initialWithdrawal × CPIJan of this year / CPIJan of the start year
5. If portfolio ≤ 0, mark as failed

Inflation-Adjusted Withdrawals

Withdrawals rise with actual inflation, using Shiller’s monthly CPI series, as in the original Cooley, Hubbard, and Walters (Trinity) study: a retiree who starts at 4% raises the dollar withdrawal each year to keep its purchasing power. The first year’s withdrawal is the starting amount; each later year’s is scaled by the ratio of January CPI to January CPI in the start year.

Withdrawalyear = InitialWithdrawal × CPIJan(start + year) / CPIJan(start)

This is why high-inflation starts hurt. A retiree starting in January 1966 at 5% saw prices roughly triple over the next 16 years while stocks, dividends included, lost about 15% after inflation. The money ran out in 1985, year 20 of 30. A flat 3% raise would have hidden that.

Over 30 years on the current data (125 starting years), 4% lasted from every starting year, 5% from 88.8% of them and 6% from 74.4% of them.

The retirement length comes from your plan: from your retirement date (today if you have already retired, your FI date if there is no retirement date) to your end year, or 30 years if the plan has no date.

Why Simplified?

The full projection engine handles income, expenses, multiple assets, drawdown strategy, contributions, stress tests, and more. The heatmap answers a simpler, more universal question: “Given just a portfolio and a withdrawal rate, what does history say?”

This is the same question the original Trinity Study answered, making the heatmap directly comparable to published 4% rule research.

9. Stress Test Models

Stress tests are either permanent transforms (applied once before simulation) or year-based modifiers (applied at specific years during the loop).

Market Crash

Type: Year-based modifier

Applies an instant drop to risky asset balances in the trigger year, followed by a linear growth penalty over the years after it:

Crash year: balance = balance × (1 - dropPercent/100)

Growth penalty = -3% × (1 - yearsAfterCrash / slowerGrowthYears)

Safe assets (cash, savings) are protected. With 5 years of slower growth the penalty is -2.4% in the first year after the crash, then -1.8%, -1.2%, -0.6% and 0%. Growth returns to normal, but the drop itself is never made back: balances stay below where they would have been.

Lower Returns

Type: Permanent transform

Subtracts a fixed percentage from every growth rate before simulation begins:

ratenew = max(0, rateoriginal - reduceBy)

The floor at 0% prevents negative growth rates in the deterministic projection (Monte Carlo can still produce negative individual months).

Income Loss

Type: Year-based modifier

Binary shutoff: if yearOffset ≥ triggerYear: income = 0

All income streams are zeroed simultaneously. This is a worst-case model, not a “one income stream drops” scenario.

Spending Spike

Type: Year-based modifier

Permanent multiplier starting at the trigger year:

expenses = expenses × (1 + increasePercent/100)

High Inflation (Cumulative Compounding)

Type: Year-based modifier (time-limited rate, permanent price level)

During the inflation window, expenses are multiplied by a cumulative factor that grows each year:

multiplier = ((1 + R/100) / (1 + D/100))yearsInSpike

Where R is the spike inflation rate, D is the default inflation rate, and yearsInSpike is 1-based (first year of spike = 1). It stays at the spike's length once the spike ends.

Why the ratio? The per-expense inflation rates (healthcare at 5%, housing at 3%, etc.) already model baseline inflation. The stress test needs to add extra inflation above baseline. The ratio captures the gap: if default is 3% and the spike is 8%, the extra factor per year is 1.08/1.03 = 1.0485, approximately 4.85% additional inflation above what’s already modeled.

Why cumulative? Inflation compounds. If prices rise 4.85% extra in year 1 and 4.85% extra in year 2, the total gap after 2 years isn’t 9.7%. It’s 1.0485² = 1.099, or 9.9%. After 5 years: 1.0485⁵ = 1.267, or 26.7% above baseline.

This matches the theoretical result: (1.085) / (1.035) = 1.469 / 1.159 = 1.267.

After the window ends, the multiplier stays at its last value (1.267 above). Per-expense rates resume from that level, so prices keep what the spike added.

Retire Earlier or Later

Type: Permanent transform

The retirement date moves by whole years: retirement = retirement + years, with years from -20 to +30 (default -2). A date still ahead never moves to before this month. Income that starts on the retirement date moves with it. On a suggested date the shift applies after the solver has found the date, since the solver replaces the retirement date with every month it tries.

10. Edge Cases & Safety

Zero Growth Rate

When r = 0, the compound growth formula degenerates to linear accumulation:

FV = P + C × n

The engine handles this naturally: Balance × (1 + 0) = Balance, and contributions simply stack.

Negative Returns in Monte Carlo

The lognormal distribution ensures 1 + R > 0, so the monthly return is always greater than -100%. However, returns can be severely negative (a 3-sigma event at 15% volatility gives a monthly return of about -12%). The engine applies these directly, with no artificial floor.

Asset Balance Floor

Display uses max(0, balance) to prevent negative values in charts and tables. Internally, balances can briefly go negative during the monthly calculation (e.g., withdrawal exceeds balance before next asset in priority kicks in), but the display floor prevents confusing visual artifacts.

Depletion Detection

In the lightweight projection (used by Monte Carlo and backtesting), depletion is detected when all asset balances are ≤ 0 during the retirement phase, or in a working month whose shortfall (outflows above income) needed them:

if ((isRetirementPhase || shortfall > 0) && totalAssets ≤ 0) → depleted

The simulation records the depletion month and marks survived = false. The main projection continues running (to show what happens after depletion), and the run does not count among the runs where the money lasts.

Data Boundary Safety

The backtester checks array bounds before accessing historical records. If the data runs out mid-simulation (shouldn’t happen with proper period selection, but defensive code), it falls back to the asset’s configured growth rate.