FIREを成功させる資産配分:モンテカルロシミュレーションで別世界5000個生成し検証

第1章：オルカン1本の夢と現実

“できるだけ働かずに生きたい” という欲望は皆さんありますか？

…そこで気になるのが、よく巷で言われる投資論。

例えば、こんなシナリオを見たことはありませんか？

• 5,000万円をオルカンに一括投資

• 年平均リターン8.5%で30年運用

• 複利の力で6億円になる！

まさに夢のような数字です。
私も最初はこの「ウハウハ説」を信じかけました。

平均リターンを鵜呑みにする危うさ

しかし、実際に生活費や税金、社会保険料を差し引きながら資産を取り崩すとどうなるでしょうか？
FIREシミュレーションをしてみると、意外な落とし穴が見えてきました。

仮定条件はこうです：

• 初期資産：7,000万円

• 取り崩し：年間4%（インフレ調整あり）

• 期間：40年

• 分布：暴落を強めに反映するため Student-t 分布（df=5）

• 成功条件：資産が一度も4,500万円（-35%のドローダウン相当）を下回らないこと

結果は驚きでした。

シミュレーション結果：オルカン100%の現実

• オルカン100%：成功率は 約50%台。

• つまり、2人に1人は途中で資産が底を割る計算になります。

中央値では数億円に育つ可能性もありますが、谷底に落ちるリスクは無視できません。

読者への問いかけ

「もしあなたが今の資産をオルカン100%に投じて、
40年間取り崩しながら生活するとしたら——
2人に1人は破産するコイン投げに人生を賭けられますか？」

第2章：株だけでは心許ない

前章で「オルカン100%はコイン投げ」という現実を見ました。
では、株式の中でも人気の QQQ（NASDAQ100） や VOO（S&P500） だけに賭けたらどうでしょうか？

世間ではよく、

果たして本当にそうでしょうか？

シミュレーション条件

• 初期資産：7,000万円

• 取り崩し：年間4%（FR4 / PW4 両ルール）

• 期間：40年

• 分布：Student-t(df=5)

• 比較対象：

  • QQQ 100%

  • VOO 100%

  • 参考：オルカン100%

結果1：成功率の比較

• QQQ100%
  成功率はオルカンよりやや上。ただし上下のブレが極端。

• VOO100%
  オルカンと似た傾向。米国株に偏る分、分散効果は薄い。

• オルカン100%
  やはり50%台のコイン投げ水準。

結果2：期末資産の分布

• QQQ100%：伸び（億超えの夢）は大きい。しかし同時に、資産が底を割るリスクも高い。

• VOO100%：オルカンと大差なく、中庸の結果。

• オルカン100%：全世界分散でも、谷底リスクは完全には消えない。

読者への問いかけ

「夢を追うならQQQ100%。
ただし、その夢と引き換えに“生活資金が尽きる”可能性も高まる。

あなたなら、夢の爆発力と破産リスクのどちらを取りますか？」

第3章：GLDとEDVを加える

前章では「株だけだと夢はあるが谷底も深い」という現実を見ました。
では、株と相関の低い資産を組み合わせたらどうなるのでしょうか？

投資の世界でよく耳にするのが「分散投資」。
特に「株と逆に動く資産を持て」と言われます。

候補に挙がるのが、ゴールド（GLD）と米国超長期国債（EDV）です。

なぜGLDとEDVなのか？

• GLD（ゴールド）

  • 株式との相関が低く、独立した値動きをする。

  • インフレ局面でも価値が保たれることが多い。

• EDV（米国超長期国債）

  • 株式と弱い逆相関（-0.1程度）。

  • 株が暴落するときに買われやすく「床」を支える役割。

この2つを組み合わせれば、株の弱点である「暴落時の取り崩し」を緩和できるはずです。

検証ポートフォリオ

• P1：QQQ 50% + GLD 50%

• P2：QQQ 50% + GLD 30% + EDV 20%

• P3：QQQ 40% + GLD 40% + EDV 10% + AGG 10%（米国総合債券も追加）

結果1：成功率

• P1：シンプルにGLDを半分入れただけでも、成功率は向上。

• P2：GLDとEDVを組み合わせるとさらに改善。成功率80%台へ。

• P3：AGGを加えることで変動がさらに抑えられ、安定度が最も高い。

結果2：期末資産分布

• P1：右方向（成長の夢）はやや削れるが、下側が浅くなる。

• P2：中央値も改善し、左の尻尾が短く「破産ゾーン」を避けやすい。

• P3：最も下振れが浅く、“資産寿命を延ばす” ポートフォリオ。

効率的フロンティアから見ても明らか

色分けを見ると、株100%の右上だけでなく、GLDや債券を組み合わせたゾーンが鮮やかに光る。
つまり「リスクを抑えつつリターン効率を高める」効果が、データ上も証明されました。

ちなみに、2つの効率的フロンティアのポートフォリオが算出できたので、ここに一応記載いたします。

読者への問いかけ

「あなたは“億超えの夢”を追って破産リスクを抱えますか？
それとも、“そこそこ増えて破産しない”現実的な道を選びますか？」

