zfin/src/analytics/projections.zig

/// Historical simulation engine for retirement projections.
///
/// Implements the FIRECalc algorithm: for each starting year in the Shiller
/// historical dataset (1871–present), simulate a retirement of `horizon` years
/// using actual market returns, bond returns, and inflation. The portfolio is
/// rebalanced annually to the target stock/bond allocation.
///
/// Key outputs:
///   - Safe withdrawal amount at a given confidence level (binary search to $1)
///   - Success rate for a given spending level
///   - Percentile bands of portfolio value at each year (for charting)
const std = @import("std");
const log = std.log.scoped(.projections);
const shiller = @import("../data/shiller.zig");
const srf = @import("srf");
const Date = @import("../models/date.zig").Date;

// ── Life events ─────────────────────────────────────────────────

/// A resolved event ready for the simulation loop. All age-based timing
/// has been converted to simulation years. The simulation functions only
/// need this — no person indices, no ages array.
pub const ResolvedEvent = struct {
    start_year: u16,
    duration: u16, // 0 = permanent
    annual_amount: f64, // positive = income, negative = expense
    inflation_adjusted: bool,

    pub fn isActive(self: *const ResolvedEvent, y: u16) bool {
        if (y < self.start_year) return false;
        if (self.duration == 0) return true;
        return y < self.start_year + self.duration;
    }

    pub fn cashFlow(self: *const ResolvedEvent, y: u16, cumulative_inflation: f64) f64 {
        if (!self.isActive(y)) return 0;
        if (self.inflation_adjusted) return self.annual_amount * cumulative_inflation;
        return self.annual_amount;
    }
};

/// A discrete cash flow event that modifies the simulation's annual
/// withdrawal. Positive amount = income (reduces withdrawal, e.g.
/// Social Security). Negative = expense (increases withdrawal, e.g.
/// college tuition).
pub const LifeEvent = struct {
    name: [max_name_len]u8 = .{0} ** max_name_len,
    name_len: u8 = 0,
    start_age: u16,
    person: u8 = 0, // 0-indexed into birthdates array
    duration: u16 = 0, // 0 = permanent (until end of horizon)
    annual_amount: f64, // positive = income, negative = expense
    inflation_adjusted: bool = true,

    const max_name_len = 48;

    pub fn getName(self: *const LifeEvent) []const u8 {
        return self.name[0..self.name_len];
    }

    /// Simulation year when this event starts, given the persons' current ages.
    /// Returns null if the person index is out of range.
    pub fn startYear(self: *const LifeEvent, current_ages: []const u16) ?u16 {
        if (self.person >= current_ages.len) return null;
        const age = current_ages[self.person];
        if (self.start_age <= age) return 0;
        return self.start_age - age;
    }

    /// Is this event active in simulation year `y`?
    pub fn isActive(self: *const LifeEvent, y: u16, current_ages: []const u16) bool {
        const start = self.startYear(current_ages) orelse return false;
        if (y < start) return false;
        if (self.duration == 0) return true; // permanent
        return y < start + self.duration;
    }

    /// Cash flow contribution for simulation year `y`.
    /// Positive = income (reduces net withdrawal), negative = expense.
    pub fn cashFlow(self: *const LifeEvent, y: u16, cumulative_inflation: f64, current_ages: []const u16) f64 {
        if (!self.isActive(y, current_ages)) return 0;
        if (self.inflation_adjusted) return self.annual_amount * cumulative_inflation;
        return self.annual_amount;
    }

    /// Resolve this event into a ResolvedEvent using the given current ages.
    /// Returns null if the person index is out of range.
    pub fn resolve(self: *const LifeEvent, current_ages: []const u16) ?ResolvedEvent {
        const start = self.startYear(current_ages) orelse return null;
        return .{
            .start_year = start,
            .duration = self.duration,
            .annual_amount = self.annual_amount,
            .inflation_adjusted = self.inflation_adjusted,
        };
    }
};

// ── User configuration (from projections.srf) ──────────────────

