資産配分完全ガイド：ポートフォリオ戦略からリバランシングまで

# pip install numpy pandas yfinance scipy matplotlib seaborn

import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from scipy.optimize import minimize
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')

# ポートフォリオの資産リスト（米国ETFベース）
TICKERS = {
    'SPY':  '米国大型株（S&P 500）',
    'VEA':  '先進国株（米国除く）',
    'VWO':  '新興国株',
    'BND':  '米国債券',
    'TLT':  '米国長期国債',
    'GLD':  '金',
    'VNQ':  '米国REIT',
    'PDBC': 'コモディティ',
}

TRADING_DAYS = 252
RISK_FREE_RATE = 0.04  # 4%を無リスク収益率として使用

def download_prices(tickers, start='2015-01-01', end='2024-01-01'):
    """調整済み終値をダウンロードする"""
    print(f"データをダウンロード中: {len(tickers)}銘柄...")
    raw = yf.download(list(tickers.keys()), start=start, end=end)
    prices = raw['Adj Close'].dropna(how='all').dropna()
    print(f"取得完了: {len(prices)}日分 ({prices.index[0].date()} 〜 {prices.index[-1].date()})")
    return prices

prices = download_prices(TICKERS)
returns = prices.pct_change().dropna()
print(f"\nリターン行列のサイズ: {returns.shape}")

def analyze_asset_performance(returns, tickers, trading_days=252, rf=0.04):
    """
    個別資産の年率化統計量を計算する
    returns: 日次リターンのDataFrame
    tickers: {ティッカー: 名称} の辞書
    """
    results = []
    for ticker in tickers:
        if ticker not in returns.columns:
            continue
        r = returns[ticker]
        ann_ret = r.mean() * trading_days
        ann_vol = r.std() * np.sqrt(trading_days)
        sharpe = (ann_ret - rf) / ann_vol

        # 最大ドローダウン計算
        cum_ret = (1 + r).cumprod()
        rolling_max = cum_ret.cummax()
        drawdown = (cum_ret - rolling_max) / rolling_max
        max_dd = drawdown.min()

        results.append({
            '資産名': tickers[ticker],
            '年率リターン': ann_ret,
            '年率変動率': ann_vol,
            'シャープ比率': sharpe,
            '最大ドローダウン': max_dd,
        })

    return pd.DataFrame(results).set_index('資産名')

stats = analyze_asset_performance(returns, TICKERS, TRADING_DAYS, RISK_FREE_RATE)

print("\n=== 個別資産パフォーマンス分析（年率） ===")
for col in ['年率リターン', '年率変動率', '最大ドローダウン']:
    stats[col] = stats[col].map('{:.2%}'.format)
stats['シャープ比率'] = stats['シャープ比率'].map('{:.3f}'.format)
print(stats.to_string())

def portfolio_stats(weights, mean_returns, cov_matrix, rf=0.04):
    """
    ポートフォリオの年率化統計量を返す
    (年率リターン, 年率変動率, シャープ比率)
    """
    w = np.array(weights)
    ret = np.dot(w, mean_returns) * TRADING_DAYS
    vol = np.sqrt(np.dot(w, np.dot(cov_matrix * TRADING_DAYS, w)))
    sharpe = (ret - rf) / vol
    return ret, vol, sharpe


def simulate_random_portfolios(returns, n=10000, rf=0.04):
    """
    ランダムポートフォリオのシミュレーションで効率的フロンティアを近似
    """
    mean_ret = returns.mean()
    cov_mat = returns.cov()
    n_assets = len(mean_ret)

    all_ret = np.zeros(n)
    all_vol = np.zeros(n)
    all_sharpe = np.zeros(n)
    all_weights = np.zeros((n, n_assets))

    np.random.seed(42)
    for i in range(n):
        # ディリクレ分布でランダムウェイトを生成（合計=1）
        w = np.random.dirichlet(np.ones(n_assets))
        all_weights[i] = w
        r, v, s = portfolio_stats(w, mean_ret, cov_mat, rf)
        all_ret[i] = r
        all_vol[i] = v
        all_sharpe[i] = s

    return all_ret, all_vol, all_sharpe, all_weights