第4章：効率的フロンティアで見えてきたこと

前章では「GLDやEDVを混ぜると破産確率が大きく下がる」ことを確認しました。
では実際に、どんな比率で組み合わせるのが効率的なのでしょうか？

投資の世界には「効率的フロンティア」という考え方があります。
これは「リスクに対してリターンが最大になる資産配分の曲線」を描き出す手法です。

効率的フロンティアとは？

• 横軸：リスク（資産のブレ幅＝標準偏差）

• 縦軸：期待リターン

• 曲線の上にある点ほど「投資効率が良い」

つまり、同じリスクを取るならリターンが高い方が良いし、同じリターンならリスクが低い方が良い。
その「最も効率の良い組み合わせ」が並んだのが、効率的フロンティアです。

実際にExcelソルバーとPythonで算出

対象資産は、これまで登場した6種類：

• QQQ（NASDAQ100）

• VOO（S&P500）

• ACWI（オルカン）

• GLD（ゴールド）

• EDV（米国超長期国債）

• AGG（米国総合債券）

30,000通り以上のランダムポートフォリオを生成し、リスクリターンを算出しました。

結果：ナスダックと金をメインとしサブで債券の組み合わせ

• 右上（高リスク・高リターン）の端にはQQQ100%が位置。夢はあるが振れ幅も大きい。

• 真ん中あたりで鮮やかに色付いたのは、QQQ＋GLD+EDV＋AGGを含むポートフォリオ。

• 株の比率を下げつつゴールド・債券を組み合わせることで、Sharpe比率（リスク調整後リターン）が改善。

意外な事実

• QQQとVOOはほぼ同じ動き → 両方を持つ意味は薄い。

• GLDは株との相関が低い → 分散効果が非常に大きい。

• EDVは株と弱い逆相関 → 暴落時のクッション役。

• AGGは値動きが小さく、全体の安定化要因。

効率的フロンティアが示すのは、株100%が最強ではないという現実です。

読者への問いかけ

「あなたのポートフォリオは“なんとなく感覚”で組んでいませんか？
データに基づくと、リスクを抑えながら効率良く増やせる組み合わせは確かに存在します。

では、実際にその配分で40年シミュレーションしたらどうなるでしょう？」

第5章：最適ポートフォリオのモンテカルロ検証

効率的フロンティアで見えた「株＋ゴールド＋債券の組み合わせが効率的」という事実。
では実際に、その配分で40年生活をシミュレーションしたらどうなるのでしょうか？

ここからは、モンテカルロ・シミュレーションを使って「破産しにくく、それでも増やせる黄金比率」を探っていきます。

シミュレーション条件

• 初期資産：7,000万円

• 取り崩しルール

  • FR4：インフレ調整後、毎年固定で4%引き出す

  • PW4：その年の資産の4%を引き出す

• 期間：40年間

• 試行回数：5,000回

• 分布：Student-t分布（df=5）＝暴落を強めに反映

• 成功基準：資産が一度も4,500万円（-35%のドローダウン相当）を下回らないこと

比較ポートフォリオ

• ACWI100%（オルカン一本）

• QQQ100%（NASDAQ100）

• VOO100%（S&P500）

• P1：QQQ50% + GLD50%

• P2：QQQ50% + GLD30% + EDV20%

• P3：QQQ40% + GLD40% + EDV10% + AGG10%

結果1：成功率

• ACWI/VOO/QQQ100%：成功率は50〜60%台。いずれも「コイン投げ水準」。

• P1：GLDを半分入れるだけで成功率が上昇。

• P2：さらにEDVを加えると80%台に改善。

• P3：AGGを組み込むことで最も安定し、FR4ルールでも8割以上の成功率。

結果2：取り崩しルールを変えると？

• PW4（資産の4%ルール）に切り替えると、全体的に成功率がさらに改善。

• P2やP3では9割近い成功率を記録。

• 支出を資産に合わせて柔軟にするだけで、「資産寿命」が大きく延びる。

結果3：期末資産の分布

• QQQ100%：上振れは大きいが、下側の尻尾も深い。

• P2/P3：中央値が安定し、谷底に落ちる確率が大幅に低下。

• 特にP3は「左の尻尾が短い」＝破産確率が極めて低い。

考察

• QQQ100%は“夢”を見られるが“谷底”も深い。

• オルカン/VOO100%は無難だが、破産リスクを完全には消せない。

• GLD＋EDV＋AGGを組み合わせたP2/P3は、破産リスクを大幅に低下させつつ、資産を増やす現実的な選択肢。

効率的フロンティアとモンテカルロの両面から、同じ結論が浮かび上がりました。

読者への問いかけ

「あなたなら、億超えを夢見てオルカンやレバナスへの“コイン投げ”に賭けますか？
それとも、“破産しない現実的な配分”で安心を取りに行きますか？」

第6章：まとめと行動提案

ここまで、オルカン一本の夢から始まり、株だけのリスク、そしてゴールドや債券を組み合わせた分散投資までシミュレーションしてきました。
データが示したのはシンプルです。