/// User-configurable projection parameters, loaded from projections.srf.
///
/// Example projections.srf (union-tagged SRF records):
///   #!srfv1
///   type::config,target_stock_pct:num:77
///   type::config,horizon:num:30
///   type::config,horizon_age:num:90           # resolves to (90 − oldest current age)
///   type::birthdate,date::1975-03-15
///   type::birthdate,date::1978-06-22,person:num:2
///   type::event,name::Social Security,start_age:num:67,amount:num:38400
pub const UserConfig = struct {
    /// Target stock allocation percentage (0-100). Used for simulation blending.
    target_stock_pct: ?f64 = null,
    /// Retirement horizons to simulate (years). Defaults to 20,30,45.
    horizons: [max_horizons]u16 = .{ 20, 30, 45 } ++ .{0} ** (max_horizons - 3),
    horizon_count: u8 = 3,
    /// Age-based horizon targets. Resolved at context-load time to
    /// `target_age − max(currentAges())` years — i.e. how long until the
    /// oldest configured person hits `target_age`. Rationale: the first
    /// person to hit the target age sets the meaningful planning horizon,
    /// because spending typically drops substantially after the first death.
    horizon_ages: [max_horizons]u16 = .{0} ** max_horizons,
    horizon_age_count: u8 = 0,
    /// Confidence levels for safe withdrawal. Always 90/95/99.
    confidence_levels: [3]f64 = .{ 0.90, 0.95, 0.99 },
    /// Birthdates for age-based event timing.
    birthdates: [max_persons]Date = .{Date.fromYmd(1970, 1, 1)} ** max_persons,
    birthdate_count: u8 = 0,
    /// Life events (income/expenses) that modify annual cash flow.
    events: [max_events]LifeEvent = undefined,
    event_count: u8 = 0,

    const max_horizons: usize = 8;
    const max_persons: usize = 4;
    pub const max_events: usize = 16;

    /// Errors that can arise when resolving age-based horizons.
    pub const ResolveError = error{
        /// `type::config,horizon_age:num:N` was specified in projections.srf
        /// but no `type::birthdate` record exists to anchor the calculation.
        HorizonAgeWithoutBirthdate,
    };

    pub fn getHorizons(self: *const UserConfig) []const u16 {
        return self.horizons[0..self.horizon_count];
    }

    pub fn getConfidenceLevels(self: *const UserConfig) []const f64 {
        return &self.confidence_levels;
    }

    pub fn getEvents(self: *const UserConfig) []const LifeEvent {
        return self.events[0..self.event_count];
    }

    /// Compute current ages (in whole years) from birthdates.
    pub fn currentAges(self: *const UserConfig) [max_persons]u16 {
        return currentAgesAsOf(self, Date.fromEpoch(std.time.timestamp()));
    }

    /// Compute ages as of a specific date (for testing).
    pub fn currentAgesAsOf(self: *const UserConfig, as_of: Date) [max_persons]u16 {
        var ages: [max_persons]u16 = .{0} ** max_persons;
        for (0..self.birthdate_count) |i| {
            const years = Date.yearsBetween(self.birthdates[i], as_of);
            ages[i] = if (years > 0) @intFromFloat(years) else 0;
        }
        return ages;
    }

    /// Resolve age-based horizons (`horizon_ages`) into year counts and
    /// append them to `horizons`. For each target age, computes
    /// `target_age − max(currentAgesAsOf(as_of))` — the number of years
    /// until the oldest configured person hits that age. Targets that are
    /// already in the past (oldest age ≥ target) are silently skipped.
    ///
    /// Errors if `horizon_ages` is non-empty but no birthdate is configured.
    /// Safe to call multiple times; subsequent calls are no-ops because
    /// `horizon_age_count` is cleared after resolution.
    pub fn resolveHorizonAges(self: *UserConfig, as_of: Date) ResolveError!void {
        if (self.horizon_age_count == 0) return;
        if (self.birthdate_count == 0) return error.HorizonAgeWithoutBirthdate;

        const ages = self.currentAgesAsOf(as_of);
        var oldest: u16 = 0;
        for (0..self.birthdate_count) |i| {
            if (ages[i] > oldest) oldest = ages[i];
        }

        for (0..self.horizon_age_count) |i| {
            const target = self.horizon_ages[i];
            if (target <= oldest) continue; // already past target, skip silently
            const years: u16 = target - oldest;
            if (self.horizon_count < max_horizons) {
                self.horizons[self.horizon_count] = years;
                self.horizon_count += 1;
            }
        }
        // Clear so a second call is a no-op.
        self.horizon_age_count = 0;
    }

    /// Resolve age-based horizons using today's date. Convenience wrapper
    /// around `resolveHorizonAges`.
    pub fn resolveHorizonAgesNow(self: *UserConfig) ResolveError!void {
        return self.resolveHorizonAges(Date.fromEpoch(std.time.timestamp()));
    }

    /// Sum all event cash flows for simulation year `y`.
    pub fn eventNetCashFlow(self: *const UserConfig, y: u16, cumulative_inflation: f64, current_ages: []const u16) f64 {
        var total: f64 = 0;
        for (self.events[0..self.event_count]) |*ev| {
            total += ev.cashFlow(y, cumulative_inflation, current_ages);
        }
        return total;
    }

    /// Resolve all events into ResolvedEvents for the simulation.
    /// Skips events with invalid person indices.
    pub fn resolveEvents(self: *const UserConfig) [max_events]ResolvedEvent {
        const ages = self.currentAges();
        return resolveEventsWithAges(self, &ages);
    }

    /// Resolve all events using pre-computed ages (for testing).
    pub fn resolveEventsWithAges(self: *const UserConfig, ages: []const u16) [max_events]ResolvedEvent {
        var resolved: [max_events]ResolvedEvent = undefined;
        for (self.events[0..self.event_count], 0..) |*ev, i| {
            resolved[i] = ev.resolve(ages) orelse .{
                .start_year = std.math.maxInt(u16), // effectively never active
                .duration = 0,
                .annual_amount = 0,
                .inflation_adjusted = true,
            };
        }
        return resolved;
    }
};

// ── SRF parse types (private) ───────────────────────────────────

const SrfConfig = struct {
    type: []const u8 = "",
    target_stock_pct: ?f64 = null,
    horizon: ?u16 = null,
    horizon_age: ?u16 = null,
};

const SrfBirthdate = struct {
    type: []const u8 = "",
    date: Date,
    person: ?u8 = null, // 1-indexed in SRF; null = sequential
};

const SrfEvent = struct {
    type: []const u8 = "",
    name: []const u8 = "",
    start_age: u16 = 0,
    person: u8 = 1, // 1-indexed in SRF
    duration: u16 = 0,
    amount: f64 = 0,
    inflation_adjusted: bool = true,
};

const SrfProjection = union(enum) {
    pub const srf_tag_field = "type";
    config: SrfConfig,
    birthdate: SrfBirthdate,
    event: SrfEvent,
};

