Factor-risk-model-with-principal-component-analysis/
├── main.py                      # メインプログラム
├── data_feed.py                 # データ取得・処理
├── portfolio_decomposition.py   # PCA実装とリスク分解
├── mappings.py                  # ファクターマッピング
└── pcaMappingsResults.csv       # 分析結果

graph TD
    A[Yahoo Finance] --> B[データ取得]
    B --> C[前処理]
    C --> D[PCA適用]
    D --> E[線形回帰]
    E --> F[ファクターマッピング]
    F --> G[リスク評価]
    G --> H[結果出力]

# リスクファクターの体系的分類
risk_factors = {
    "市場リスク": {
        "米国株式": ["^NDX"],
        "欧州株式": ["^STOXX50E"],
        "アジア株式": ["^HSI"]
    },
    "金利リスク": {
        "短期金利": ["SHY"],
        "中期金利": ["IEF"]
    },
    "信用リスク": {
        "投資適格": ["LQD"],
        "ハイイールド": ["HYG"]
    },
    "その他のリスク": {
        "ボラティリティ": ["VIXY"],
        "インフレーション": ["TIP"],
        "国際分散": ["IBND"]
    }
}

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import numpy as np
import pandas as pd

class PCAFactorModel:
    def __init__(self, n_components=None):
        """
        PCAファクターモデルの初期化

        Args:
            n_components: 保持する主成分の数
        """
        self.n_components = n_components
        self.pca = PCA(n_components=n_components)
        self.scaler = StandardScaler()

    def fit(self, returns_data):
        """
        リターンデータにPCAを適用

        Args:
            returns_data: 証券リターンのDataFrame
        """
        # データの標準化
        scaled_returns = self.scaler.fit_transform(returns_data)

        # PCAの適用
        self.pca_factors = self.pca.fit_transform(scaled_returns)

        # 説明分散比率
        self.explained_variance_ratio = self.pca.explained_variance_ratio_

        return self

    def get_factor_loadings(self):
        """ファクターローディングの取得"""
        return pd.DataFrame(
            self.pca.components_.T,
            columns=[f'PC{i+1}' for i in range(self.pca.n_components_)],
            index=self.returns_data.columns
        )

    def calculate_systematic_risk_ratio(self, security_returns, n_factors=3):
        """
        システマティックリスク比率の計算

        Args:
            security_returns: 個別証券のリターン
            n_factors: 使用する主成分の数
        """
        # 主要ファクターによる回帰
        factors = self.pca_factors[:, :n_factors]

        # 線形回帰
        from sklearn.linear_model import LinearRegression
        model = LinearRegression()
        model.fit(factors, security_returns)

        # R²値がシステマティックリスク比率
        systematic_risk_ratio = model.score(factors, security_returns)

        return systematic_risk_ratio

def map_securities_to_factors(security_returns, factor_returns):
    """
    証券をリスクファクターにマッピング

    Args:
        security_returns: 個別証券のリターン時系列
        factor_returns: リスクファクターのリターン時系列

    Returns:
        mapping_results: マッピング係数と統計情報
    """
    results = {}

    for security in security_returns.columns:
        # 線形回帰モデル
        model = LinearRegression()
        X = factor_returns.values
        y = security_returns[security].values

        # 欠損値の処理
        mask = ~(np.isnan(X).any(axis=1) | np.isnan(y))
        X_clean = X[mask]
        y_clean = y[mask]

        if len(X_clean) > 0:
            model.fit(X_clean, y_clean)

            # マッピング係数（ファクターエクスポージャー）
            coefficients = model.coef_

            # 決定係数（R²）
            r_squared = model.score(X_clean, y_clean)

            # 予測値
            y_pred = model.predict(X_clean)

            # 相関係数
            correlation = np.corrcoef(y_clean, y_pred)[0, 1]

            results[security] = {
                'coefficients': dict(zip(factor_returns.columns, coefficients)),
                'r_squared': r_squared,
                'correlation': correlation,
                'systematic_risk_ratio': r_squared  # R²はシステマティックリスクの比率
            }

    return results

import yfinance as yf
from datetime import datetime

def fetch_market_data(symbols, start_date, end_date):
    """
    Yahoo Financeから市場データを取得

    Args:
        symbols: ティッカーシンボルのリスト
        start_date: 開始日
        end_date: 終了日

    Returns:
        price_data: 調整済み終値のDataFrame
    """
    data = {}

    for symbol in symbols:
        try:
            ticker = yf.Ticker(symbol)
            hist = ticker.history(start=start_date, end=end_date)
            data[symbol] = hist['Close']
        except Exception as e:
            print(f"Error fetching {symbol}: {e}")

    # DataFrameに結合
    price_data = pd.DataFrame(data)

    return price_data