これまでの検証で分かったこと

1. オルカン100%は夢があるが、現実はコイン投げ

   • 成功率は50%前後。2人に1人は途中で資産が底を割る。

2. QQQ100%は爆発力があるが“谷底”も深い

   • 億超えも夢じゃないが、破産リスクを大きく抱える。

3. ゴールド＋債券を加えると破産確率が激減

   • P2（QQQ50/GLD30/EDV20）やP3（QQQ40/GLD40/EDV10/AGG10）は成功率80〜90%に改善。

4. 取り崩しルールを工夫するだけで生存率がさらに上がる

   • 固定額（FR4）より、資産の割合で決める（PW4）の方が寿命が長い。

FIREを考える人への教訓

• 「平均リターン8.5%で6億円！」というのは幻想。

• 投資は未来を保証するものではなく、確率のゲーム。

• 大事なのは「どれだけ資産を増やすか」よりも、いかに破産しないか。

行動提案

あなたが今日からできることはシンプルです。

1. 小さく始める

   • まずは少額で、自分の資産でシミュレーションしてみる。

   • 「机上の理論」を「自分の数字」に落とし込むこと。

2. 区分を意識する

   • 生活費、余剰資金、投資資金を分ける。

   • 生活費までマーケットに賭けないことが安心につながる。

3. 複線を作る

   • 投資一本足打法ではなく、副業・年金・配当など複数のキャッシュフローを用意する。

   • FIREは「投資の勝ち負け」ではなく、「収入源を多様化したライフデザイン」。

最後に

夢を追うのは自由です。
けれど、データが示す現実を直視しなければ「働かない人生」が逆に「働かざるを得ない人生」へと転落します。

あなたは——
夢の爆発力を取りますか？ それとも破産回避の現実を取りますか？

選択は、いまこの記事を読んでいるあなたの手の中にあります。

 計算の仕組み

Pythonのこの記事より簡易的にしたサンプルコード

"""
Monte Carlo Retirement Simulation — Full Spec (Internal AI Hand-off)
Author: <your name or org>
Python: 3.10+
Libraries: numpy, pandas, matplotlib

Unit convention: x10k JPY (1 unit = 10,000 JPY)
"""

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter

# ==============================
# Global Settings (Spec-Compliant)
# ==============================
SEED = 20250829  # reproducible
rng = np.random.default_rng(SEED)

N_TRIALS = 30000
ASSETS = ["QQQ", "GLD", "VOO", "EDV", "BND"]
YEARS_A = 20   # Part A horizon
YEARS_B = 50   # Part B horizon
UNIT_NOTE = "x10k JPY"  # 1 unit = 10,000 JPY

# Axis number format with comma
def comma_fmt(x, pos=None):
    return f"{int(x):,}"
fmt_comma = FuncFormatter(comma_fmt)

# Y-axis clipping (0 to 4億円)
Y_MIN, Y_MAX = 0, 40000

# Output directory
OUTDIR = "outputs"
os.makedirs(OUTDIR, exist_ok=True)

# ==============================
# Asset Parameters
# ==============================
means = np.array([0.159, 0.078, 0.092, 0.039, 0.040])
stds  = np.array([0.325, 0.150, 0.250, 0.245, 0.092])

# Correlation matrix (symmetrize, SPD-fix if needed)
corr_raw = np.array([
    [ 1.00, 0.06, 0.94, -0.11, 0.23],
    [ 0.16, 1.00, 0.18,  0.08, 0.21],
    [ 0.92, 0.08, 1.00, -0.11, 0.22],
    [-0.12, 0.26,-0.12,  1.00, 0.33],
    [ 0.57, 0.37, 0.61,  0.70, 1.00],
], dtype=float)

corr = 0.5 * (corr_raw + corr_raw.T)
np.fill_diagonal(corr, 1.0)

# Higham (2002) nearest SPD
def nearest_spd(A: np.ndarray) -> np.ndarray:
    B = (A + A.T) / 2
    U, s, Vt = np.linalg.svd(B)
    H = Vt.T @ np.diag(s) @ Vt
    A2 = (B + H) / 2
    A3 = (A2 + A2.T) / 2

    def is_pd(X):
        try:
            np.linalg.cholesky(X)
            return True
        except np.linalg.LinAlgError:
            return False

    if is_pd(A3):
        return A3

    spacing = np.spacing(np.linalg.norm(A))
    I = np.eye(A.shape[0])
    k = 1
    while not is_pd(A3):
        min_eig = np.min(np.real(np.linalg.eigvals(A3)))
        A3 += I * (-min_eig * k**2 + spacing)
        k += 1
    return A3

# Ensure SPD
try:
    Lcorr = np.linalg.cholesky(corr)
except np.linalg.LinAlgError:
    corr = nearest_spd(corr)
    Lcorr = np.linalg.cholesky(corr)

cov = np.outer(stds, stds) * corr

# ==============================
# Heavy-tail Multivariate t (df=5)
# ==============================
NU = 5  # degrees of freedom

