import pandas as pd
import yfinance as yf

#code = "7337.T"  # ひろぎんホールディングス
code = "7186.T"  # 横浜フィナンシャルグループ

# Tickerオブジェクトの作成
ticker = yf.Ticker(code)

# 過去6年間の日次データを取得（period='6y', interval='1d'）
# 補正前（調整前）の株価を取得するため、auto_adjust=False を指定
df = ticker.history(period="6y", interval="1d", auto_adjust=False)

# index（DatetimeIndex）をカラムとしてリセットし、元のコードと同じ列構成にする
df = df.reset_index()

# 必要な列のみを抽出して並び替え
# yfinanceはデフォルトで 'Date', 'Open', 'High', 'Low', 'Close', 'Volume' などを返します
df = df[['Date', 'Open', 'High', 'Low', 'Close', 'Volume']]

# タイムゾーンを日本標準時（Asia/Tokyo）に変換
df['Date'] = df['Date'].dt.tz_convert('Asia/Tokyo')

# Date列のフォーマットを日付のみに変更する（時分秒をカットする）
df['Date'] = df['Date'].dt.strftime('%Y-%m-%d')

# csv出力。小数点以下は２桁までとする。インデクス（通番）は出力しないようにする
df.to_csv(f"{code}.csv", float_format='%.2f', index=False)

Date,Open,High,Low,Close,Volume
2020-10-02,690.00,700.00,650.00,650.00,1645800
2020-10-05,660.00,663.00,627.00,631.00,773700
2020-10-06,621.00,624.00,603.00,605.00,1001200
2020-10-07,602.00,608.00,584.00,599.00,827200
2020-10-08,603.00,614.00,601.00,609.00,599100
2020-10-09,612.00,617.00,605.00,613.00,328600
...

stockA = pd.read_csv("7186.T.csv", parse_dates=["Date"], index_col="Date")
stockB = pd.read_csv("7387.T.csv", parse_dates=["Date"], index_col="Date")
combined_data = pd.concat([stockA["Close"], stockB["Close"]], axis=1).dropna()
combined_data.columns = ["stockA", "stockB"]
combined_data = combined_data.asfreq("B").ffill()

_, p_value, _ = coint(window_data["NYK"], window_data["MOL"])
eg_exists = p_value < 0.05

jo_result = coint_johansen(window_data, det_order=0, k_ar_diff=1)
trace_stat = jo_result.lr1[0]    # トレース統計量（r=0 の検定）
crit_val   = jo_result.cvt[0, 1] # 5% 臨界値（= 15.49）
jo_exists  = trace_stat > crit_val

window_results.append(
    (window_data.index[0], window_data.index[-1],
     eg_exists, jo_exists, p_value, trace_stat)
)

roll_cov = log_data["stockA"].rolling(window_len).cov(log_data["stockB"])
roll_var = log_data["stockB"].rolling(window_len).var()
beta_wf  = roll_cov / roll_var

alpha_wf = (
    log_data["stockA"].rolling(window_len).mean()
    - beta_wf * log_data["stockB"].rolling(window_len).mean()
)

spread_wf = log_data["stockA"] - beta_wf * log_data["stockB"] - alpha_wf

valid_mask = spread_wf.notna()
scaler_spread = MinMaxScaler(feature_range=(-1, 1))
_norm = scaler_spread.fit_transform(spread_wf[valid_mask].values.reshape(-1, 1)).flatten()
stationary_series = pd.Series(np.nan, index=combined_data.index)
stationary_series[valid_mask] = _norm

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from statsmodels.tsa.stattools import adfuller, coint
from statsmodels.tsa.vector_ar.vecm import coint_johansen
import japanize_matplotlib
from sklearn.preprocessing import MinMaxScaler
from numpy.linalg import LinAlgError

# --- 1. データの読み込み---
stockA = pd.read_csv("7186.T.csv", parse_dates=["Date"], index_col="Date")
stockB = pd.read_csv("7337.T.csv", parse_dates=["Date"], index_col="Date")
nameA = "横浜フィナンシャルグループ"
nameB = "ひろぎんホールディングス"
# 終値のみを使用し、結合と欠損値処理（全期間を学習期間とする）
combined_data = pd.concat([stockA["Close"], stockB["Close"]], axis=1).dropna()
combined_data.columns = ["stockA", "stockB"]
combined_data = combined_data.asfreq("B").ffill()