mean_ret = returns.mean()
cov_mat = returns.cov()

p_ret, p_vol, p_sharpe, p_weights = simulate_random_portfolios(
    returns, n=10000, rf=RISK_FREE_RATE
)

print(f"\nシミュレーション完了: {len(p_ret):,}ポートフォリオ")
print(f"最大シャープ比率: {p_sharpe.max():.3f}")
best_idx = np.argmax(p_sharpe)
print(f"最大シャープ時のリターン: {p_ret[best_idx]:.2%}")
print(f"最大シャープ時の変動率: {p_vol[best_idx]:.2%}")

def maximize_sharpe_ratio(returns, rf=0.04):
    """
    scipy.optimize.minimizeを使ってシャープ比率を最大化する
    (負のシャープ比率を最小化する問題に変換)
    """
    mean_ret = returns.mean()
    cov_mat = returns.cov()
    n = len(mean_ret)
    tickers = returns.columns.tolist()

    def neg_sharpe(w):
        ret, vol, _ = portfolio_stats(w, mean_ret, cov_mat, rf)
        return -((ret - rf) / vol)

    # 制約条件: 重みの合計 = 1
    constraints = [{'type': 'eq', 'fun': lambda w: w.sum() - 1.0}]
    # 境界条件: 各資産の重みは0以上1以下（空売りなし）
    bounds = [(0.0, 1.0)] * n
    # 初期値: 等ウェイト
    w0 = np.ones(n) / n

    result = minimize(
        neg_sharpe, w0,
        method='SLSQP',
        bounds=bounds,
        constraints=constraints,
        options={'maxiter': 2000, 'ftol': 1e-9}
    )

    return result.x, tickers


opt_w, opt_tickers = maximize_sharpe_ratio(returns, RISK_FREE_RATE)
opt_ret, opt_vol, opt_sharpe = portfolio_stats(opt_w, mean_ret, cov_mat, RISK_FREE_RATE)

print("\n=== シャープ比率最大化ポートフォリオ ===")
print(f"期待年率リターン: {opt_ret:.2%}")
print(f"年率変動率:       {opt_vol:.2%}")
print(f"シャープ比率:     {opt_sharpe:.3f}")
print("\n資産配分:")
for ticker, w in zip(opt_tickers, opt_w):
    if w > 0.001:
        label = TICKERS.get(ticker, ticker)
        print(f"  {label:25s}: {w:.1%}")

def minimize_portfolio_variance(returns):
    """最小分散ポートフォリオを求める"""
    mean_ret = returns.mean()
    cov_mat = returns.cov()
    n = len(mean_ret)

    def portfolio_vol(w):
        return np.sqrt(np.dot(w, np.dot(cov_mat * TRADING_DAYS, w)))

    constraints = [{'type': 'eq', 'fun': lambda w: w.sum() - 1.0}]
    bounds = [(0.0, 1.0)] * n
    w0 = np.ones(n) / n

    result = minimize(
        portfolio_vol, w0,
        method='SLSQP',
        bounds=bounds,
        constraints=constraints
    )
    return result.x


mv_w = minimize_portfolio_variance(returns)
mv_ret, mv_vol, mv_sharpe = portfolio_stats(mv_w, mean_ret, cov_mat, RISK_FREE_RATE)

print("\n=== 最小分散ポートフォリオ ===")
print(f"期待年率リターン: {mv_ret:.2%}")
print(f"年率変動率:       {mv_vol:.2%}")
print(f"シャープ比率:     {mv_sharpe:.3f}")
print("\n資産配分:")
for ticker, w in zip(opt_tickers, mv_w):
    if w > 0.001:
        label = TICKERS.get(ticker, ticker)
        print(f"  {label:25s}: {w:.1%}")