def draw_asset_returns(years: int, trials: int) -> np.ndarray:
    """
    Generate (years, trials, assets) returns using multivariate-t with df=NU,
    correlation structure 'corr' and per-asset means/stds.
    Variance matches specified stds via scaling sqrt((nu-2)/W).
    """
    Z = rng.standard_normal(size=(years, trials, len(ASSETS)))
    Y = Z @ Lcorr.T  # correlate across assets
    W = rng.chisquare(df=NU, size=(years, trials))
    scale = np.sqrt((NU - 2) / W)  # variance preserving
    Y_t = Y * scale[..., None]
    R = means[None, None, :] + Y_t * stds[None, None, :]
    # Physical bound: avoid wealth flips
    R = np.clip(R, -0.999, None)
    return R

# ==============================
# Inflation Processes
# ==============================
def draw_inflation(mean: float, sd: float, years: int, trials: int) -> np.ndarray:
    """Draw (years, trials) inflation with floor 0%."""
    infl = rng.normal(loc=mean, scale=sd, size=(years, trials))
    infl = np.clip(infl, 0.0, None)
    return infl

infl_A = draw_inflation(0.02,   0.01,  YEARS_A, N_TRIALS)
infl_B = draw_inflation(0.0394, 0.0131, YEARS_B, N_TRIALS)

def make_spending_schedule_B(infl_mat: np.ndarray,
                             years: int = YEARS_B,
                             base1: float = 150.0,
                             base2: float = 250.0) -> np.ndarray:
    """
    Spending is post-tax requirement.
    Year1 uses base1; each next year multiplies by (1 + last year's inflation).
    At Year6 (index=5) onward, base resets to base2 then continues inflation chaining.
    Returns (years, trials) of required spending in x10k JPY.
    """
    req = np.zeros_like(infl_mat)
    req[0, :] = base1
    for t in range(1, min(5, years)):
        req[t, :] = req[t-1, :] * (1.0 + infl_mat[t-1, :])
    if years > 5:
        req[5, :] = base2 * (1.0 + infl_mat[4, :])
        for t in range(6, years):
            req[t, :] = req[t-1, :] * (1.0 + infl_mat[t-1, :])
    return req

spend_B = make_spending_schedule_B(infl_B, YEARS_B, base1=150.0, base2=250.0)

# ==============================
# Portfolios
# ==============================
portfolios = {
    "Portfolio1 QQQ+GLD": np.array([0.50, 0.50, 0.00, 0.00, 0.00]),
    "Portfolio2 Mix":     np.array([0.45, 0.45, 0.00, 0.05, 0.05]),
    "VOO 100%":           np.array([0.00, 0.00, 1.00, 0.00, 0.00]),
}

# Single-asset corners for scatter overlay
single_asset_ports = {
    "QQQ 100%": np.array([1,0,0,0,0]),
    "GLD 100%": np.array([0,1,0,0,0]),
    "VOO 100% (corner)": np.array([0,0,1,0,0]),
    "EDV 100%": np.array([0,0,0,1,0]),
    "BND 100%": np.array([0,0,0,0,1]),
}

# ==============================
# Helpers
# ==============================
def portfolio_returns(asset_returns: np.ndarray, weights: np.ndarray) -> np.ndarray:
    """Combine asset returns -> portfolio returns (years, trials)."""
    return np.einsum('yta,a->yt', asset_returns, weights)

def stats_quantiles(arr: np.ndarray, quantiles) -> dict:
    """Return dict of quantiles with keys PXX (zero-padded)."""
    qs = np.quantile(arr, quantiles)
    out = {}
    for q, v in zip(quantiles, qs):
        pct = int(round(q*100))
        out[f"P{pct:02d}"] = float(v)
    return out

# ==============================
# Part A Simulation (Accumulation)
# ==============================
ASSET_RET_A = draw_asset_returns(YEARS_A, N_TRIALS)
TARGET_A = 8000.0
INIT_TOTAL_A = 3000.0
INIT_INVEST_A = 2300.0
INIT_CASH_A = 700.0
FLOOR_A = 0.65 * INIT_TOTAL_A  # 65% floor

def simulate_part_a(weights: np.ndarray, name: str) -> dict:
    r_port = portfolio_returns(ASSET_RET_A, weights)  # years x trials
    invest_paths = np.empty_like(r_port)
    invest_paths[0, :] = INIT_INVEST_A * (1.0 + r_port[0, :])
    for t in range(1, YEARS_A):
        invest_paths[t, :] = invest_paths[t-1, :] * (1.0 + r_port[t, :])
    total_paths = invest_paths + INIT_CASH_A
    final_total = total_paths[-1, :]

    min_over_time = np.min(total_paths, axis=0)
    fail_floor = (min_over_time < FLOOR_A)
    hit_target = (final_total >= TARGET_A)

    qvals = stats_quantiles(final_total, [0.70, 0.50, 0.15, 0.05])
    return {
        "name": name,
        "final_total": final_total,     # (trials,)
        "total_paths": total_paths,     # (years, trials)
        "hit_target": hit_target,       # (trials,)
        "fail_floor": fail_floor,       # (trials,)
        "hit_prob": float(hit_target.mean()),
        "floor_fail_rate": float(fail_floor.mean()),
        "qvals": qvals,
    }