/// Parse a projections.srf file into a UserConfig.
/// Returns default config if data is null or unparseable.
///
/// Format (union-tagged SRF records):
///   type::config,target_stock_pct:num:80
///   type::config,horizon:num:30
///   type::birthdate,date::1975-03-15
///   type::event,name::Social Security,start_age:num:67,amount:num:38400
pub fn parseProjectionsConfig(data: ?[]const u8) UserConfig {
    var config = UserConfig{};
    const raw = data orelse return config;
    if (raw.len == 0) return config;

    var reader = std.Io.Reader.fixed(raw);
    var it = srf.iterator(&reader, std.heap.smp_allocator, .{ .alloc_strings = false }) catch return config;
    defer it.deinit();

    var saw_horizon = false;
    var birthdate_seq: u8 = 0;

    while (it.next() catch null) |field_it| {
        const rec = field_it.to(SrfProjection) catch continue;
        switch (rec) {
            .config => |c| {
                if (c.target_stock_pct) |val| {
                    config.target_stock_pct = val;
                }
                if (c.horizon) |h| {
                    if (!saw_horizon) {
                        config.horizon_count = 0;
                        saw_horizon = true;
                    }
                    if (config.horizon_count < UserConfig.max_horizons and h > 0) {
                        config.horizons[config.horizon_count] = h;
                        config.horizon_count += 1;
                    }
                }
                if (c.horizon_age) |age| {
                    // Age-based horizons are stored raw and resolved later
                    // via `UserConfig.resolveHorizonAges(as_of)` once the
                    // view layer knows the projection date. They also count
                    // as "saw_horizon" so a file containing only
                    // `horizon_age` records replaces the default {20,30,45}
                    // once resolved.
                    if (!saw_horizon) {
                        config.horizon_count = 0;
                        saw_horizon = true;
                    }
                    if (config.horizon_age_count < UserConfig.max_horizons and age > 0) {
                        config.horizon_ages[config.horizon_age_count] = age;
                        config.horizon_age_count += 1;
                    }
                }
            },
            .birthdate => |b| {
                // person is 1-indexed in SRF; convert to 0-indexed.
                // If not specified, assign sequentially.
                const idx: u8 = if (b.person) |p| p -| 1 else birthdate_seq;
                if (idx < UserConfig.max_persons) {
                    config.birthdates[idx] = b.date;
                    if (idx >= config.birthdate_count) config.birthdate_count = idx + 1;
                }
                birthdate_seq += 1;
            },
            .event => |e| {
                if (config.event_count < UserConfig.max_events and e.start_age > 0) {
                    var ev = LifeEvent{
                        .start_age = e.start_age,
                        .person = e.person -| 1, // 1-indexed → 0-indexed
                        .duration = e.duration,
                        .annual_amount = e.amount,
                        .inflation_adjusted = e.inflation_adjusted,
                    };
                    const len = @min(e.name.len, LifeEvent.max_name_len);
                    @memcpy(ev.name[0..len], e.name[0..len]);
                    ev.name_len = @intCast(len);
                    config.events[config.event_count] = ev;
                    config.event_count += 1;
                }
            },
        }
    }

    return config;
}

// ── Configuration ──────────────────────────────────────────────

/// Conservative return estimation defaults.
/// These are module-level constants that will eventually move to projections.srf.
pub const default_return_cap: ?f64 = null; // no cap currently
pub const default_exclude_1y_from_min: bool = true; // use MIN(3Y, 5Y, 10Y), skip 1Y

pub const ProjectionConfig = struct {
    /// Current total portfolio value in dollars.
    portfolio_value: f64,
    /// Stock allocation as a fraction (0.0–1.0). Remainder goes to bonds.
    stock_pct: f64,
    /// Retirement time horizons to simulate (in years).
    horizons: []const u16,
    /// Confidence levels for safe withdrawal (e.g. 0.90, 0.95, 0.99).
    confidence_levels: []const f64,
    /// Pre-resolved life events for the simulation.
    events: []const ResolvedEvent = &.{},
};

// ── Results ────────────────────────────────────────────────────

pub const WithdrawalResult = struct {
    /// Confidence level (e.g. 0.99 = 99%).
    confidence: f64,
    /// Maximum annual withdrawal that achieves this confidence.
    annual_amount: f64,
    /// As a fraction of starting portfolio value.
    withdrawal_rate: f64,
};

pub const YearPercentiles = struct {
    /// Year offset from retirement start (0 = start, 1 = after year 1, etc.)
    year: u16,
    p10: f64,
    p25: f64,
    p50: f64,
    p75: f64,
    p90: f64,
};

pub const SimulationResult = struct {
    horizon: u16,
    /// Number of historical cycles simulated.
    num_cycles: usize,
    /// Safe withdrawal amounts at each requested confidence level.
    withdrawals: []WithdrawalResult,
    /// Portfolio value percentiles at each year (for charting).
    /// Length = horizon + 1 (includes year 0 = starting value).
    percentile_bands: []YearPercentiles,

    pub fn deinit(self: *SimulationResult, allocator: std.mem.Allocator) void {
        allocator.free(self.withdrawals);
        allocator.free(self.percentile_bands);
    }
};

// ── Core simulation ────────────────────────────────────────────

/// Simulate a single retirement cycle starting at `start_index` in the
/// Shiller dataset, lasting `horizon` years, with the given annual spending
/// (inflation-adjusted) and stock/bond allocation.
///
/// Follows the FIRECalc methodology:
///   1. Withdraw spending (year 1 = base amount, subsequent years CPI-adjusted)
///   2. Apply market returns to the remainder
///   3. Repeat
///
/// Returns the portfolio values at each year (length = horizon + 1).
fn simulateCycle(
    buf: []f64,
    start_index: usize,
    horizon: u16,
    initial_value: f64,
    annual_spending: f64,
    stock_pct: f64,
    events: []const ResolvedEvent,
) void {
    const data = shiller.annual_returns;
    var portfolio = initial_value;
    buf[0] = portfolio;

    var cumulative_inflation: f64 = 1.0;

    for (0..horizon) |y| {
        const di = start_index + y;
        if (di >= data.len) {
            for (y + 1..@as(usize, horizon) + 1) |remaining| {
                buf[remaining] = portfolio;
            }
            return;
        }

        const yr = data[di];

        // Step 1: Withdraw spending, offset by life event cash flows
        var event_net: f64 = 0;
        for (events) |*ev| {
            event_net += ev.cashFlow(@intCast(y), cumulative_inflation);
        }
        portfolio -= annual_spending * cumulative_inflation - event_net;

        // Step 2: Apply market returns (nominal stock + nominal bond via GS10 yield)
        const blended_return = stock_pct * yr.sp500_total_return + (1.0 - stock_pct) * yr.bond_total_return;
        portfolio *= (1.0 + blended_return);

        // Step 3: Update cumulative inflation for next year's spending
        cumulative_inflation *= (1.0 + yr.cpi_inflation);

        buf[y + 1] = portfolio;
    }
}

