ある混雑度ファクターの紹介(Comomentum)

import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from tqdm import tqdm # 進捗バー表示用

class ComomentumAnalyzer:
    """
    Lou & Polk (2013) に基づく Comomentum (混雑度ファクター) を計算するクラス
    """
    def __init__(self, prices_df, market_series, window_mom=252, window_corr=252):
        """
        :param prices_df: 株価のデータフレーム (index: 日付, columns: 銘柄)
        :param market_series: ベンチマーク(S&P500等)の価格シリーズ
        :param window_mom: モメンタム計測期間 (例: 252営業日 = 1年)
        :param window_corr: 相関計測期間 (例: 252営業日)
        """
        self.prices = prices_df
        self.market = market_series
        self.returns = prices_df.pct_change().dropna(how='all')
        self.market_returns = market_series.pct_change().dropna()
        self.window_mom = window_mom
        self.window_corr = window_corr

    def _get_momentum_winners(self, date, top_n_percent=0.1):
        """
        指定された日付時点でのモメンタム上位銘柄（Winner）を特定する
        ※ 12ヶ月リターン (直近1ヶ月除外) を使用するのが一般的
        """
        # 指定日から過去の期間をスライス
        past_prices = self.prices.loc[:date].iloc[-self.window_mom:]

        if len(past_prices) < self.window_mom:
            return []

        # リターン計算（t-252 ~ t-21 : 直近1ヶ月をスキップする慣習）
        lag = 21
        mom_score = (past_prices.iloc[-lag] / past_prices.iloc[0]) - 1

        # 上位N%を抽出
        n_stocks = int(len(mom_score.dropna()) * top_n_percent)
        if n_stocks < 2:
            return []

        winners = mom_score.nlargest(n_stocks).index.tolist()
        return winners

    def _calculate_residuals(self, stock_returns, market_returns):
        """
        市場要因を除去した「残差」を計算する (CAPMモデル: r_i = alpha + beta * r_m + epsilon)
        """
        residuals = pd.DataFrame(index=stock_returns.index, columns=stock_returns.columns)

        # 市場リターンの定数項追加 (alpha用)
        X = sm.add_constant(market_returns)

        for stock in stock_returns.columns:
            y = stock_returns[stock]
            # 欠損除去して回帰
            valid_idx = y.dropna().index.intersection(X.index)
            if len(valid_idx) < len(market_returns) * 0.8: # データ不足はスキップ
                continue

            model = sm.OLS(y.loc[valid_idx], X.loc[valid_idx]).fit()
            residuals.loc[valid_idx, stock] = model.resid

        return residuals.astype(float)

    def compute_history(self, rebalance_freq='ME'):
        """
        時系列でComomentumを計算するメイン関数
        :param rebalance_freq: 計算頻度 (例: 'ME' = 月末ごと)
        """
        # dates = self.returns.resample(rebalance_freq).last().index
        dates = self.returns.index
        comomentum_history = []
        valid_dates = []

        print("Calculating Comomentum History...")
        for date in tqdm(dates):
            if date not in self.returns.index:
                continue

            # 1. その時点での「勝ち組（Winner）」を特定
            winners = self._get_momentum_winners(date)
            if not winners:
                continue

            # 2. Winner銘柄と市場の直近リターンを取得
            window_start = date - pd.Timedelta(days=365) # 近似的に1年前

            subset_returns = self.returns.loc[window_start:date, winners]
            subset_market = self.market_returns.loc[window_start:date]

            # データ長合わせ
            common_idx = subset_returns.index.intersection(subset_market.index)
            subset_returns = subset_returns.loc[common_idx]
            subset_market = subset_market.loc[common_idx]

            if len(subset_market) < self.window_corr * 0.8:
                continue

            # 3. 残差（Residuals）の計算 (直交化プロセス)
            residuals = self._calculate_residuals(subset_returns, subset_market)

            # 4. 残差間の相関行列を計算
            corr_matrix = residuals.corr()

            # 5. 対角成分を除いた平均相関を計算 (Comomentumスコア)
            # 上三角行列の要素のみを取得
            mask = np.triu(np.ones_like(corr_matrix, dtype=bool), k=1)
            avg_corr = corr_matrix.where(mask).stack().mean()

            comomentum_history.append(avg_corr)
            valid_dates.append(date)

        return pd.Series(comomentum_history, index=valid_dates, name='Comomentum')

# --- 以下、使用例（ダミーデータ生成） ---

def generate_dummy_data(n_stocks=50, n_days=1000):
    """検証用のランダムデータ生成"""
    np.random.seed(42)
    dates = pd.date_range(start='2020-01-01', periods=n_days, freq='B')

    # 市場データ
    market = pd.Series(np.cumprod(1 + np.random.normal(0.0005, 0.01, n_days)), index=dates)

    # 個別株データ（市場連動 + 固有ノイズ）
    prices_data = {}
    for i in range(n_stocks):
        beta = np.random.uniform(0.5, 1.5)
        noise = np.random.normal(0, 0.02, n_days)
        # わざとらしい「混雑」を作るため、一部期間に共通ノイズを混ぜる
        if i < 10: # 一部の銘柄群
            noise[500:600] += np.random.normal(0, 0.05, 100) # 暴れる期間

        ret = (market.pct_change().fillna(0) * beta) + noise
        prices_data[f'Stock_{i}'] = np.cumprod(1 + ret)

    return pd.DataFrame(prices_data, index=dates), market

# --- 実行セクション ---
if __name__ == "__main__":
    # 1. データ準備
    prices, market = generate_dummy_data() # 本番は株式とTOPIXなどのインデックスデータを入れる。

    # 2. インスタンス化と計算
    analyzer = ComomentumAnalyzer(prices, market)
    comomentum_series = analyzer.compute_history(rebalance_freq='BME') # 月次計算

    # 3. プロット
    plt.figure(figsize=(10, 6))

    # --- 修正ポイント: データを日次に揃えて（Reindex）、前方の値を埋める（ffill） ---
    # これにより、Comomentumも2020年1月からのデータ枠を持つ（値はNaNだが枠はある）
    daily_comomentum = comomentum_series.reindex(market.index).ffill()

    # 上段: Comomentum
    ax1 = plt.subplot(2, 1, 1)
    daily_comomentum.plot(ax=ax1, color='darkred', lw=2)
    ax1.set_title('Comomentum Factor (Crowding Score)')
    ax1.set_ylabel('Avg Pairwise Correlation of Residuals')
    ax1.grid(True, alpha=0.3)

    # 下段: 市場価格
    ax2 = plt.subplot(2, 1, 2, sharex=ax1)
    market.plot(ax=ax2, color='navy')
    ax2.set_title('Market Index')
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()