# ==============================
# Part B Simulation (Decumulation with Taxes)
# ==============================
ASSET_RET_B = draw_asset_returns(YEARS_B, N_TRIALS)

def simulate_part_b(weights: np.ndarray, name: str,
                    start_wealth: float = 8000.0,
                    start_cost_basis: float | None = None,
                    tax_rate: float = 0.2542) -> dict:
    """
    Decumulation with tax-aware withdrawals (Part B tax rate 25.42%).
    Cost basis defaults to start wealth (rebased at B start).
    """
    if start_cost_basis is None:
        start_cost_basis = start_wealth

    r_port = portfolio_returns(ASSET_RET_B, weights)  # years x trials
    wealth = np.full((YEARS_B, N_TRIALS), np.nan, dtype=float)
    basis  = np.full((YEARS_B, N_TRIALS), np.nan, dtype=float)

    cur_w = np.full(N_TRIALS, start_wealth, dtype=float)
    cur_b = np.full(N_TRIALS, start_cost_basis, dtype=float)
    floor_level = 0.65 * start_wealth

    for t in range(YEARS_B):
        # Apply return
        cur_w = cur_w * (1.0 + r_port[t, :])

        # Effective tax rate from unrealized gain ratio
        eps = 1e-9
        cur_b_safe = np.maximum(cur_b, eps)
        gain_ratio = np.maximum((cur_w - cur_b_safe) / cur_b_safe, 0.0)
        eff_tax = (gain_ratio / (1.0 + gain_ratio)) * tax_rate  # <= tax_rate

        # Required net spending for year t (already last-year inflation linked)
        net_need = spend_B[t, :]

        # Gross sale to cover net_need after tax
        gross_sale = net_need / np.maximum(1.0 - eff_tax, eps)
        gross_sale = np.minimum(gross_sale, cur_w)  # cap at wealth
        tax_paid = gross_sale * eff_tax
        principal_reduction = gross_sale - tax_paid  # reduces basis

        # Update wealth and basis
        cur_w = cur_w - gross_sale
        cur_b = np.maximum(cur_b - principal_reduction, 0.0)

        wealth[t, :] = cur_w
        basis[t, :]  = cur_b

    min_w = np.min(wealth, axis=0)
    success = (min_w >= floor_level)
    final_w = wealth[-1, :]

    p_under_principal = float((final_w < start_wealth).mean())
    p_over_principal  = float((final_w >= start_wealth).mean())

    qvals = {
        "P70": float(np.quantile(final_w, 0.70)),
        "P50": float(np.quantile(final_w, 0.50)),
        "P35": float(np.quantile(final_w, 0.35)),
        "P20": float(np.quantile(final_w, 0.20)),
        "P15": float(np.quantile(final_w, 0.15)),
        "P10": float(np.quantile(final_w, 0.10)),
        "P05": float(np.quantile(final_w, 0.05)),
    }

    return {
        "name": name,
        "wealth_paths": wealth,
        "basis_paths": basis,
        "success": success,
        "success_rate": float(success.mean()),
        "final_w": final_w,
        "p_under_principal": p_under_principal,
        "p_over_principal": p_over_principal,
        "qvals": qvals,
        "start_wealth": start_wealth,
    }

# ==============================
# Linked (A -> B)
# ==============================
def simulate_linked(name: str, weights: np.ndarray,
                    Ares: dict) -> dict:
    hit_idx = np.where(Ares["hit_target"])[0]
    if hit_idx.size == 0:
        return {
            "name": name,
            "A_hit_prob": float(Ares["hit_prob"]),
            "B_success_in_hit": np.nan,
            "linked_success": 0.0,
            "hit_count": 0,
            "indices": hit_idx,
            "wealth_paths": None,
        }

    # Start B with actual A final wealth, basis rebased to start
    start_wealths = Ares["final_total"][hit_idx]
    r_port = portfolio_returns(ASSET_RET_B, weights)[:, hit_idx]  # years x selected
    req = spend_B[:, hit_idx]

    T = YEARS_B
    M = hit_idx.size
    wealth = np.full((T, M), np.nan, dtype=float)
    basis  = np.full((T, M), np.nan, dtype=float)

    cur_w = start_wealths.copy()
    cur_b = start_wealths.copy()
    floor_level = 0.65 * start_wealths  # per-trial floor based on each start

    for t in range(T):
        cur_w = cur_w * (1.0 + r_port[t, :])
        eps = 1e-9
        cur_b_safe = np.maximum(cur_b, eps)
        gain_ratio = np.maximum((cur_w - cur_b_safe) / cur_b_safe, 0.0)
        tax_rate = 0.2542
        eff_tax = (gain_ratio / (1.0 + gain_ratio)) * tax_rate
        net_need = req[t, :]
        gross_sale = net_need / np.maximum(1.0 - eff_tax, eps)
        gross_sale = np.minimum(gross_sale, cur_w)
        tax_paid = gross_sale * eff_tax
        principal_reduction = gross_sale - tax_paid
        cur_w = cur_w - gross_sale
        cur_b = np.maximum(cur_b - principal_reduction, 0.0)
        wealth[t, :] = cur_w
        basis[t, :]  = cur_b

    min_w = np.min(wealth, axis=0)
    success = (min_w >= floor_level)
    B_success_rate_in_hit = float(success.mean())
    linked_success = float(Ares["hit_prob"] * B_success_rate_in_hit)

    return {
        "name": name,
        "A_hit_prob": float(Ares["hit_prob"]),
        "B_success_in_hit": B_success_rate_in_hit,
        "linked_success": linked_success,
        "hit_count": int(M),
        "indices": hit_idx,
        "wealth_paths": wealth,
    }