/// Run the full historical simulation: for each possible starting year,
/// simulate a retirement of `horizon` years. Returns the number of cycles
/// where the portfolio survived (never went to zero).
///
/// `all_paths` is a 2D buffer: [cycle_index][year] = portfolio value.
/// Must be pre-allocated with dimensions [num_cycles][horizon + 1].
fn runAllCycles(
    all_paths: [][]f64,
    horizon: u16,
    initial_value: f64,
    annual_spending: f64,
    stock_pct: f64,
    events: []const ResolvedEvent,
) usize {
    const num_cycles = shiller.maxCycles(horizon);
    var survived: usize = 0;

    for (0..num_cycles) |cycle| {
        simulateCycle(
            all_paths[cycle],
            cycle,
            horizon,
            initial_value,
            annual_spending,
            stock_pct,
            events,
        );

        // Check if portfolio survived the full horizon
        var failed = false;
        for (1..@as(usize, horizon) + 1) |y| {
            if (all_paths[cycle][y] <= 0) {
                // Zero out remaining years once the portfolio is depleted
                for (y..@as(usize, horizon) + 1) |z| {
                    all_paths[cycle][z] = 0;
                }
                failed = true;
                break;
            }
        }
        if (!failed) survived += 1;
    }

    return survived;
}

/// Compute the success rate (fraction of cycles that survived) for a
/// given spending level. Lightweight version that doesn't store full paths.
fn successRate(
    horizon: u16,
    initial_value: f64,
    annual_spending: f64,
    stock_pct: f64,
    events: []const ResolvedEvent,
) f64 {
    const data = shiller.annual_returns;
    const num_cycles = shiller.maxCycles(horizon);
    if (num_cycles == 0) return 0.0;

    var survived: usize = 0;

    for (0..num_cycles) |cycle| {
        var portfolio = initial_value;
        var cumulative_inflation: f64 = 1.0;
        var failed = false;

        for (0..horizon) |y| {
            const di = cycle + y;
            if (di >= data.len) break;

            const yr = data[di];

            // Withdraw spending, offset by life event cash flows
            var event_net: f64 = 0;
            for (events) |*ev| {
                event_net += ev.cashFlow(@intCast(y), cumulative_inflation);
            }
            portfolio -= annual_spending * cumulative_inflation - event_net;
            if (portfolio <= 0) {
                failed = true;
                break;
            }

            // Then grow
            const blended_return = stock_pct * yr.sp500_total_return + (1.0 - stock_pct) * yr.bond_total_return;
            portfolio *= (1.0 + blended_return);

            // Update inflation for next year
            cumulative_inflation *= (1.0 + yr.cpi_inflation);
        }

        if (!failed) survived += 1;
    }

    return @as(f64, @floatFromInt(survived)) / @as(f64, @floatFromInt(num_cycles));
}

// ── Safe withdrawal search ─────────────────────────────────────

/// Find the maximum annual withdrawal amount (in today's dollars) such that
/// the portfolio survives `horizon` years in at least `confidence` fraction
/// of all historical cycles.
///
/// Uses binary search with $1 precision, seeded with a 4%-rule estimate
/// to narrow the search band (~10 iterations instead of ~23).
pub fn findSafeWithdrawal(
    horizon: u16,
    initial_value: f64,
    stock_pct: f64,
    confidence: f64,
    events: []const ResolvedEvent,
) WithdrawalResult {
    // Seed from the 4% rule, adjusted for horizon and confidence.
    // Base ~4% for 30yr/95%. Shorter horizons allow more; longer less.
    // Higher confidence requires less.
    const base_rate = 0.04;
    const horizon_adj = 30.0 / @as(f64, @floatFromInt(horizon)); // >1 for short, <1 for long
    const conf_adj = (1.0 - confidence) / 0.05; // 1.0 at 95%, 0.2 at 99%, 2.0 at 90%
    const estimate = initial_value * base_rate * @sqrt(horizon_adj) * @sqrt(conf_adj);

    // Search band: ±50% of estimate, clamped to [0, initial_value]
    var lo: f64 = @max(estimate * 0.5, 0);
    var hi: f64 = @min(estimate * 1.5, initial_value);

    // Verify bounds bracket the answer; widen if not
    if (successRate(horizon, initial_value, lo, stock_pct, events) < confidence) {
        log.debug("findSafeWithdrawal: estimate too high, widening lo to 0 (horizon={d}, conf={d:.2})", .{ horizon, confidence });
        lo = 0;
    }
    if (successRate(horizon, initial_value, hi, stock_pct, events) >= confidence) {
        log.debug("findSafeWithdrawal: estimate too low, widening hi to portfolio value (horizon={d}, conf={d:.2})", .{ horizon, confidence });
        hi = initial_value;
    }

    // Binary search to $1 precision
    while (hi - lo > 1.0) {
        const mid = @floor((lo + hi) / 2.0);
        const rate = successRate(horizon, initial_value, mid, stock_pct, events);

        if (rate >= confidence) {
            lo = mid;
        } else {
            hi = mid;
        }
    }

    return .{
        .confidence = confidence,
        .annual_amount = lo,
        .withdrawal_rate = lo / initial_value,
    };
}

// ── Percentile bands ───────────────────────────────────────────

/// Compute percentile bands from all simulated paths for a given horizon
/// and spending level. Allocates the result.
pub fn computePercentileBands(
    allocator: std.mem.Allocator,
    horizon: u16,
    initial_value: f64,
    annual_spending: f64,
    stock_pct: f64,
    events: []const ResolvedEvent,
) ![]YearPercentiles {
    const num_cycles = shiller.maxCycles(horizon);
    if (num_cycles == 0) return &.{};

    const years: usize = @as(usize, horizon) + 1;

    // Allocate path storage: num_cycles rows of (horizon+1) f64s
    const path_data = try allocator.alloc(f64, num_cycles * years);
    defer allocator.free(path_data);

    const paths = try allocator.alloc([]f64, num_cycles);
    defer allocator.free(paths);

    for (0..num_cycles) |i| {
        paths[i] = path_data[i * years .. (i + 1) * years];
    }

    _ = runAllCycles(paths, horizon, initial_value, annual_spending, stock_pct, events);

    // For each year, sort the values across all cycles and extract percentiles
    const bands = try allocator.alloc(YearPercentiles, years);

    // Temporary buffer for sorting one year's values
    const sort_buf = try allocator.alloc(f64, num_cycles);
    defer allocator.free(sort_buf);

    for (0..years) |y| {
        // Collect values for year y across all cycles
        for (0..num_cycles) |c| {
            sort_buf[c] = paths[c][y];
        }

        // Sort
        std.mem.sort(f64, sort_buf, {}, std.sort.asc(f64));

        bands[y] = .{
            .year = @intCast(y),
            .p10 = percentile(sort_buf, 0.10),
            .p25 = percentile(sort_buf, 0.25),
            .p50 = percentile(sort_buf, 0.50),
            .p75 = percentile(sort_buf, 0.75),
            .p90 = percentile(sort_buf, 0.90),
        };
    }

    return bands;
}

