Mathematical Foundations: Fire Planner Engine

Every formula derived, every assumption justified

firenum team · Last updated February 10, 2026

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² + … + arn-1 is a(rn - 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.

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%:

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 triples. Housing doubles. Discretionary barely doubles. A flat 3% inflation rate for all categories would understate healthcare costs by 38% and overstate discretionary costs 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.

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.

Tax-Efficient Sequential Depletion

Assets are sorted into priority tiers:

  1. Taxable (brokerage, savings, cash) — no penalty at any age
  2. Traditional retirement (401k, IRA, HSA) — tax-deferred
  3. Roth (Roth IRA, Roth 401k) — tax-free growth
  4. 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 z1 = u × mul and cache z2 = 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.

The default seed is the current timestamp, so each “Run Simulation” produces a different sequence. Fix the seed for reproducibility.

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.

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 1871 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:

8. Survival Heatmap

The heatmap uses a simplified Trinity Study model—deliberately not the full projection engine.

The Model

For each combination of withdrawal rate and historical starting year:

  1. Start with a fixed portfolio (default: $1,000,000)
  2. For each year, compute the annual return from 12 monthly Shiller records
  3. Grow the portfolio: portfolio × yearlyMultiplier
  4. Subtract the inflation-adjusted withdrawal: initialWithdrawal × (1.03)year
  5. If portfolio ≤ 0, mark as failed

Inflation-Adjusted Withdrawals

Withdrawals grow at 3% annually, matching the Trinity Study methodology. The original Cooley, Hubbard, and Walters study tested inflation-adjusted withdrawals: a retiree who starts at 4% increases their dollar withdrawal each year to maintain purchasing power.

Withdrawalyear = InitialWithdrawal × (1.03)year

The 3% inflation rate matches the project’s default inflation rate. With a $1M portfolio at 4% withdrawal rate:

By year 30, you’re withdrawing nearly 10% of the original portfolio annually. This is why survival rates at higher withdrawal rates drop significantly—the compounding withdrawal catches up to portfolio growth.

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 during recovery:

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

Recovery penalty = -3% × (1 - yearsIntoRecovery / recoveryYears)

Safe assets (cash, savings) are protected. The recovery penalty ramps from -3% in the first recovery year back to 0% at the end, modeling the gradual return to normal market conditions.

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.

Expense 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)

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).

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.04855 = 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 drops back to 1.0 and per-expense rates resume as normal.

Retire Earlier/Later

Type: Permanent transform

Shifts the retirement date: retirementDate.setFullYear(retirementDate.getFullYear() + years)

This changes when the engine switches from accumulation (contributions + surplus deposits) to withdrawal mode. Moving retirement earlier shortens the saving period and lengthens the drawdown period—a double hit.

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—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:

if (isRetirementPhase && totalAssets ≤ 0) → depleted

The simulation records the depletion month and marks survived = false. The main projection continues running (to show what happens after depletion), but the success rate counts it as a failure.

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.

FAQ

Why use lognormal distributions for Monte Carlo instead of normal distributions?

Stock returns are multiplicative, not additive. A -50% loss followed by a +50% gain leaves you at 75%, not 100%. Normal distributions allow impossible returns below -100%. Lognormal distributions model this multiplicative behavior correctly—the logarithm of the growth factor is normally distributed, ensuring returns never drop below -100%.

What is the drift correction and why does it matter?

The drift correction (-σ²/2) compensates for the mathematical property that E[exp(X)] > exp(E[X]) (Jensen’s inequality). Without it, a Monte Carlo simulation with 7% expected return and 15% volatility would produce an average return of about 8.1%—systematically too optimistic. The correction ensures the expected compound return matches your configured growth rate.

Why does the backtester use nominal returns instead of inflation-adjusted returns?

The projection engine already applies inflation to each expense category individually (healthcare at 5%, housing at 3%, etc.). If the backtester also used CPI-adjusted returns, inflation would be subtracted twice—once from the returns and once from the expense growth. Nominal returns avoid this double-counting.

How accurate is the monthly compounding approximation?

The engine divides annual rates by 12 for monthly compounding (nominal rate convention). At 7% annual, this produces an effective annual rate of 7.23%. The alternative—computing the exact 12th root—would give exactly 7.00% effective annual rate. The difference is small and the nominal convention is standard in financial planning tools. The key is consistency: all rates in the system use the same convention.

Can the proportional withdrawal strategy over-withdraw from any account?

No. The withdrawal from each account is capped at that account’s balance via min(balance, netExpenses × proportion). When the total eligible balance exceeds net expenses, the proportions sum to exactly 1 and each withdrawal is less than or equal to the account balance. When total balance is less than net expenses, each account is drained to zero.

→ Open the Fire Planner
→ How Projections Work (quicker overview)
→ Practical guide: setting up and using the tracker
→ The 4% rule explained