# 対数株価（統計処理の基本系列）— 可視化のみ生株価を使用
log_data = np.log(combined_data)

# --- 2. 単位根検定（ADF検定）--- ※対数株価で実施
adf_stockA = adfuller(log_data["stockA"])
adf_stockB = adfuller(log_data["stockB"])
print(
    f"{nameA} ADF 検定統計量: {adf_stockA[0]:.4f}, p-値: {adf_stockA[1]:.4f} -> "
    f"{'単位根あり(非定常)' if adf_stockA[1] > 0.05 else '単位根なし(定常)'}"
)
print(
    f"{nameB} ADF 検定統計量: {adf_stockB[0]:.4f}, p-値: {adf_stockB[1]:.4f} -> "
    f"{'単位根あり(非定常)' if adf_stockB[1] > 0.05 else '単位根なし(定常)'}"
)
print("-" * 60)

# --- 3. ローリング（ウィンドウ）共和分検定 ---
window_size = 6  # 6か月
step = 1  # 1か月
window_len = window_size * 20
step_len = step * 20

JOHANSEN_CRIT = 15.49  # Johansen 5% 臨界値

count_coint_eg = 0
count_coint_jo = 0
window_results = []  # (start, end, eg_exists, jo_exists, p_value, trace_stat)

for i in range(0, len(log_data) - window_len + 1, step_len):
    window_data = log_data.iloc[i : i + window_len]  # 対数株価ウィンドウ

    # Engle-Granger検定（対数株価）
    _, p_value, _ = coint(window_data["stockA"], window_data["stockB"])
    eg_exists = p_value < 0.05
    if eg_exists:
        count_coint_eg += 1

    # Johansen検定（対数株価）
    trace_stat = np.nan
    try:
        jo_result = coint_johansen(window_data, det_order=0, k_ar_diff=1)
        trace_stat = jo_result.lr1[0]
        crit_val = jo_result.cvt[0, 1]
        jo_exists = trace_stat > crit_val
        if jo_exists:
            count_coint_jo += 1
    except LinAlgError as e:
        jo_exists = False
        print(f"エラー: {e} - ウィンドウ: {window_data.index[0].strftime('%Y-%m')}")

    window_results.append(
        (
            window_data.index[0],
            window_data.index[-1],
            eg_exists,
            jo_exists,
            p_value,
            trace_stat,
        )
    )

    print(
        f"期間 {window_data.index[0].strftime('%Y-%m')} - {window_data.index[-1].strftime('%Y-%m')}: "
        f"EG={'あり' if eg_exists else 'なし'}(p={p_value:.4f})  "
        f"Johansen={'あり' if jo_exists else 'なし'}"
        f"({'N/A' if np.isnan(trace_stat) else f'stat={trace_stat:.2f}'})"
    )

print("-" * 60)
print(f"総ウィンドウ数: {len(window_results)}")
print(
    f"共和分検出回数 -> E.G.検定: {count_coint_eg}回, Johansen検定: {count_coint_jo}回"
)
print("-" * 60)

# --- 4. WF（ローリング）スプレッド計算（Panel 2 表示用） ---
# 対数株価を使用: spread = log(A) - β×log(B) - α
# β = Cov(log A, log B) / Var(log B) で弾力性（elasticity）として推定する
# 先頭 window_len-1 点は NaN となる
roll_cov = log_data["stockA"].rolling(window_len).cov(log_data["stockB"])
roll_var = log_data["stockB"].rolling(window_len).var()
beta_wf  = roll_cov / roll_var                                     # 時変ヘッジ比率（弾力性）
# 切片 α = mean(log A) - β × mean(log B) を各窓で計算し、残差のみを取り出す
alpha_wf = (
    log_data["stockA"].rolling(window_len).mean()
    - beta_wf * log_data["stockB"].rolling(window_len).mean()
)
spread_wf = log_data["stockA"] - beta_wf * log_data["stockB"] - alpha_wf