/// Linear interpolation percentile on a sorted slice.
fn percentile(sorted: []const f64, p: f64) f64 {
    if (sorted.len == 0) return 0;
    if (sorted.len == 1) return sorted[0];

    const n = @as(f64, @floatFromInt(sorted.len - 1));
    const idx = p * n;
    const lo_idx: usize = @intFromFloat(@floor(idx));
    const hi_idx: usize = @min(lo_idx + 1, sorted.len - 1);
    const frac = idx - @floor(idx);

    return sorted[lo_idx] * (1.0 - frac) + sorted[hi_idx] * frac;
}

// ── High-level API ─────────────────────────────────────────────

/// Run the full projection analysis for one horizon: compute safe withdrawal
/// at each confidence level and percentile bands using the median withdrawal.
pub fn runProjection(
    allocator: std.mem.Allocator,
    config: ProjectionConfig,
    horizon: u16,
) !SimulationResult {
    const num_cycles = shiller.maxCycles(horizon);

    // Compute safe withdrawal at each confidence level
    const withdrawals = try allocator.alloc(WithdrawalResult, config.confidence_levels.len);
    for (config.confidence_levels, 0..) |conf, i| {
        withdrawals[i] = findSafeWithdrawal(
            horizon,
            config.portfolio_value,
            config.stock_pct,
            conf,
            config.events,
        );
    }

    // Use the median confidence level's withdrawal for the percentile chart
    const median_idx = config.confidence_levels.len / 2;
    const chart_spending = if (withdrawals.len > 0) withdrawals[median_idx].annual_amount else 0;

    const bands = try computePercentileBands(
        allocator,
        horizon,
        config.portfolio_value,
        chart_spending,
        config.stock_pct,
        config.events,
    );

    return .{
        .horizon = horizon,
        .num_cycles = num_cycles,
        .withdrawals = withdrawals,
        .percentile_bands = bands,
    };
}

/// Run projections for all configured horizons.
pub fn runAllProjections(
    allocator: std.mem.Allocator,
    config: ProjectionConfig,
) ![]SimulationResult {
    const results = try allocator.alloc(SimulationResult, config.horizons.len);
    errdefer {
        for (results) |*r| r.deinit(allocator);
        allocator.free(results);
    }

    for (config.horizons, 0..) |h, i| {
        results[i] = try runProjection(allocator, config, h);
    }

    return results;
}

// ── Tests ──────────────────────────────────────────────────────

test "successRate with zero spending is 100%" {
    const rate = successRate(30, 1_000_000, 0, 0.75, &.{});
    try std.testing.expectApproxEqAbs(@as(f64, 1.0), rate, 0.001);
}

test "successRate with excessive spending is 0%" {
    // Spending the entire portfolio in year 1 should fail every cycle
    const rate = successRate(30, 1_000_000, 1_000_000, 0.75, &.{});
    try std.testing.expectApproxEqAbs(@as(f64, 0.0), rate, 0.001);
}

test "successRate decreases with higher spending" {
    const rate_low = successRate(30, 1_000_000, 20_000, 0.75, &.{});
    const rate_mid = successRate(30, 1_000_000, 40_000, 0.75, &.{});
    const rate_high = successRate(30, 1_000_000, 60_000, 0.75, &.{});

    try std.testing.expect(rate_low >= rate_mid);
    try std.testing.expect(rate_mid >= rate_high);
}

test "findSafeWithdrawal produces reasonable results" {
    const result = findSafeWithdrawal(30, 1_000_000, 0.75, 0.95, &.{});

    try std.testing.expect(result.annual_amount >= 10_000);
    try std.testing.expect(result.annual_amount <= 60_000);
    try std.testing.expect(result.withdrawal_rate >= 0.01);
    try std.testing.expect(result.withdrawal_rate <= 0.06);
}

test "higher confidence means lower withdrawal" {
    const r90 = findSafeWithdrawal(30, 1_000_000, 0.75, 0.90, &.{});
    const r95 = findSafeWithdrawal(30, 1_000_000, 0.75, 0.95, &.{});
    const r99 = findSafeWithdrawal(30, 1_000_000, 0.75, 0.99, &.{});

    try std.testing.expect(r90.annual_amount >= r95.annual_amount);
    try std.testing.expect(r95.annual_amount >= r99.annual_amount);
}

test "longer horizon means lower withdrawal" {
    const r20 = findSafeWithdrawal(20, 1_000_000, 0.75, 0.95, &.{});
    const r30 = findSafeWithdrawal(30, 1_000_000, 0.75, 0.95, &.{});
    const r45 = findSafeWithdrawal(45, 1_000_000, 0.75, 0.95, &.{});

    try std.testing.expect(r20.annual_amount >= r30.annual_amount);
    try std.testing.expect(r30.annual_amount >= r45.annual_amount);
}

test "computePercentileBands basic properties" {
    const allocator = std.testing.allocator;
    const bands = try computePercentileBands(allocator, 30, 1_000_000, 30_000, 0.75, &.{});
    defer allocator.free(bands);

    // Should have horizon + 1 entries
    try std.testing.expectEqual(@as(usize, 31), bands.len);

    // Year 0 should be the starting value for all percentiles
    try std.testing.expectApproxEqAbs(@as(f64, 1_000_000), bands[0].p50, 1.0);

    // Percentiles should be ordered at each year
    for (bands) |b| {
        try std.testing.expect(b.p10 <= b.p25);
        try std.testing.expect(b.p25 <= b.p50);
        try std.testing.expect(b.p50 <= b.p75);
        try std.testing.expect(b.p75 <= b.p90);
    }
}