def clean_price_data(price_data, missing_threshold=0.1):
    """
    価格データのクリーニング

    Args:
        price_data: 価格データのDataFrame
        missing_threshold: 欠損値の許容割合

    Returns:
        cleaned_data: クリーニング済みデータ
    """
    # 欠損値の割合を計算
    missing_ratio = price_data.isnull().sum() / len(price_data)

    # 閾値を超える列を削除
    valid_columns = missing_ratio[missing_ratio < missing_threshold].index
    cleaned_data = price_data[valid_columns].copy()

    # 線形補間
    cleaned_data = cleaned_data.interpolate(method='linear', limit_direction='both')

    # 残りのNaNを前方埋め
    cleaned_data = cleaned_data.fillna(method='ffill')

    # リターンの計算
    returns = cleaned_data.pct_change().dropna()

    return returns

def calculate_var(returns, confidence_level=0.95, method='historical'):
    """
    Value at Riskの計算

    Args:
        returns: リターンの時系列
        confidence_level: 信頼水準
        method: 計算方法（'historical' or 'parametric'）

    Returns:
        var: VaR値
    """
    if method == 'historical':
        # ヒストリカルVaR
        var = np.percentile(returns, (1 - confidence_level) * 100)

    elif method == 'parametric':
        # パラメトリックVaR（正規分布仮定）
        mean = returns.mean()
        std = returns.std()
        from scipy.stats import norm
        var = norm.ppf(1 - confidence_level, mean, std)

    return var

def calculate_ewma_var(returns, lambda_param=0.94, confidence_level=0.95):
    """
    指数加重移動平均（EWMA）を使用したVaRの計算

    Args:
        returns: リターンの時系列
        lambda_param: 減衰パラメータ
        confidence_level: 信頼水準
    """
    # EWMA分散の計算
    ewma_variance = returns.ewm(alpha=1-lambda_param).var()
    ewma_std = np.sqrt(ewma_variance.iloc[-1])

    # VaRの計算
    from scipy.stats import norm
    var = norm.ppf(1 - confidence_level) * ewma_std

    return var

import matplotlib.pyplot as plt
import seaborn as sns

def visualize_risk_decomposition(mapping_results, top_n=10):
    """
    リスク分解の可視化

    Args:
        mapping_results: ファクターマッピングの結果
        top_n: 表示する上位証券数
    """
    # システマティックリスク比率でソート
    sorted_securities = sorted(
        mapping_results.items(),
        key=lambda x: x[1]['systematic_risk_ratio'],
        reverse=True
    )[:top_n]

    # プロット作成
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

    # システマティックリスク比率
    securities = [s[0] for s in sorted_securities]
    systematic_risks = [s[1]['systematic_risk_ratio'] for s in sorted_securities]
    idiosyncratic_risks = [1 - s[1]['systematic_risk_ratio'] for s in sorted_securities]

    # 積み上げ棒グラフ
    ax1.bar(securities, systematic_risks, label='Systematic Risk')
    ax1.bar(securities, idiosyncratic_risks, bottom=systematic_risks, 
            label='Idiosyncratic Risk', alpha=0.7)
    ax1.set_xlabel('Securities')
    ax1.set_ylabel('Risk Proportion')
    ax1.set_title('Risk Decomposition by Security')
    ax1.legend()
    ax1.set_xticklabels(securities, rotation=45)

    # ファクターエクスポージャーのヒートマップ
    exposure_data = []
    for security, results in sorted_securities:
        exposure_data.append(list(results['coefficients'].values()))

    exposure_df = pd.DataFrame(
        exposure_data,
        index=securities,
        columns=list(sorted_securities[0][1]['coefficients'].keys())
    )

    sns.heatmap(exposure_df, cmap='RdBu_r', center=0, ax=ax2,
                cbar_kws={'label': 'Factor Exposure'})
    ax2.set_title('Factor Exposure Heatmap')
    ax2.set_xlabel('Risk Factors')

    plt.tight_layout()
    plt.show()

def run_factor_risk_analysis():
    """ファクターリスク分析の実行"""

    # 1. パラメータ設定
    start_date = datetime(2015, 1, 1)
    end_date = datetime(2020, 1, 1)

    # リスクファクター
    factor_symbols = ['^NDX', '^STOXX50E', '^HSI', 'IEF', 'SHY', 
                      'VIXY', 'LQD', 'HYG', 'IBND', 'TIP']

    # 分析対象証券
    security_symbols = ['AAPL', 'MSFT', 'GOOGL', 'AMZN', 'FB', 
                        'BRK-B', 'JPM', 'JNJ', 'V', 'PG']

    # 2. データ取得
    print("Fetching factor data...")
    factor_data = fetch_market_data(factor_symbols, start_date, end_date)

    print("Fetching security data...")
    security_data = fetch_market_data(security_symbols, start_date, end_date)

    # 3. データクリーニング
    factor_returns = clean_price_data(factor_data)
    security_returns = clean_price_data(security_data)

    # 4. PCA適用
    print("Applying PCA...")
    pca_model = PCAFactorModel(n_components=5)
    pca_model.fit(factor_returns)

    print(f"Explained variance ratio: {pca_model.explained_variance_ratio}")

    # 5. ファクターマッピング
    print("Mapping securities to factors...")
    mapping_results = map_securities_to_factors(security_returns, factor_returns)

    # 6. 結果の出力
    results_df = pd.DataFrame(mapping_results).T
    results_df.to_csv('factor_mapping_results.csv')

    # 7. 可視化
    visualize_risk_decomposition(mapping_results)

    return mapping_results

