Monte-Carlo-Simulation (Roadmap Nr. 19, Stufe A)

Button in der Planansicht -> Dialog: Erklaerung, Eingaben, Lauf, Ergebnis + Faecher.
Statt einer festen Rendite/Inflation werden tausende Zufallspfade gerechnet und die
Erfolgs-/Ruinwahrscheinlichkeit des Plans ausgewiesen.

Kern (computePlan-Nahtstelle):
- computePlan(plan, sample?) nimmt optional pro Jahr Inflation und pro Element/Jahr
  eine Rendite. Ohne sample bitgenau wie bisher (durch Golden Tests abgesichert).
- Dazu Inflation auf ein kumulatives Deflator-Array umgestellt (statt (1+i)^t), damit
  sie pro Jahr variieren kann. Deterministisch identisch.
- Reine Funktion ohne Server-Deps -> Simulation laeuft komplett im Browser, null
  Serverlast. ~10'000 Laeufe in ~1 s, Fortschritt alle 500 Laeufe (kein Freeze).

Statistik (montecarlo.ts):
- Fettschwaenzig (standardisierte Student-t, nu=5): Extremcrashs realistisch haeufig,
  eine Normalverteilung wuerde sie stark unterschaetzen.
- Gemeinsamer Marktschock (rho=0.7): riskante Anlagen fallen im Crash zusammen,
  nicht gegeneinander.
- Boeden: 0 % fuer PK/3a, -100 % sonst. Seedbar (reproduzierbar).
- Zwei Renditezahlen: geplante (Zielbalken) vs. historische (Streu-Mittelpunkt) --
  sonst waere P(>= geplantes Endvermoegen) immer ~50 %.

UI: Erklaerung, pro Element/Inflation historischer Oe (Pflicht, kein Default) +
Streuungsstufe (recherchierte sigma-Werte) + Anzahl Laeufe. Ergebnis: Ruin-/
Erfolgswahrscheinlichkeit, Faecher (10/50/90) + deterministische Linie.

7 MC-Tests inkl. deterministischer Aequivalenz, Vol-Drag, Reproduzierbarkeit, Boden,
Ruin (41 -> 48). Keine DB-Aenderung. Spezifikation auf v0.7 (Kap. 4.12, 9.15).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

src/lib/calculations.ts

@@ -210,7 +210,16 @@ function personByRole<T extends { role: PersonRole }>(persons: T[], role: string
   return persons.find((p) => p.role === role) ?? null;
 }
 