def plot_frontier(p_ret, p_vol, p_sharpe, opt_w, mv_w,
                  mean_ret, cov_mat, tickers_dict, rf=0.04):
    """効率的フロンティアと主要ポートフォリオを描画"""
    fig, ax = plt.subplots(figsize=(13, 9))

    # ランダムポートフォリオの散布図（シャープ比率で色付け）
    sc = ax.scatter(p_vol, p_ret, c=p_sharpe, cmap='plasma',
                    alpha=0.4, s=8, zorder=2)
    cbar = plt.colorbar(sc, ax=ax)
    cbar.set_label('シャープ比率', fontsize=12)

    # 最大シャープ比率ポートフォリオ
    r_o, v_o, _ = portfolio_stats(opt_w, mean_ret, cov_mat, rf)
    ax.scatter(v_o, r_o, marker='*', color='red', s=700, zorder=6,
               label=f'最大シャープ比率 ({r_o:.1%} / {v_o:.1%})')

    # 最小分散ポートフォリオ
    r_m, v_m, _ = portfolio_stats(mv_w, mean_ret, cov_mat, rf)
    ax.scatter(v_m, r_m, marker='D', color='gold', s=250, zorder=6,
               label=f'最小分散 ({r_m:.1%} / {v_m:.1%})')

    # 個別資産のプロット
    for ticker in returns.columns:
        a_ret = mean_ret[ticker] * TRADING_DAYS
        a_vol = returns[ticker].std() * np.sqrt(TRADING_DAYS)
        ax.scatter(a_vol, a_ret, marker='s', color='lightblue',
                   edgecolors='navy', s=100, zorder=5)
        ax.annotate(
            tickers_dict.get(ticker, ticker),
            (a_vol, a_ret),
            textcoords='offset points',
            xytext=(7, 5),
            fontsize=8
        )

    ax.xaxis.set_major_formatter(mticker.PercentFormatter(xmax=1, decimals=0))
    ax.yaxis.set_major_formatter(mticker.PercentFormatter(xmax=1, decimals=0))
    ax.set_xlabel('年率変動率（リスク）', fontsize=13)
    ax.set_ylabel('年率期待リターン', fontsize=13)
    ax.set_title('効率的フロンティア（モンテカルロ・シミュレーション）', fontsize=15)
    ax.legend(fontsize=11)
    plt.tight_layout()
    plt.savefig('efficient_frontier_ja.png', dpi=150)
    plt.show()
    print("効率的フロンティアを保存しました。")


plot_frontier(p_ret, p_vol, p_sharpe, opt_w, mv_w,
              mean_ret, cov_mat, TICKERS, RISK_FREE_RATE)

def simulate_retirement(
    initial_amount: float,
    monthly_contrib: float,
    weights: np.ndarray,
    mean_ret: pd.Series,
    cov_mat: pd.DataFrame,
    years: int = 30,
    n_sims: int = 1000,
    annual_withdrawal: float = 0.0,
) -> dict:
    """
    退職資産のモンテカルロ・シミュレーション
    積立フェーズ（monthly_contrib > 0）と
    取崩フェーズ（annual_withdrawal > 0）の両方に対応
    """
    n_days = years * TRADING_DAYS
    port_daily_mean = np.dot(weights, mean_ret)
    port_daily_std = np.sqrt(np.dot(weights, np.dot(cov_mat, weights)))
    daily_withdrawal = annual_withdrawal / TRADING_DAYS

    paths = np.zeros((n_days + 1, n_sims))
    paths[0] = initial_amount

    np.random.seed(0)
    for sim in range(n_sims):
        for day in range(1, n_days + 1):
            daily_ret = np.random.normal(port_daily_mean, port_daily_std)
            monthly_add = monthly_contrib if (day % 21 == 0) else 0
            paths[day, sim] = (
                paths[day - 1, sim] * (1 + daily_ret)
                + monthly_add
                - daily_withdrawal
            )
            if paths[day, sim] < 0:
                paths[day, sim] = 0  # 資産枯渇

    final = paths[-1]
    ruin_rate = (final == 0).mean()

    return {
        'paths': paths,
        'final': final,
        'median': np.median(final),
        'p25': np.percentile(final, 25),
        'p75': np.percentile(final, 75),
        'p5': np.percentile(final, 5),
        'p95': np.percentile(final, 95),
        'ruin_rate': ruin_rate,
    }