# 正規化: NaN を除いた有効データで MinMaxScaler を学習し、全インデックスに適用
valid_mask = spread_wf.notna()
scaler_spread = MinMaxScaler(feature_range=(-1, 1))
_norm = scaler_spread.fit_transform(spread_wf[valid_mask].values.reshape(-1, 1)).flatten()
stationary_series = pd.Series(np.nan, index=combined_data.index)
stationary_series[valid_mask] = _norm

print(f"ローリング WF ヘッジ比率 β（最新窓）: {beta_wf.iloc[-1]:.4f}  切片 α: {alpha_wf.iloc[-1]:.4f}")
print(f"β の推移: {beta_wf.min():.4f} 〜 {beta_wf.max():.4f}（平均 {beta_wf.mean():.4f}）")

# ================================================================
# ローリング窓データの準備（Panel 3 用）
# ================================================================
win_starts = np.array([r[0] for r in window_results])
win_ends = np.array([r[1] for r in window_results])
win_mids = np.array([r[0] + (r[1] - r[0]) / 2 for r in window_results])
eg_passes = np.array([r[2] for r in window_results])
jo_passes = np.array([r[3] for r in window_results])
eg_pvals = np.array([r[4] for r in window_results], dtype=float)
jo_traces = np.array([r[5] for r in window_results], dtype=float)

# --- 4. 可視化プロット ---
fig, axs = plt.subplots(
    3,
    1,
    figsize=(14, 13),
    gridspec_kw={"height_ratios": [3, 1.2, 2]},
)

# ================================================================
# Panel 1: 生価格
# ================================================================
ax1 = axs[0]
ax1.plot(
    combined_data.index,
    combined_data["stockA"],
    label=nameA,
    color="blue",
)
ax1.plot(
    combined_data.index,
    combined_data["stockB"],
    label=nameB,
    color="green",
)

# 各ローリング窓の共和分成立状況を背景に薄く着色
for start, end, eg, jo, *_ in window_results:
    if eg and jo:
        bg = "#c8f7c5"  # 薄緑: 両検定合格
    elif eg or jo:
        bg = "#fef9c3"  # 薄黄: 片方合格
    else:
        bg = "#ffd6d6"  # 薄ピンク: 両検定不合格
    ax1.axvspan(start, end, alpha=0.15, color=bg, lw=0)

ax1.set_title("二銘柄の株価推移（生株価・参考用）", fontsize=11)
ax1.set_ylabel("株価（円）")
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.4)

# ================================================================
# Panel 2: WF スプレッド（ルックアヘッドなし）
# ================================================================
ax2 = axs[1]
ax2r = ax2.twinx()  # 右軸: ローリングヘッジ比率 β

# 左軸: 正規化スプレッド（赤）
ax2.plot(
    combined_data.index,
    stationary_series,
    label=f"WF スプレッド（ローリング {window_size}ヶ月 β）",
    color="red",
    lw=1.2,
    zorder=3,
)
ax2.axhline(0, color="red", linestyle=":", lw=0.8, alpha=0.5)

# 右軸: ローリングβ（灰色破線）— βの時変性を確認するための補助情報
ax2r.plot(
    combined_data.index,
    beta_wf,
    label="ヘッジ比率 β（右軸）",
    color="gray",
    lw=0.9,
    linestyle="--",
    alpha=0.6,
    zorder=2,
)
ax2r.set_ylabel("ヘッジ比率 β", color="gray", fontsize=8)
ax2r.tick_params(axis="y", labelcolor="gray")

ax2.set_title(
    f"二銘柄の対数価格乖離度 ── WF スプレッド（切片補正済み・ルックアヘッドなし）\n"
    f"spread = log(A) − β×log(B) − α　各時点の β・α は過去 {window_size}ヶ月データのみで推定",
    fontsize=10,
)
ax2.set_ylabel("乖離度 ([-1,1] 正規化)")
lines_l, labels_l = ax2.get_legend_handles_labels()
lines_r, labels_r = ax2r.get_legend_handles_labels()
ax2.legend(lines_l + lines_r, labels_l + labels_r, fontsize=8)
ax2.grid(True, alpha=0.4)