test "runProjection produces valid results" {
    const allocator = std.testing.allocator;
    const config = ProjectionConfig{
        .portfolio_value = 1_000_000,
        .stock_pct = 0.75,
        .horizons = &.{ 20, 30 },
        .confidence_levels = &.{ 0.90, 0.95, 0.99 },
    };

    var result = try runProjection(allocator, config, 30);
    defer result.deinit(allocator);

    try std.testing.expectEqual(@as(u16, 30), result.horizon);
    try std.testing.expectEqual(@as(usize, 3), result.withdrawals.len);
    try std.testing.expectEqual(@as(usize, 31), result.percentile_bands.len);

    // Withdrawals should be ordered: 90% > 95% > 99%
    try std.testing.expect(result.withdrawals[0].annual_amount >= result.withdrawals[1].annual_amount);
    try std.testing.expect(result.withdrawals[1].annual_amount >= result.withdrawals[2].annual_amount);
}

test "percentile interpolation" {
    const data = [_]f64{ 10, 20, 30, 40, 50 };
    try std.testing.expectApproxEqAbs(@as(f64, 10.0), percentile(&data, 0.0), 0.01);
    try std.testing.expectApproxEqAbs(@as(f64, 30.0), percentile(&data, 0.5), 0.01);
    try std.testing.expectApproxEqAbs(@as(f64, 50.0), percentile(&data, 1.0), 0.01);
    try std.testing.expectApproxEqAbs(@as(f64, 20.0), percentile(&data, 0.25), 0.01);
}

test "realistic portfolio safe withdrawal" {
    // Approximate real portfolio: ~$8.34M, ~82.5% stocks
    const portfolio = 8_340_000;
    const stock_pct = 0.825;

    const r99_45 = findSafeWithdrawal(45, portfolio, stock_pct, 0.99, &.{});
    const r95_45 = findSafeWithdrawal(45, portfolio, stock_pct, 0.95, &.{});
    const r99_30 = findSafeWithdrawal(30, portfolio, stock_pct, 0.99, &.{});

    // 95% should be higher than 99%
    try std.testing.expect(r95_45.annual_amount > r99_45.annual_amount);
    // 30yr should be higher than 45yr at same confidence
    try std.testing.expect(r99_30.annual_amount > r99_45.annual_amount);
    // Should produce $290K+ at 99%/45yr based on FIRECalc reference (~$305K on $7.7M)
    try std.testing.expect(r99_45.annual_amount >= 290_000);
    try std.testing.expect(r99_45.annual_amount <= 350_000);
    try std.testing.expect(r99_45.withdrawal_rate >= 0.03);
    try std.testing.expect(r99_45.withdrawal_rate <= 0.05);
}

test "simulateCycle produces correct year-0 value" {
    var buf: [31]f64 = undefined;
    simulateCycle(&buf, 0, 30, 1_000_000, 0, 0.75, &.{});
    try std.testing.expectApproxEqAbs(@as(f64, 1_000_000), buf[0], 0.01);
}

test "simulateCycle with zero spending grows portfolio" {
    var buf: [31]f64 = undefined;
    simulateCycle(&buf, 0, 30, 1_000_000, 0, 0.75, &.{});
    // Over any 30-year period in history, zero spending should grow the portfolio
    try std.testing.expect(buf[30] > 1_000_000);
}

test "parseProjectionsConfig defaults" {
    const config = parseProjectionsConfig(null);
    try std.testing.expect(config.target_stock_pct == null);
    try std.testing.expectEqual(@as(u8, 3), config.horizon_count);
    try std.testing.expectEqual(@as(u16, 20), config.getHorizons()[0]);
    try std.testing.expectEqual(@as(u16, 30), config.getHorizons()[1]);
    try std.testing.expectEqual(@as(u16, 45), config.getHorizons()[2]);
}

test "parseProjectionsConfig from SRF" {
    const data =
        \\#!srfv1
        \\type::config,target_stock_pct:num:77
        \\type::config,horizon:num:25
        \\type::config,horizon:num:35
        \\type::config,horizon:num:50
    ;
    const config = parseProjectionsConfig(data);
    try std.testing.expectApproxEqAbs(@as(f64, 77.0), config.target_stock_pct.?, 0.01);
    try std.testing.expectEqual(@as(u8, 3), config.horizon_count);
    try std.testing.expectEqual(@as(u16, 25), config.getHorizons()[0]);
    try std.testing.expectEqual(@as(u16, 35), config.getHorizons()[1]);
    try std.testing.expectEqual(@as(u16, 50), config.getHorizons()[2]);
}

test "parseProjectionsConfig partial" {
    const data = "#!srfv1\ntype::config,target_stock_pct:num:82.5\n";
    const config = parseProjectionsConfig(data);
    try std.testing.expectApproxEqAbs(@as(f64, 82.5), config.target_stock_pct.?, 0.01);
    // Horizons should remain default
    try std.testing.expectEqual(@as(u8, 3), config.horizon_count);
    try std.testing.expectEqual(@as(u16, 20), config.getHorizons()[0]);
}

test "parseProjectionsConfig empty string" {
    const config = parseProjectionsConfig("");
    try std.testing.expect(config.target_stock_pct == null);
    try std.testing.expectEqual(@as(u8, 3), config.horizon_count);
}

test "parseProjectionsConfig invalid data" {
    const config = parseProjectionsConfig("not valid srf");
    try std.testing.expect(config.target_stock_pct == null);
}

test "parseProjectionsConfig horizon_age parsed raw" {
    const data =
        \\#!srfv1
        \\type::config,horizon_age:num:90
        \\type::config,horizon_age:num:95
        \\type::birthdate,date::1975-03-15
    ;
    const config = parseProjectionsConfig(data);
    // horizon_ages are stored raw; not yet resolved into horizons.
    try std.testing.expectEqual(@as(u8, 2), config.horizon_age_count);
    try std.testing.expectEqual(@as(u16, 90), config.horizon_ages[0]);
    try std.testing.expectEqual(@as(u16, 95), config.horizon_ages[1]);
    // A horizon_age record counts as "saw_horizon", so the default
    // {20,30,45} is cleared. horizon_count is 0 until resolution.
    try std.testing.expectEqual(@as(u8, 0), config.horizon_count);
}