# ==============================
# Visualization Helpers
# ==============================
def save_fig(fig, filename: str):
    fig.tight_layout()
    path = os.path.join(OUTDIR, filename)
    fig.savefig(path, dpi=150, bbox_inches='tight')
    plt.close(fig)
    return path

def stat_box_text_A(res: dict) -> str:
    q = res["qvals"]
    text = (
        f"Hit Prob (>=80M): {res['hit_prob']*100:.2f}%\n"
        f"P70: {q['P70']:,.0f}  {UNIT_NOTE}\n"
        f"P50: {q['P50']:,.0f}  {UNIT_NOTE}\n"
        f"P15: {q['P15']:,.0f}  {UNIT_NOTE}\n"
        f"P05: {q['P05']:,.0f}  {UNIT_NOTE}"
    )
    return text

def stat_box_text_B(res: dict) -> str:
    q = res["qvals"]
    text = (
        f"Success Rate: {res['success_rate']*100:.2f}%\n"
        f"Final P70: {q['P70']:,.0f}\n"
        f"Final P50: {q['P50']:,.0f}\n"
        f"Final P15: {q['P15']:,.0f}\n"
        f"Final P05: {q['P05']:,.0f}  ({UNIT_NOTE})"
    )
    return text

# ==============================
# Plots per Spec
# ==============================
def plot_partA_hist(name: str, res: dict):
    fig = plt.figure(figsize=(8,5))
    plt.hist(res["final_total"], bins=60, alpha=0.75)
    plt.axvline(TARGET_A, linestyle='--', linewidth=1.5, label="Target 80M")
    plt.gca().xaxis.set_major_formatter(fmt_comma)
    plt.title(f"Part A Final Wealth Distribution — {name} (unit: {UNIT_NOTE})")
    plt.xlabel(f"Final Wealth ({UNIT_NOTE})")
    plt.ylabel("Count")
    plt.xlim(Y_MIN, Y_MAX)
    plt.legend(loc='lower right')
    plt.text(0.98, 0.98, stat_box_text_A(res), transform=plt.gca().transAxes,
             ha='right', va='top', bbox=dict(boxstyle="round", facecolor="white", alpha=0.8), fontsize=9)
    return save_fig(fig, f"PartA_hist_{name.replace(' ','_')}.png")

def plot_partA_paths(name: str, res: dict):
    paths = res["total_paths"]
    years_axis = np.arange(1, YEARS_A+1)
    subset = rng.choice(N_TRIALS, size=5000, replace=False)
    sample = rng.choice(subset, size=150, replace=False)

    fig = plt.figure(figsize=(8,5))
    for idx in sample:
        plt.plot(years_axis, np.clip(paths[:, idx], Y_MIN, Y_MAX), linewidth=0.5, alpha=0.3)
    median_path = np.median(paths, axis=1)
    plt.plot(years_axis, np.clip(median_path, Y_MIN, Y_MAX), linestyle='--', linewidth=2, label="Median")
    plt.axhline(INIT_TOTAL_A, color='black', linewidth=1.5, label="Initial 30M")
    plt.axhline(TARGET_A, linestyle='--', linewidth=1.5, label="Target 80M")
    plt.ylim(Y_MIN, Y_MAX)
    plt.gca().yaxis.set_major_formatter(fmt_comma)
    plt.title(f"Part A Wealth Paths (150 of 5,000) — {name}")
    plt.xlabel("Year")
    plt.ylabel(f"Wealth ({UNIT_NOTE})")
    plt.legend(loc='lower right')
    plt.text(0.98, 0.98, stat_box_text_A(res), transform=plt.gca().transAxes,
             ha='right', va='top', bbox=dict(boxstyle="round", facecolor="white", alpha=0.8), fontsize=9)
    return save_fig(fig, f"PartA_paths_{name.replace(' ','_')}.png")