# 積立シミュレーション：初期100万円 + 毎月3万円積立 × 30年
result = simulate_retirement(
    initial_amount=1_000_000,
    monthly_contrib=30_000,
    weights=opt_w,
    mean_ret=mean_ret,
    cov_mat=cov_mat,
    years=30,
    n_sims=1000,
)

print("\n=== 30年間の資産積立シミュレーション結果 ===")
print(f"初期投資: 100万円 + 毎月3万円積立")
print(f"中央値 (50th):      {result['median']:>15,.0f} 円")
print(f"楽観 (75th):        {result['p75']:>15,.0f} 円")
print(f"悲観 (25th):        {result['p25']:>15,.0f} 円")
print(f"最悪ケース (5th):   {result['p5']:>15,.0f} 円")
print(f"最良ケース (95th):  {result['p95']:>15,.0f} 円")

# 結果の可視化
fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# シミュレーション軌跡
ax1 = axes[0]
years_axis = np.linspace(0, 30, result['paths'].shape[0])
for i in range(min(150, result['paths'].shape[1])):
    ax1.plot(years_axis, result['paths'][:, i] / 1e6,
             alpha=0.06, color='steelblue', lw=0.8)
ax1.plot(years_axis, np.median(result['paths'], axis=1) / 1e6,
         color='red', lw=2.5, label='中央値')
ax1.plot(years_axis, np.percentile(result['paths'], 5, axis=1) / 1e6,
         color='orange', lw=1.5, linestyle='--', label='5パーセンタイル')
ax1.set_xlabel('年数', fontsize=12)
ax1.set_ylabel('ポートフォリオ価値（百万円）', fontsize=12)
ax1.set_title('30年間のモンテカルロ・シミュレーション', fontsize=13)
ax1.legend()

# 最終資産価値の分布
ax2 = axes[1]
ax2.hist(result['final'] / 1e6, bins=60,
         color='steelblue', edgecolor='white', alpha=0.8)
ax2.axvline(result['median'] / 1e6, color='red',
            linestyle='--', lw=2, label='中央値')
ax2.axvline(result['p5'] / 1e6, color='orange',
            linestyle='--', lw=2, label='5パーセンタイル')
ax2.set_xlabel('30年後の資産価値（百万円）', fontsize=12)
ax2.set_ylabel('頻度', fontsize=12)
ax2.set_title('最終資産価値の分布', fontsize=13)
ax2.legend()

plt.tight_layout()
plt.savefig('monte_carlo_ja.png', dpi=150)
plt.show()

fig, ax = plt.subplots(figsize=(10, 8))
corr = returns.corr()
labels = [TICKERS.get(t, t) for t in corr.columns]
corr_df = corr.copy()
corr_df.columns = labels
corr_df.index = labels

sns.heatmap(
    corr_df,
    annot=True,
    fmt='.2f',
    cmap='RdYlGn',
    vmin=-1,
    vmax=1,
    center=0,
    ax=ax,
    linewidths=0.5,
    annot_kws={'size': 9},
    square=True
)
ax.set_title('資産クラス間の相関係数ヒートマップ', fontsize=14)
plt.tight_layout()
plt.savefig('correlation_heatmap_ja.png', dpi=150)
plt.show()
print("相関ヒートマップを保存しました。")