test "resolveHorizonAges uses oldest birthdate (first-to-hit semantics)" {
    // Person 1: born 1975, age 50 as of 2025. Person 2: born 1980, age 45.
    // Target age 90 → 90 − 50 = 40 years (first to hit 90 is the older).
    var config = parseProjectionsConfig(
        \\#!srfv1
        \\type::config,horizon_age:num:90
        \\type::birthdate,date::1975-06-15
        \\type::birthdate,date::1980-06-15,person:num:2
    );
    const as_of = Date.fromYmd(2025, 6, 15);
    try config.resolveHorizonAges(as_of);
    try std.testing.expectEqual(@as(u8, 1), config.horizon_count);
    try std.testing.expectEqual(@as(u16, 40), config.horizons[0]);
    // Resolved; horizon_age_count cleared to make resolve idempotent.
    try std.testing.expectEqual(@as(u8, 0), config.horizon_age_count);
}

test "resolveHorizonAges errors without a birthdate" {
    var config = parseProjectionsConfig(
        \\#!srfv1
        \\type::config,horizon_age:num:90
    );
    const as_of = Date.fromYmd(2025, 1, 1);
    try std.testing.expectError(error.HorizonAgeWithoutBirthdate, config.resolveHorizonAges(as_of));
}

test "resolveHorizonAges skips targets already in the past" {
    // Oldest age is 60 as of 2025; target 40 is already past — skipped.
    var config = parseProjectionsConfig(
        \\#!srfv1
        \\type::config,horizon_age:num:40
        \\type::config,horizon_age:num:90
        \\type::birthdate,date::1965-01-01
    );
    const as_of = Date.fromYmd(2025, 6, 15);
    try config.resolveHorizonAges(as_of);
    // Only age 90 resolves (90 − 60 = 30).
    try std.testing.expectEqual(@as(u8, 1), config.horizon_count);
    try std.testing.expectEqual(@as(u16, 30), config.horizons[0]);
}

test "resolveHorizonAges mixes with explicit horizon records" {
    var config = parseProjectionsConfig(
        \\#!srfv1
        \\type::config,horizon:num:30
        \\type::config,horizon_age:num:95
        \\type::birthdate,date::1975-06-15
    );
    const as_of = Date.fromYmd(2025, 6, 15);
    try config.resolveHorizonAges(as_of);
    // Explicit 30 from `horizon`, then appended 95 − 50 = 45 from `horizon_age`.
    try std.testing.expectEqual(@as(u8, 2), config.horizon_count);
    try std.testing.expectEqual(@as(u16, 30), config.horizons[0]);
    try std.testing.expectEqual(@as(u16, 45), config.horizons[1]);
}

test "resolveHorizonAges is a no-op when nothing to resolve" {
    var config = parseProjectionsConfig(
        \\#!srfv1
        \\type::config,horizon:num:30
    );
    // No birthdate, no horizon_age → should succeed, not error.
    const as_of = Date.fromYmd(2025, 1, 1);
    try config.resolveHorizonAges(as_of);
    try std.testing.expectEqual(@as(u8, 1), config.horizon_count);
    try std.testing.expectEqual(@as(u16, 30), config.horizons[0]);
}

test "UserConfig getHorizons default" {
    const config = UserConfig{};
    const horizons = config.getHorizons();
    try std.testing.expectEqual(@as(usize, 3), horizons.len);
    try std.testing.expectEqual(@as(u16, 20), horizons[0]);
    try std.testing.expectEqual(@as(u16, 30), horizons[1]);
    try std.testing.expectEqual(@as(u16, 45), horizons[2]);
}

test "UserConfig getConfidenceLevels" {
    const config = UserConfig{};
    const levels = config.getConfidenceLevels();
    try std.testing.expectEqual(@as(usize, 3), levels.len);
    try std.testing.expectApproxEqAbs(@as(f64, 0.90), levels[0], 0.001);
    try std.testing.expectApproxEqAbs(@as(f64, 0.95), levels[1], 0.001);
    try std.testing.expectApproxEqAbs(@as(f64, 0.99), levels[2], 0.001);
}

test "runAllProjections produces results for each horizon" {
    const allocator = std.testing.allocator;
    const config = ProjectionConfig{
        .portfolio_value = 1_000_000,
        .stock_pct = 0.75,
        .horizons = &.{ 20, 30 },
        .confidence_levels = &.{ 0.95, 0.99 },
    };
    const results = try runAllProjections(allocator, config);
    defer {
        for (results) |*r| r.deinit(allocator);
        allocator.free(results);
    }

    try std.testing.expectEqual(@as(usize, 2), results.len);
    try std.testing.expectEqual(@as(u16, 20), results[0].horizon);
    try std.testing.expectEqual(@as(u16, 30), results[1].horizon);
    try std.testing.expectEqual(@as(usize, 2), results[0].withdrawals.len);
    try std.testing.expectEqual(@as(usize, 2), results[1].withdrawals.len);
}

test "LifeEvent.startYear basic" {
    const ev = LifeEvent{ .start_age = 67, .person = 0, .annual_amount = 38400 };
    const ages = [_]u16{50};
    try std.testing.expectEqual(@as(?u16, 17), ev.startYear(&ages));
}

test "LifeEvent.startYear already active" {
    const ev = LifeEvent{ .start_age = 40, .person = 0, .annual_amount = 38400 };
    const ages = [_]u16{50};
    try std.testing.expectEqual(@as(?u16, 0), ev.startYear(&ages));
}

test "LifeEvent.startYear person out of range" {
    const ev = LifeEvent{ .start_age = 67, .person = 5, .annual_amount = 38400 };
    const ages = [_]u16{50};
    try std.testing.expectEqual(@as(?u16, null), ev.startYear(&ages));
}