def plot_partB_paths(name: str, res: dict):
    wealth = res["wealth_paths"]
    years_axis = np.arange(1, YEARS_B+1)
    subset = rng.choice(N_TRIALS, size=5000, replace=False)
    sample = rng.choice(subset, size=150, replace=False)

    fig = plt.figure(figsize=(8,5))
    for idx in sample:
        plt.plot(years_axis, np.clip(wealth[:, idx], Y_MIN, Y_MAX), linewidth=0.5, alpha=0.3)
    median_path = np.median(wealth, axis=1)
    plt.plot(years_axis, np.clip(median_path, Y_MIN, Y_MAX), linestyle='--', linewidth=2, label="Median")
    plt.axhline(8000, color='black', linewidth=1.5, label="Initial 80M")
    plt.ylim(Y_MIN, Y_MAX)
    plt.gca().yaxis.set_major_formatter(fmt_comma)
    plt.title(f"Part B Withdrawal Paths (150 of 5,000) — {name}")
    plt.xlabel("Year")
    plt.ylabel(f"Wealth ({UNIT_NOTE})")
    plt.legend(loc='lower right')
    plt.text(0.98, 0.98, stat_box_text_B(res), transform=plt.gca().transAxes,
             ha='right', va='top', bbox=dict(boxstyle="round", facecolor="white", alpha=0.8), fontsize=9)
    return save_fig(fig, f"PartB_paths_{name.replace(' ','_')}.png")

def plot_linked_paths(name: str, link: dict):
    if link["hit_count"] <= 0 or link["wealth_paths"] is None:
        return None
    wealth = link["wealth_paths"]
    years_axis = np.arange(1, YEARS_B+1)
    m = wealth.shape[1]
    k = min(150, m)
    sample = rng.choice(m, size=k, replace=False)

    fig = plt.figure(figsize=(8,5))
    for idx in sample:
        plt.plot(years_axis, np.clip(wealth[:, idx], Y_MIN, Y_MAX), linewidth=0.5, alpha=0.3)
    median_path = np.median(wealth, axis=1)
    plt.plot(years_axis, np.clip(median_path, Y_MIN, Y_MAX), linestyle='--', linewidth=2, label="Median")
    plt.ylim(Y_MIN, Y_MAX)
    plt.gca().yaxis.set_major_formatter(fmt_comma)
    plt.title(f"Linked B Paths (A-hit only) — {name}")
    plt.xlabel("Year")
    plt.ylabel(f"Wealth ({UNIT_NOTE})")
    plt.legend(loc='lower right')
    text = (f"A Hit Prob: {link['A_hit_prob']*100:.2f}%\n"
            f"B Success | Hit: {link['B_success_in_hit']*100:.2f}%\n"
            f"Linked Success: {link['linked_success']*100:.2f}%")
    plt.text(0.98, 0.98, text, transform=plt.gca().transAxes,
             ha='right', va='top', bbox=dict(boxstyle="round", facecolor="white", alpha=0.8), fontsize=9)
    return save_fig(fig, f"LinkedB_paths_{name.replace(' ','_')}.png")

def plot_success_rate_bars(results_A: dict, results_B: dict):
    labels = list(portfolios.keys())
    x = np.arange(len(labels))
    A_rates = [results_A[n]["hit_prob"]*100 for n in labels]
    B_rates = [results_B[n]["success_rate"]*100 for n in labels]

    fig = plt.figure(figsize=(8,5))
    width = 0.35
    plt.bar(x - width/2, A_rates, width=width, label="Part A Target-Hit (%)")
    plt.bar(x + width/2, B_rates, width=width, label="Part B Success (%)")
    plt.xticks(x, labels, rotation=0)
    plt.ylabel("Rate (%)")
    plt.title("Success Rate Comparison by Portfolio")
    plt.legend(loc='lower right')
    return save_fig(fig, "SuccessRate_Comparison.png")

def plot_B_terminal_overlay(results_B: dict):
    labels = list(portfolios.keys())
    fig = plt.figure(figsize=(8,5))
    bins = np.linspace(Y_MIN, Y_MAX, 80)
    for name in labels:
        arr = np.clip(results_B[name]["final_w"], Y_MIN, Y_MAX)
        plt.hist(arr, bins=bins, alpha=0.4, label=name)
    plt.gca().xaxis.set_major_formatter(fmt_comma)
    plt.title(f"Part B Final Wealth Distributions (unit: {UNIT_NOTE})")
    plt.xlabel(f"Final Wealth ({UNIT_NOTE})")
    plt.ylabel("Count")
    plt.xlim(Y_MIN, Y_MAX)
    plt.legend(loc='lower right')
    return save_fig(fig, "PartB_Final_Distributions.png")

def plot_risk_return_map(portfolios: dict):
    # Random 30,000 portfolios
    N_PORTS = 30000
    W = rng.dirichlet(alpha=np.ones(len(ASSETS)), size=N_PORTS)
    exp_ret = W @ means  # arithmetic
    exp_vol = np.sqrt(np.einsum('ij,jk,ik->i', W, cov, W))

    # Required CAGR for 30M -> 80M in 20y
    required_cagr = (TARGET_A / INIT_TOTAL_A) ** (1/YEARS_A) - 1
    score = exp_ret / required_cagr

    fig = plt.figure(figsize=(7.2,5.4))
    sc = plt.scatter(exp_vol, exp_ret, c=score, s=8, alpha=0.6)
    plt.colorbar(sc, label="Return / Required CAGR")
    plt.xlabel("Expected Volatility (σ)")
    plt.ylabel("Expected Return (μ)")
    plt.title("Random 30,000 Portfolios — Risk/Return Map")

    # Overlay named portfolios
    for name, w in portfolios.items():
        mu = float(w @ means)
        sig = float(np.sqrt(w @ cov @ w))
        plt.scatter([sig], [mu], marker='*', s=180, edgecolor='k', label=name)

    # Overlay single-asset corners
    for name, w in single_asset_ports.items():
        mu = float(w @ means)
        sig = float(np.sqrt(w @ cov @ w))
        plt.scatter([sig], [mu], marker='^', s=80, edgecolor='k', label=name)

    plt.legend(loc='lower right')
    return save_fig(fig, "Portfolio_Scatter_30000.png")