# 実行
if __name__ == "__main__":
    results = run_factor_risk_analysis()

def interpret_results(mapping_results):
    """分析結果の解釈"""

    for security, results in mapping_results.items():
        print(f"\n=== {security} ===")
        print(f"Systematic Risk Ratio: {results['systematic_risk_ratio']:.2%}")
        print(f"Correlation with factors: {results['correlation']:.3f}")

        # 主要なファクターエクスポージャー
        exposures = sorted(
            results['coefficients'].items(),
            key=lambda x: abs(x[1]),
            reverse=True
        )

        print("\nTop Factor Exposures:")
        for factor, exposure in exposures[:3]:
            print(f"  {factor}: {exposure:.3f}")

class PortfolioRiskManager:
    def __init__(self, mapping_results):
        self.mapping_results = mapping_results

    def calculate_portfolio_risk(self, weights, factor_covariance):
        """ポートフォリオリスクの計算"""
        # ファクターエクスポージャーの集計
        portfolio_exposure = np.zeros(len(factor_covariance))

        for security, weight in weights.items():
            if security in self.mapping_results:
                exposures = self.mapping_results[security]['coefficients']
                for i, (factor, exposure) in enumerate(exposures.items()):
                    portfolio_exposure[i] += weight * exposure

        # ポートフォリオ分散
        portfolio_variance = portfolio_exposure @ factor_covariance @ portfolio_exposure.T
        portfolio_risk = np.sqrt(portfolio_variance)

        return portfolio_risk

    def optimize_risk_parity(self, factor_covariance):
        """リスクパリティ最適化"""
        from scipy.optimize import minimize

        n_assets = len(self.mapping_results)

        # 目的関数：各資産のリスク寄与を均等化
        def risk_parity_objective(weights):
            portfolio_risk = self.calculate_portfolio_risk(
                dict(zip(self.mapping_results.keys(), weights)),
                factor_covariance
            )

            # 各資産のリスク寄与
            marginal_contributions = []
            for i in range(n_assets):
                w_temp = weights.copy()
                w_temp[i] += 0.0001
                risk_up = self.calculate_portfolio_risk(
                    dict(zip(self.mapping_results.keys(), w_temp)),
                    factor_covariance
                )
                marginal_contribution = (risk_up - portfolio_risk) / 0.0001
                marginal_contributions.append(weights[i] * marginal_contribution)

            # リスク寄与の分散を最小化
            return np.var(marginal_contributions)

        # 制約条件
        constraints = [
            {'type': 'eq', 'fun': lambda x: np.sum(x) - 1},  # 合計1
            {'type': 'ineq', 'fun': lambda x: x}  # 非負制約
        ]

        # 初期値
        x0 = np.ones(n_assets) / n_assets

        # 最適化
        result = minimize(
            risk_parity_objective,
            x0,
            method='SLSQP',
            constraints=constraints
        )

        return dict(zip(self.mapping_results.keys(), result.x))

def stress_test_portfolio(portfolio_weights, mapping_results, stress_scenarios):
    """
    ポートフォリオのストレステスト

    Args:
        portfolio_weights: ポートフォリオの構成比
        mapping_results: ファクターマッピング結果
        stress_scenarios: ストレスシナリオ

    Returns:
        stress_results: 各シナリオでの損失
    """
    stress_results = {}

    for scenario_name, factor_shocks in stress_scenarios.items():
        portfolio_loss = 0

        for security, weight in portfolio_weights.items():
            if security in mapping_results:
                security_loss = 0
                exposures = mapping_results[security]['coefficients']

                for factor, shock in factor_shocks.items():
                    if factor in exposures:
                        security_loss += exposures[factor] * shock

                portfolio_loss += weight * security_loss

        stress_results[scenario_name] = portfolio_loss

    return stress_results

# ストレスシナリオの例
stress_scenarios = {
    "Market Crash": {
        "^NDX": -0.20,
        "^STOXX50E": -0.25,
        "^HSI": -0.30,
        "VIXY": 0.50
    },
    "Interest Rate Shock": {
        "IEF": -0.10,
        "SHY": -0.05,
        "LQD": -0.15,
        "HYG": -0.20
    },
    "Credit Crisis": {
        "LQD": -0.10,
        "HYG": -0.30,
        "^NDX": -0.15
    }
}