test "LifeEvent.isActive permanent" {
    const ev = LifeEvent{ .start_age = 60, .person = 0, .duration = 0, .annual_amount = 38400 };
    const ages = [_]u16{50};
    try std.testing.expect(!ev.isActive(9, &ages)); // before start (year 10)
    try std.testing.expect(ev.isActive(10, &ages)); // start year
    try std.testing.expect(ev.isActive(30, &ages)); // well after
}

test "LifeEvent.isActive with duration" {
    const ev = LifeEvent{ .start_age = 53, .person = 0, .duration = 4, .annual_amount = -60000 };
    const ages = [_]u16{50};
    try std.testing.expect(!ev.isActive(2, &ages)); // before start
    try std.testing.expect(ev.isActive(3, &ages)); // year 3 (age 53)
    try std.testing.expect(ev.isActive(6, &ages)); // year 6 (age 56, last year)
    try std.testing.expect(!ev.isActive(7, &ages)); // year 7 (age 57, past duration)
}

test "LifeEvent.cashFlow inflation adjusted" {
    const ev = LifeEvent{ .start_age = 50, .person = 0, .annual_amount = 10000, .inflation_adjusted = true };
    const ages = [_]u16{50};
    try std.testing.expectApproxEqAbs(@as(f64, 12000), ev.cashFlow(0, 1.2, &ages), 0.01);
}

test "LifeEvent.cashFlow nominal" {
    const ev = LifeEvent{ .start_age = 50, .person = 0, .annual_amount = 10000, .inflation_adjusted = false };
    const ages = [_]u16{50};
    try std.testing.expectApproxEqAbs(@as(f64, 10000), ev.cashFlow(0, 1.2, &ages), 0.01);
}

test "LifeEvent.cashFlow inactive returns zero" {
    const ev = LifeEvent{ .start_age = 67, .person = 0, .annual_amount = 38400 };
    const ages = [_]u16{50};
    try std.testing.expectApproxEqAbs(@as(f64, 0), ev.cashFlow(5, 1.0, &ages), 0.01);
}

test "parseProjectionsConfig birthdates and events" {
    const data =
        \\#!srfv1
        \\type::config,target_stock_pct:num:80
        \\type::config,horizon:num:30
        \\type::birthdate,date::1975-03-15
        \\type::birthdate,date::1978-06-22,person:num:2
        \\type::event,name::Social Security,start_age:num:67,person:num:1,amount:num:38400
        \\type::event,name::College,start_age:num:53,duration:num:4,amount:num:-60000,inflation_adjusted:bool:false
    ;
    const config = parseProjectionsConfig(data);
    try std.testing.expectApproxEqAbs(@as(f64, 80.0), config.target_stock_pct.?, 0.01);
    try std.testing.expectEqual(@as(u8, 1), config.horizon_count);
    try std.testing.expectEqual(@as(u8, 2), config.birthdate_count);
    try std.testing.expectEqual(@as(i16, 1975), config.birthdates[0].year());
    try std.testing.expectEqual(@as(i16, 1978), config.birthdates[1].year());
    try std.testing.expectEqual(@as(u8, 2), config.event_count);
    // First event: Social Security
    const ev0 = config.events[0];
    try std.testing.expectEqualStrings("Social Security", ev0.getName());
    try std.testing.expectEqual(@as(u16, 67), ev0.start_age);
    try std.testing.expectEqual(@as(u8, 0), ev0.person);
    try std.testing.expectEqual(@as(u16, 0), ev0.duration);
    try std.testing.expectApproxEqAbs(@as(f64, 38400), ev0.annual_amount, 0.01);
    try std.testing.expect(ev0.inflation_adjusted);
    // Second event: College
    const ev1 = config.events[1];
    try std.testing.expectEqualStrings("College", ev1.getName());
    try std.testing.expectEqual(@as(u16, 53), ev1.start_age);
    try std.testing.expectEqual(@as(u16, 4), ev1.duration);
    try std.testing.expectApproxEqAbs(@as(f64, -60000), ev1.annual_amount, 0.01);
    try std.testing.expect(!ev1.inflation_adjusted);
}

test "income event increases safe withdrawal" {
    // With a permanent $20K/yr income event starting immediately,
    // the safe withdrawal should be higher than without.
    const no_events = findSafeWithdrawal(30, 1_000_000, 0.75, 0.95, &.{});
    const income_event = [_]ResolvedEvent{.{
        .start_year = 0,
        .duration = 0,
        .annual_amount = 20_000,
        .inflation_adjusted = true,
    }};
    const with_income = findSafeWithdrawal(30, 1_000_000, 0.75, 0.95, &income_event);
    try std.testing.expect(with_income.annual_amount > no_events.annual_amount);
    // The increase should be roughly $20K (the income offsets withdrawal)
    const diff = with_income.annual_amount - no_events.annual_amount;
    try std.testing.expect(diff >= 15_000 and diff <= 25_000);
}

test "expense event decreases safe withdrawal" {
    const no_events = findSafeWithdrawal(30, 1_000_000, 0.75, 0.95, &.{});
    const expense_event = [_]ResolvedEvent{.{
        .start_year = 0,
        .duration = 5,
        .annual_amount = -20_000,
        .inflation_adjusted = true,
    }};
    const with_expense = findSafeWithdrawal(30, 1_000_000, 0.75, 0.95, &expense_event);
    try std.testing.expect(with_expense.annual_amount < no_events.annual_amount);
}

test "UserConfig.eventNetCashFlow sums active events" {
    var config = UserConfig{};
    config.birthdate_count = 1;
    config.birthdates[0] = Date.fromYmd(1975, 1, 1);
    config.events[0] = .{ .start_age = 50, .person = 0, .annual_amount = 30000 };
    config.events[1] = .{ .start_age = 55, .person = 0, .annual_amount = 10000 };
    config.event_count = 2;
    const ages = [_]u16{50};
    // At year 0: only first event active (age 50)
    try std.testing.expectApproxEqAbs(@as(f64, 30000), config.eventNetCashFlow(0, 1.0, &ages), 0.01);
    // At year 5: both active (ages 55, 55)
    try std.testing.expectApproxEqAbs(@as(f64, 40000), config.eventNetCashFlow(5, 1.0, &ages), 0.01);
}