-export function computePlan(plan: PlanInput): PlanComputed {
+// Ein Zufalls-Szenario fuer die Monte-Carlo-Simulation: liefert je Jahr eine Inflation und
+// je Element/Jahr eine Rendite. Ohne Sample rechnet computePlan rein deterministisch (die
+// geplanten Annahmen), mit Sample einen einzelnen simulierten Pfad. Jahr ist 1-basiert
+// (ab Planbeginn); der Inflations-Index ist 0-basiert (inflation[0] = Jahr 1).
+export interface PlanSample {
+  inflation: number[];
+  assetReturn: (elementId: string, year: number) => number;
+}
+
+export function computePlan(plan: PlanInput, sample?: PlanSample): PlanComputed {
   const phases = [...plan.phases].sort((a, b) => a.sequenceNumber - b.sequenceNumber);
   const persons = plan.persons;
   const personA = persons.find((p) => p.role === "PERSON_A") ?? persons[0];
@@ -218,6 +227,15 @@ export function computePlan(plan: PlanInput): PlanComputed {
   const retirementAge = new Map<string, number>();
   for (const p of persons) retirementAge.set(p.id, p.retirementAge);
 
+  // Kumulierter Inflations-Deflator je Jahr (cumInfl[0] = 1, cumInfl[k] = Kaufkraftfaktor nach
+  // k Jahren). Deterministisch identisch zur bisherigen (1+infl)^k-Formel; mit Sample variiert
+  // die Inflation pro Jahr. Ersetzt die frueheren geschlossenen Potenz-Ausdruecke.
+  const totalYears = phases.reduce((s, p) => s + Math.max(1, p.durationYears), 0);
+  const inflationOfYear = (year: number) =>
+    sample ? sample.inflation[year - 1] ?? plan.inflationRateDefault : plan.inflationRateDefault;
+  const cumInfl: number[] = [1];
+  for (let y = 1; y <= totalYears; y++) cumInfl[y] = cumInfl[y - 1] * (1 + inflationOfYear(y) / 100);
+
   const gapYearsByPerson = new Map<string, number>();
   // AHV-Beitragskarriere je Person: reales Einkommen x Beitragsjahre, und Beitragsjahre.
   const ahvIncomeAccum = new Map<string, number>();
@@ -247,9 +265,7 @@ export function computePlan(plan: PlanInput): PlanComputed {
     const nextPhase = phases[i + 1];
     const isFirstPhase = i === 0;
     const duration = Math.max(1, phase.durationYears);
-    // Inflation ist plan-weit (V5): keine Phasen-Ueberschreibung mehr.
-    const infl = plan.inflationRateDefault;
-    const cumInflStart = cumulativeInflation; // Kaufkraft-Deflator zu Phasenbeginn
+    const cumInflStart = cumInfl[yearsBefore]; // Kaufkraft-Deflator zu Phasenbeginn
 
     const personInfos: PersonPhaseInfo[] = persons.map((p) => {
       const ra = retirementAge.get(p.id)!;
@@ -404,7 +420,10 @@ export function computePlan(plan: PlanInput): PlanComputed {
             const attributed =
               owner ?? (plan.householdType === "SINGLE" && e.ownerRole === "HOUSEHOLD" ? personA : null);
             if (attributed) {
-              const avgRealGross = avgRealFlow(basis, idx, infl, duration, cumInflStart) * AHV_GROSS_FROM_NET_FACTOR;
+              // Die AHV-Karriere ist eine Real-Groesse auf Planungsbasis -- bewusst mit der
+              // festen Plan-Inflation, nicht der (evtl. gewuerfelten) Sample-Inflation.
+              const avgRealGross =
+                avgRealFlow(basis, idx, plan.inflationRateDefault, duration, cumInflStart) * AHV_GROSS_FROM_NET_FACTOR;
               phaseRealIncomeByPerson.set(
                 attributed.id,
                 (phaseRealIncomeByPerson.get(attributed.id) ?? 0) + avgRealGross
@@ -553,7 +572,7 @@ export function computePlan(plan: PlanInput): PlanComputed {
       let incomeFlow = renteTotal;
       for (const inc of incomes) incomeFlow += inc.basis * Math.pow(1 + inc.idx / 100, t - 1);
       // Ausgaben: real (Basis x (1+reale Mehrausgabe)^(t-1)); nominal = real x kumul. Inflation.
-      const inflFactor = cumInflStart * Math.pow(1 + infl / 100, t - 1);
+      const inflFactor = cumInfl[yearsBefore + t - 1];
       let expenseRealBase = 0;
       for (const exp of expenses) expenseRealBase += exp.basis * Math.pow(1 + exp.idx / 100, t - 1);
 
@@ -592,7 +611,8 @@ export function computePlan(plan: PlanInput): PlanComputed {
       // Vermoegen verzinsen + Sparbeitrag; Bezugsrate entnehmen (gekappt am Bestand) und ins Cash.
       let cashFromWithdraw = 0;
       for (const a of assets) {
-        const grown = a.value * (1 + a.r / 100) + a.rate;
+        const r = sample ? sample.assetReturn(a.ec.elementId, yearsBefore + t) : a.r;
+        const grown = a.value * (1 + r / 100) + a.rate;
         const w = Math.min(a.withdrawal, Math.max(0, grown));
         a.value = grown - w;
         cashFromWithdraw += w;
@@ -607,7 +627,8 @@ export function computePlan(plan: PlanInput): PlanComputed {
         debtRates += pay;
         // Wertsteigerung wirkt auf die LIEGENSCHAFT, nicht auf das Eigenkapital -- das ist der
         // Hebel: 1 % von 1 Mio sind 10'000, also 10 % eines Eigenkapitals von 100'000.
-        re.value *= 1 + re.growth / 100;
+        const g = sample ? sample.assetReturn(re.ec.elementId, yearsBefore + t) : re.growth;
+        re.value *= 1 + g / 100;
       }
       for (const d of debts) {
         const pay = Math.min(d.repay, d.owed);
@@ -629,7 +650,7 @@ export function computePlan(plan: PlanInput): PlanComputed {
 
     // Flow-Deflator fuer den Endwert (Jahr `duration`): eine Kaufkraft-Stufe weniger als der
     // Bestands-Deflator am Phasenende.
-    const flowDeflatorEnd = cumInflStart * Math.pow(1 + infl / 100, duration - 1);
+    const flowDeflatorEnd = cumInfl[yearsBefore + duration - 1];
 
     // Endwerte je Element setzen (Einkommen/Ausgaben nominal; Ausgaben-Nominal = real x Infl.).
     for (const inc of incomes) {
@@ -663,7 +684,7 @@ export function computePlan(plan: PlanInput): PlanComputed {
     const cashEnd = Math.round(cash);
     const startWealthNominal = Math.round(wealthStart + cashStart);
     const endWealthNominal = Math.round(wealthEnd + cashEnd);
-    cumulativeInflation = cumInflStart * Math.pow(1 + infl / 100, duration);
+    cumulativeInflation = cumInfl[yearsBefore + duration];
 
     result.push({
       id: phase.id,

src/lib/montecarlo.test.ts

@@ -0,0 +1,107 @@
+import { describe, it, expect } from "vitest";
+import { computePlan } from "@/lib/calculations";
+import { runMonteCarlo, defaultVolatilityLevel, RETURN_VOLATILITY_LEVELS, type MonteCarloParams } from "@/lib/montecarlo";
+import type { PlanInput } from "@/lib/types";
+
+// Plan: 40-jaehrig, 1 Phase 10 Jahre, ein Sonstiges Vermoegen 100'000 @ 5 %, 2 % Inflation.
+// Deterministisch: 100'000 x 1.05^10 = 162'889 Endvermoegen, kein Ruin.
+function basePlan(): PlanInput {
+  return {
+    id: "p", name: "T", householdType: "SINGLE", inflationRateDefault: 2, initialCash: 0,
+    persons: [{ id: "A", role: "PERSON_A", name: null, age: 40, retirementAge: 70 }],
+    phases: [{ id: "p1", sequenceNumber: 1, name: "p1", durationYears: 10, cashTransition: {} }],
+    elements: [
+      {
+        id: "asset", category: "OTHER_ASSET", name: "ETF", ownerRole: "HOUSEHOLD", orderIndex: 1,
+        phaseValues: { p1: { startValue: 100000, expectedReturn: 5, annualContribution: 0 } },
+        transitionValues: {},
+      },
+    ],
+  };
+}
+
+function params(over: Partial<MonteCarloParams> = {}): MonteCarloParams {
+  return {
+    runs: 2000,
+    inflationMean: 2,
+    inflationSigma: 0,
+    elements: { asset: { mean: 5, sigma: 0, floor: -100 } },
+    target: 0,
+    seed: 12345,
+    ...over,
+  };
+}
+
+const detEnd = () => computePlan(basePlan()).phases[0].endWealthNominal;
+
+describe("Monte Carlo", () => {
+  it("Streuung 0 reproduziert exakt das deterministische Ergebnis", async () => {
+    const r = await runMonteCarlo(basePlan(), params());
+    const det = detEnd();
+    expect(r.finalWealthMedian).toBe(det);
+    expect(r.finalWealthP10).toBe(det);
+    expect(r.finalWealthP90).toBe(det); // Baender kollabieren auf die deterministische Linie
+    expect(r.ruinProbability).toBe(0);
+    expect(r.bands[r.bands.length - 1].p50).toBe(det);
+  });
+
+  it("Volatilitaet spreizt den Faecher (p90 > p10)", async () => {
+    const r = await runMonteCarlo(basePlan(), params({ elements: { asset: { mean: 5, sigma: 20, floor: -100 } } }));
+    expect(r.finalWealthP90).toBeGreaterThan(r.finalWealthP10);
+    // Median bleibt in der Naehe des deterministischen Werts (leicht darunter wegen Vol-Drag).
+    expect(r.finalWealthMedian).toBeLessThan(r.finalWealthP90);
+    expect(r.finalWealthMedian).toBeGreaterThan(r.finalWealthP10);
+  });
+
+  it("Erfolgswahrscheinlichkeit sinkt mit steigendem Zielbetrag", async () => {
+    const opts = { elements: { asset: { mean: 5, sigma: 20, floor: -100 } } };
+    const low = await runMonteCarlo(basePlan(), params({ ...opts, target: 50000 }));
+    const high = await runMonteCarlo(basePlan(), params({ ...opts, target: 500000 }));
+    expect(low.successProbability).toBeGreaterThan(high.successProbability);
+    expect(low.successProbability).toBeGreaterThan(0.9); // 50k ist fast sicher erreicht
+    expect(high.successProbability).toBeLessThan(0.1); // 500k praktisch unerreichbar
+  });
+
+  it("gleicher Seed -> identisches Ergebnis (reproduzierbar)", async () => {
+    const opts = params({ elements: { asset: { mean: 5, sigma: 20, floor: -100 } }, seed: 777 });
+    const a = await runMonteCarlo(basePlan(), opts);
+    const b = await runMonteCarlo(basePlan(), opts);
+    expect(a.finalWealthMedian).toBe(b.finalWealthMedian);
+    expect(a.ruinProbability).toBe(b.ruinProbability);
+  });
+
+  it("Boden 0 % (PK/3a): Rendite nie negativ -> Endvermoegen nie unter dem Startwert", async () => {
+    // Ohne Beitraege kann ein bei 0 % gebodetes Asset nur wachsen oder gleich bleiben.
+    const r = await runMonteCarlo(
+      basePlan(),
+      params({ elements: { asset: { mean: 0, sigma: 80, floor: 0 } } })
+    );
+    expect(r.finalWealthP10).toBeGreaterThanOrEqual(100000);
+  });
+
+  it("Ruinwahrscheinlichkeit: sicherer Verzehr fuehrt immer in den Ruin", async () => {
+    // Rente 20k, Ausgaben 60k, kleines Vermoegen -> deterministisch Ruin, auch ohne Streuung.
+    const plan: PlanInput = {
+      id: "p", name: "T", householdType: "SINGLE", inflationRateDefault: 0, initialCash: 0,
+      persons: [{ id: "A", role: "PERSON_A", name: null, age: 65, retirementAge: 65 }],
+      phases: [{ id: "p1", sequenceNumber: 1, name: "p1", durationYears: 20, cashTransition: {} }],
+      elements: [
+        { id: "inc", category: "INCOME", name: "Rente", ownerRole: "PERSON_A", orderIndex: 1, phaseValues: { p1: { amount: 20000, teuerungsausgleich: 0 } }, transitionValues: {} },
+        { id: "exp", category: "EXPENSE", name: "Ausgaben", ownerRole: "HOUSEHOLD", orderIndex: 2, phaseValues: { p1: { amount: 60000, teuerungsausgleich: 0 } }, transitionValues: {} },
+        { id: "asset", category: "OTHER_ASSET", name: "V", ownerRole: "HOUSEHOLD", orderIndex: 3, phaseValues: { p1: { startValue: 100000, expectedReturn: 0, annualContribution: 0 } }, transitionValues: {} },
+      ],
+    };
+    const r = await runMonteCarlo(plan, {
+      runs: 1000, inflationMean: 0, inflationSigma: 0,
+      elements: { asset: { mean: 0, sigma: 10, floor: -100 } }, target: 0, seed: 1,
+    });
+    expect(r.ruinProbability).toBe(1); // Verzehr uebersteigt Rente + Vermoegen in jedem Pfad
+  });
+
+  it("Hilfsfunktionen: Default-Stufen und σ-Tabelle", () => {
+    expect(defaultVolatilityLevel("PENSION_FUND")).toBe("sehr_niedrig");
+    expect(defaultVolatilityLevel("OTHER_ASSET")).toBe("moderat");
+    expect(RETURN_VOLATILITY_LEVELS.moderat).toBe(15);
+    expect(RETURN_VOLATILITY_LEVELS.sehr_hoch).toBe(55);
+  });
+});

src/lib/montecarlo.ts

@@ -0,0 +1,221 @@
+// Monte-Carlo-Simulation (Roadmap Nr. 19, Stufe A). Laeuft vollstaendig im Browser, weil
+// computePlan eine reine Funktion ohne Server-Abhaengigkeiten ist.
+//
+// Modell: Pro Jahr wird EIN gemeinsamer Marktschock gezogen (damit riskante Anlagen zusammen
+// fallen, nicht gegeneinander). Jede Rendite = historischer Mittelwert + Standardabweichung x
+// (Marktanteil x Marktschock + Eigenanteil x Eigenrauschen). Beide Schocks sind fettschwaenzig
+// (standardisierte Student-t), damit Extremcrashs realistisch haeufig auftreten -- eine
+// Normalverteilung wuerde sie stark unterschaetzen. Boeden: 0 % fuer PK/3a, -100 % sonst.
+
+import { computePlan, type PlanSample } from "@/lib/calculations";
+import type { ElementCategory } from "@/lib/elements";
+import type { PlanInput } from "@/lib/types";
+
+// Standardabweichungen (annualisiert, in %) hinter den Streuungsstufen. Recherchiert und
+// gerundet: Anleihen ~6 %, globale Aktien ~15-18 %, Schweizer Immobilien(fonds) ~2 %,
+// Bitcoin ~54 %. Quellen: siehe SPEZIFIKATION Kap. 4.
+export const RETURN_VOLATILITY_LEVELS = {
+  sehr_niedrig: 3,
+  niedrig: 6,
+  moderat: 15,
+  hoch: 25,
+  sehr_hoch: 55,
+} as const;
+
+// Inflation ist in der Schweiz historisch stabil (Standardabweichung der letzten 20 Jahre
+// ~1 %). Hoehere Stufen ergaeben Hyperinflations-Annahmen -- deshalb nur zwei Stufen.
+export const INFLATION_VOLATILITY_LEVELS = {
+  sehr_niedrig: 1,
+  niedrig: 2,
+} as const;
+
+export type ReturnVolatilityLevel = keyof typeof RETURN_VOLATILITY_LEVELS | "manuell";
+export type InflationVolatilityLevel = keyof typeof INFLATION_VOLATILITY_LEVELS | "manuell";
+
+// Default-Stufe je Element-Typ (damit die Simulation ohne Eingabe eine plausible Streuung hat).
+export function defaultVolatilityLevel(category: ElementCategory): ReturnVolatilityLevel {
+  switch (category) {
+    case "PENSION_FUND":
+      return "sehr_niedrig";
+    case "REAL_ESTATE":
+    case "PILLAR_3A":
+      return "niedrig";
+    default:
+      return "moderat";
+  }
+}
+
+// Kategorien mit Marktrendite -- nur diese brauchen Monte-Carlo-Parameter.
+export const RETURN_BEARING: ElementCategory[] = ["PENSION_FUND", "PILLAR_3A", "OTHER_ASSET", "REAL_ESTATE"];
+
+// PK/3a schreiben dem Versicherten keine negative Rendite gut -> Boden bei 0 %.
+export function floorFor(category: ElementCategory): number {
+  return category === "PENSION_FUND" || category === "PILLAR_3A" ? 0 : -100;
+}
+
+export interface ElementMcParams {
+  mean: number; // historische Durchschnittsrendite (%/Jahr)
+  sigma: number; // Standardabweichung (%/Jahr)
+  floor: number; // 0 fuer PK/3a, -100 sonst
+}
+
+export interface MonteCarloParams {
+  runs: number;
+  inflationMean: number;
+  inflationSigma: number;
+  elements: Record<string, ElementMcParams>; // key = elementId
+  target: number; // Zielbetrag fuer die Erfolgswahrscheinlichkeit (nominal)
+  seed?: number;
+}
+
+export interface MonteCarloResult {
+  runs: number;
+  ruinProbability: number; // P(Vermoegen faellt vor Planende unter 0)
+  successProbability: number; // P(Endvermoegen >= Zielbetrag)
+  finalWealthP10: number;
+  finalWealthMedian: number;
+  finalWealthP90: number;
+  // Faecher ueber das Alter: je Alterspunkt der pessimistische/mittlere/optimistische Wert.
+  bands: { age: number; p10: number; p50: number; p90: number }[];
+}
+
+// --- Zufallszahlen (seedbar, damit ein Lauf reproduzierbar ist) ---
+
+function mulberry32(seed: number): () => number {
+  let a = seed >>> 0;
+  return () => {
+    a |= 0;
+    a = (a + 0x6d2b79f5) | 0;
+    let t = Math.imul(a ^ (a >>> 15), 1 | a);
+    t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
+    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
+  };
+}
+
+function normal(rng: () => number): number {
+  // Box-Muller. Kleiner Schutz gegen log(0).
+  const u = Math.max(rng(), 1e-12);
+  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * rng());
+}
+
+// Standardisierte Student-t mit nu Freiheitsgraden (Einheitsvarianz -> die eingegebene
+// Standardabweichung bleibt die tatsaechliche). Fette Raender: nu = 5.
+const NU = 5;
+function studentT(rng: () => number): number {
+  let chi2 = 0;
+  for (let i = 0; i < NU; i++) {
+    const z = normal(rng);
+    chi2 += z * z;
+  }
+  const t = normal(rng) / Math.sqrt(chi2 / NU);
+  return t * Math.sqrt((NU - 2) / NU); // auf Einheitsvarianz standardisieren
+}
+
+// Ladefaktor auf den gemeinsamen Marktschock: rho^2 ist die Korrelation zweier riskanter
+// Anlagen. rho = 0.7 -> ~0.5. Diversifikation hilft etwas, rettet aber nicht im Crash.
+const RHO = 0.7;
+const RHO_IDIO = Math.sqrt(1 - RHO * RHO);
+
+function percentile(sortedAsc: number[], p: number): number {
+  if (sortedAsc.length === 0) return 0;
+  const idx = Math.min(sortedAsc.length - 1, Math.max(0, Math.round(p * (sortedAsc.length - 1))));
+  return sortedAsc[idx];
+}
+
+// Ein Alterspunkt je Phasenanfang plus das Planende (wie im Vermoegensverlauf-Chart).
+function agePointsOf(plan: PlanInput): number[] {
+  const phases = [...plan.phases].sort((a, b) => a.sequenceNumber - b.sequenceNumber);
+  const startAge = plan.persons.find((p) => p.role === "PERSON_A")?.age ?? plan.persons[0]?.age ?? 0;
+  const points: number[] = [startAge];
+  let acc = 0;
+  for (const p of phases) {
+    acc += Math.max(1, p.durationYears);
+    points.push(startAge + acc);
+  }
+  return points;
+}
+
+// Wert je Alterspunkt fuer EINEN Durchlauf (Start-/Endvermoegen der Phasen, nominal).
+function wealthTrajectory(computed: ReturnType<typeof computePlan>): number[] {
+  const phases = computed.phases;
+  if (phases.length === 0) return [];
+  const pts = [phases[0].startWealthNominal];
+  for (const p of phases) pts.push(p.endWealthNominal);
+  return pts;
+}
+
+export async function runMonteCarlo(
+  plan: PlanInput,
+  params: MonteCarloParams,
+  onProgress?: (done: number, total: number) => void
+): Promise<MonteCarloResult> {
+  const phases = [...plan.phases].sort((a, b) => a.sequenceNumber - b.sequenceNumber);
+  const totalYears = phases.reduce((s, p) => s + Math.max(1, p.durationYears), 0);
+  const rng = mulberry32(params.seed ?? (Math.random() * 2 ** 32) >>> 0);
+
+  const agePoints = agePointsOf(plan);
+  const nPoints = agePoints.length;
+  const wealthByPoint: number[][] = Array.from({ length: nPoints }, () => []);
+  const finalWealth: number[] = [];
+  let ruinCount = 0;
+  let successCount = 0;
+
+  const elementIds = Object.keys(params.elements);
+
+  for (let run = 0; run < params.runs; run++) {
+    // Pro Jahr ein Marktschock; pro Element/Jahr eine Rendite (fette Raender, gebodet).
+    const inflation: number[] = new Array(totalYears);
+    for (let y = 0; y < totalYears; y++) {
+      inflation[y] = params.inflationMean + params.inflationSigma * studentT(rng);
+    }
+
+    const returns = new Map<string, Float64Array>();
+    for (const id of elementIds) returns.set(id, new Float64Array(totalYears + 1));
+    for (let y = 1; y <= totalYears; y++) {
+      const market = studentT(rng);
+      for (const id of elementIds) {
+        const ep = params.elements[id];
+        const shock = RHO * market + RHO_IDIO * studentT(rng);
+        const r = ep.mean + ep.sigma * shock;
+        returns.get(id)![y] = Math.max(ep.floor, r);
+      }
+    }
+
+    const sample: PlanSample = {
+      inflation,
+      assetReturn: (elementId, year) => returns.get(elementId)?.[year] ?? 0,
+    };
+
+    const computed = computePlan(plan, sample);
+    if (computed.ruinAge !== null) ruinCount++;
+
+    const traj = wealthTrajectory(computed);
+    for (let i = 0; i < nPoints; i++) wealthByPoint[i].push(traj[i] ?? 0);
+    const fw = traj[traj.length - 1] ?? 0;
+    finalWealth.push(fw);
+    if (fw >= params.target) successCount++;
+
+    // Alle ~500 Laeufe die Kontrolle abgeben, damit die Oberflaeche nicht einfriert.
+    if (run % 500 === 499) {
+      onProgress?.(run + 1, params.runs);
+      await new Promise((r) => setTimeout(r, 0));
+    }
+  }
+  onProgress?.(params.runs, params.runs);
+
+  finalWealth.sort((a, b) => a - b);
+  const bands = agePoints.map((age, i) => {
+    const col = wealthByPoint[i].sort((a, b) => a - b);
+    return { age, p10: percentile(col, 0.1), p50: percentile(col, 0.5), p90: percentile(col, 0.9) };
+  });
+
+  return {
+    runs: params.runs,
+    ruinProbability: ruinCount / params.runs,
+    successProbability: successCount / params.runs,
+    finalWealthP10: percentile(finalWealth, 0.1),
+    finalWealthMedian: percentile(finalWealth, 0.5),
+    finalWealthP90: percentile(finalWealth, 0.9),
+    bands,
+  };
+}