# ==============================
# Check Memo
# ==============================
def check_memo(weights: np.ndarray, name: str, hit_prob: float) -> str:
    expected_arith = float(weights @ means)
    required = (TARGET_A / INIT_TOTAL_A) ** (1/YEARS_A) - 1
    return (f"{name}: Expected arithmetic return ≈ {expected_arith*100:.2f}%/yr, "
            f"required CAGR to reach 80M in 20 years ≈ {required*100:.2f}%/yr; "
            f"simulated target-hit probability was {hit_prob*100:.2f}%.")

# ==============================
# Main
# ==============================
def main():
    # ---- Part A
    results_A = {name: simulate_part_a(w, name) for name, w in portfolios.items()}

    # Part A tables -> CSV
    rows_A = []
    for name, res in results_A.items():
        q = res["qvals"]
        rows_A.append({
            "Portfolio": name,
            "Hit Prob (>=80M) [%]": round(res["hit_prob"]*100, 2),
            "Floor Fail Rate [%]": round(res["floor_fail_rate"]*100, 2),
            "Final P70": round(q["P70"], 1),
            "Final P50": round(q["P50"], 1),
            "Final P15": round(q["P15"], 1),
            "Final P05": round(q["P05"], 1),
        })
        # Plots for A
        plot_partA_hist(name, res)
        plot_partA_paths(name, res)

    df_A = pd.DataFrame(rows_A)
    df_A.to_csv(os.path.join(OUTDIR, "PartA_Summary.csv"), index=False)

    # ---- Part B
    results_B = {name: simulate_part_b(w, name, start_wealth=8000.0, start_cost_basis=8000.0, tax_rate=0.2542)
                 for name, w in portfolios.items()}

    rows_B = []
    for name, res in results_B.items():
        q = res["qvals"]
        rows_B.append({
            "Portfolio": name,
            "Success Rate [%]": round(res["success_rate"]*100, 2),
            "Final < 8000 [%]": round(res["p_under_principal"]*100, 2),
            "Final >= 8000 [%]": round(res["p_over_principal"]*100, 2),
            "Final P70": round(q["P70"], 1),
            "Final P50": round(q["P50"], 1),
            "Final P35": round(q["P35"], 1),
            "Final P20": round(q["P20"], 1),
            "Final P15": round(q["P15"], 1),
            "Final P10": round(q["P10"], 1),
        })
        plot_partB_paths(name, res)

    df_B = pd.DataFrame(rows_B)
    df_B.to_csv(os.path.join(OUTDIR, "PartB_Summary.csv"), index=False)

    # Additional Part B overlays
    plot_B_terminal_overlay(results_B)

    # ---- Linked (A->B)
    results_linked = {name: simulate_linked(name, w, results_A[name]) for name, w in portfolios.items()}

    rows_L = []
    for name, link in results_linked.items():
        rows_L.append({
            "Portfolio": name,
            "A Hit Prob (>=80M) [%]": round(results_A[name]["hit_prob"]*100, 2),
            "B Success | Hit [%]": round((0 if np.isnan(link["B_success_in_hit"]) else link["B_success_in_hit"])*100, 2),
            "Linked Success [%]": round((0 if np.isnan(link["linked_success"]) else link["linked_success"])*100, 2),
            "Hit Trials (count)": int(link["hit_count"]),
        })
        plot_linked_paths(name, link)

    df_L = pd.DataFrame(rows_L)
    df_L.to_csv(os.path.join(OUTDIR, "Linked_AtoB_Summary.csv"), index=False)

    # ---- Success rate comparison bar
    plot_success_rate_bars(results_A, results_B)

    # ---- Risk/Return map (30k ports)
    plot_risk_return_map(portfolios)

    # ---- Check memo
    memos = [check_memo(w, n, results_A[n]['hit_prob']) for n, w in portfolios.items()]
    with open(os.path.join(OUTDIR, "check_memo.txt"), "w") as f:
        f.write("\n".join(memos))

    # Print quick console summary
    print("=== Part A Summary ===")
    print(df_A.to_string(index=False))
    print("\n=== Part B Summary ===")
    print(df_B.to_string(index=False))
    print("\n=== Linked A→B Summary ===")
    print(df_L.to_string(index=False))
    print("\n=== Check Memo ===")
    print("\n".join(memos))
    print(f"\nAll outputs saved in: {os.path.abspath(OUTDIR)}")

if __name__ == "__main__":
    main()