# ================================================================
# Panel 3 (NEW): ローリング共和分検定結果
# ================================================================
ax3 = axs[2]
ax3r = ax3.twinx()  # 右軸: Johansen 統計量

# 背景色: 両検定の合否組み合わせ
for i in range(len(win_mids)):
    eg, jo = eg_passes[i], jo_passes[i]
    if eg and jo:
        bg = "#d4edda"
    elif eg or jo:
        bg = "#fff3cd"
    else:
        bg = "#f8d7da"
    ax3.axvspan(win_starts[i], win_ends[i], alpha=0.25, color=bg, lw=0, zorder=1)

# EG p値 (左軸, 青)
ax3.step(
    win_mids, eg_pvals, where="mid", color="steelblue", lw=1.5, label="EG p値", zorder=3
)
ax3.axhline(
    0.05, color="steelblue", linestyle="--", lw=1, alpha=0.8, label="EG 閾値 (p=0.05)"
)
ax3.scatter(
    win_mids[eg_passes],
    eg_pvals[eg_passes],
    color="green",
    s=50,
    zorder=5,
    marker="o",
    label="EG 合格",
)
ax3.scatter(
    win_mids[~eg_passes],
    eg_pvals[~eg_passes],
    color="crimson",
    s=50,
    zorder=5,
    marker="x",
    label="EG 不合格",
)
ax3.set_ylabel("EG p値", color="steelblue", fontsize=9)
ax3.tick_params(axis="y", labelcolor="steelblue")
ax3.set_ylim(-0.05, 1.05)

# Johansen 統計量 (右軸, オレンジ)
valid = ~np.isnan(jo_traces)
ax3r.step(
    win_mids[valid],
    jo_traces[valid],
    where="mid",
    color="darkorange",
    lw=1.5,
    label="J統計量",
    zorder=3,
)
ax3r.axhline(
    JOHANSEN_CRIT,
    color="darkorange",
    linestyle="--",
    lw=1,
    alpha=0.8,
    label=f"Johansen 閾値 ({JOHANSEN_CRIT})",
)
ax3r.scatter(
    win_mids[jo_passes & valid],
    jo_traces[jo_passes & valid],
    color="orange",
    s=50,
    zorder=5,
    marker="o",
    label="Johansen 合格",
)
ax3r.scatter(
    win_mids[~jo_passes & valid],
    jo_traces[~jo_passes & valid],
    color="red",
    s=50,
    zorder=5,
    marker="x",
    label="Johansen 不合格",
)
ax3r.set_ylabel("Johansen 統計量", color="darkorange", fontsize=9)
ax3r.tick_params(axis="y", labelcolor="darkorange")

# 凡例の統合
lines_l, labels_l = ax3.get_legend_handles_labels()
lines_r, labels_r = ax3r.get_legend_handles_labels()
patch_both = mpatches.Patch(color="#d4edda", alpha=0.8, label="両検定 合格")
patch_one = mpatches.Patch(color="#fff3cd", alpha=0.8, label="片方 合格")
patch_none = mpatches.Patch(color="#f8d7da", alpha=0.8, label="両検定 不合格")
ax3.legend(
    handles=lines_l + lines_r + [patch_both, patch_one, patch_none],
    labels=labels_l + labels_r + ["両検定 合格", "片方 合格", "両検定 不合格"],
    loc="upper right",
    fontsize=7.5,
    ncol=2,
)

ax3.set_title(
    f"ローリング共和分検定 結果（窓={window_size}ヶ月 / ステップ={step}ヶ月 / "
    f"総{len(window_results)}窓）  "
    f"EG合格: {count_coint_eg}/{len(window_results)}  "
    f"Johansen合格: {count_coint_jo}/{len(window_results)}",
    fontsize=9,
)
ax3.set_xlabel("日付")
ax3.grid(True, alpha=0.3)

# 表示期間: 全ウィンドウが含まれるようデータ全期間に設定
start_plot_date = combined_data.index[0]
end_plot_date = combined_data.index[-1]
for ax in axs:
    ax.set_xlim(start_plot_date, end_plot_date)
ax2r.set_xlim(start_plot_date, end_plot_date)
ax3r.set_xlim(start_plot_date, end_plot_date)

plt.tight_layout(h_pad=1.5)
plt.show()