"""
Harvey, Liu & Zhu(2016)의 "새 팩터는 t값 3을 넘어야 한다"는 주장을 공개 데이터로 작게 확인한다
  1) 미국 팩터(시장·SMB·HML·모멘텀) 월평균 수익률의 t값을 논문 표본 기간과 그 뒤로 나눠 본다
  2) 수익이 전혀 없는 가짜 팩터 316개를 만들어, t값 2와 3을 넘는 것이 몇 개 나오는지 센다
  3) 검정 횟수에 따라 본페로니 기준이 어떻게 올라가는지 계산한다
  4) 규모 × B/M 25개 포트폴리오의 CAPM 절편 25개에 본페로니·홀름·BHY 보정을 적용한다

사용법:
    pip install pandas numpy matplotlib
    python factor_hurdle.py
    python factor_hurdle.py --synthetic     # 난수 데이터로 동작 확인

데이터: 프렌치 교수 데이터 라이브러리의 3팩터·모멘텀 팩터·25개 포트폴리오(시가총액 가중) 월별 수익률, 공개 저장소의 사본.
  - 3팩터와 25개 포트폴리오 사본은 2023년 12월, 모멘텀 사본은 2024년 2월에 끝난다. 모두 2023년 12월까지만 쓴다.
  - t값은 월평균 수익률 / (월 수익률 표준편차 / √개월 수)다. p값은 개월 수가 많아 정규분포로 근사한다.
결과: output/factor_hurdle.txt 와 차트 PNG
"""
import argparse
import os
from statistics import NormalDist

import numpy as np
import pandas as pd

HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "output")
BASE = "https://raw.githubusercontent.com/a91quaini/"
FF3_URL = BASE + "reproduceTFRP/main/data-raw/F-F_Research_Data_Factors.CSV"
P25_URL = BASE + "reproduceTFRP/main/data-raw/25_Portfolios_5x5.CSV"
MOM_URL = BASE + "intrinsicFRP/main/data-raw/F-F_Momentum_Factor.CSV"
END = "2023-12"
N_FAKE, N_SIM, SEED = 316, 2000, 7
# (팩터, 논문 표본 시작, 끝, 그 뒤 시작): SMB·HML은 Fama & French(1993), 모멘텀은 Jegadeesh & Titman(1993)의 표본
SPLITS = [("Mkt-RF", "1963-07", "1991-12", "1992-01"), ("SMB", "1963-07", "1991-12", "1992-01"),
          ("HML", "1963-07", "1991-12", "1992-01"), ("Mom", "1965-01", "1989-12", "1990-01")]
N01 = NormalDist()
lines = []


def log(s=""):
    print(s)
    lines.append(str(s))


def read_monthly(url):
    df = pd.read_csv(url)
    df.columns = [c.strip() for c in df.columns]
    date = df["Date"].astype(str).str.strip()
    df = df[date.str.len() == 6].copy()                    # 월별 행만 (연간 표가 붙은 사본 대비)
    df.index = pd.to_datetime(df.pop("Date").astype(str).str.strip(), format="%Y%m") + pd.offsets.MonthEnd(0)
    return df.astype(float) / 100


def load(synthetic=False):
    """(팩터 표: Mkt-RF·SMB·HML·RF·Mom, 25개 포트폴리오) 월별, 소수, 2023년 12월까지."""
    if synthetic:
        rng = np.random.default_rng(1)
        idx = pd.date_range("1926-07-31", "2023-12-31", freq="ME")
        f = pd.DataFrame(rng.normal([0.006, 0.002, 0.003, 0.003, 0.006], [0.045, 0.03, 0.03, 0.002, 0.045],
                                    (len(idx), 5)), index=idx, columns=["Mkt-RF", "SMB", "HML", "RF", "Mom"])
        p = pd.DataFrame(f[["Mkt-RF"]].values * rng.uniform(0.8, 1.3, 25) + f[["RF"]].values
                         + rng.normal(0.001, 0.02, (len(idx), 25)), index=idx, columns=[f"P{i}" for i in range(25)])
        return f, p
    f = read_monthly(FF3_URL)
    f["Mom"] = read_monthly(MOM_URL)["Mom"]
    return f.loc[:END], read_monthly(P25_URL).loc[:END]


def tstat(x):
    x = np.asarray(x, dtype=float)
    return x.mean() / (x.std(ddof=1) / np.sqrt(len(x)))


def p_two_sided(t):
    return 2 * (1 - N01.cdf(abs(t)))


def bonferroni_cut(n, alpha=0.05):
    """검정 n번, 양쪽 유의수준 alpha에서 하나라도 잘못 걸릴 확률을 alpha로 묶는 t값 기준."""
    return N01.inv_cdf(1 - alpha / (2 * n))


def holm(p, alpha=0.05):
    """홀름: 작은 p값부터 alpha/(m-k+1)과 비교해 처음 실패하는 곳에서 멈춘다. 통과 여부 배열."""
    m, order = len(p), np.argsort(p)
    ok = np.zeros(m, bool)
    for k, i in enumerate(order):
        if p[i] > alpha / (m - k):
            break
        ok[i] = True
    return ok


def bhy(p, alpha=0.05):
    """BHY(Benjamini-Hochberg-Yekutieli): 거짓 발견 비율을 alpha로 묶는다. 검정 사이 상관이 있어도 성립하는 보수적 판."""
    m, order = len(p), np.argsort(p)
    c = np.sum(1 / np.arange(1, m + 1))
    passed = [k for k in range(m) if p[order[k]] <= (k + 1) * alpha / (m * c)]
    ok = np.zeros(m, bool)
    if passed:
        ok[order[:max(passed) + 1]] = True
    return ok


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--synthetic", action="store_true")
    a = ap.parse_args()
    os.makedirs(OUT, exist_ok=True)
    f, p = load(a.synthetic)
    log(f"데이터: 팩터 {f.index[0]:%Y-%m} ~ {f.index[-1]:%Y-%m}, 25개 포트폴리오 {p.index[0]:%Y-%m} ~ {p.index[-1]:%Y-%m}")

    log("\n[1] 팩터 월평균 수익률의 t값: 논문 표본, 그 뒤, 1963-07 이후 전체")
    rows = {}
    for name, s, e, after in SPLITS:
        x = f[name].dropna()
        ins, post, full = x.loc[s:e], x.loc[after:], x.loc["1963-07":]
        rows[name] = {"표본": f"{s}~{e}", "표본 월평균(%)": ins.mean() * 100, "표본 t": tstat(ins),
                      "그 뒤": f"{after}~{END}", "그 뒤 월평균(%)": post.mean() * 100, "그 뒤 t": tstat(post),
                      "전체 t": tstat(full)}
    t1 = pd.DataFrame(rows).T
    pd.set_option("display.width", 200)
    log(t1.to_string(float_format=lambda v: f"{v:.2f}"))

    log(f"\n[2] 수익이 0인 가짜 팩터 {N_FAKE}개 (개월 수는 1963-07 ~ {END}와 같게), {N_SIM}번 반복")
    T = len(f.loc["1963-07":])
    rng = np.random.default_rng(SEED)
    c2, c3, mx = [], [], []
    for _ in range(N_SIM):
        r = rng.normal(0, 0.03, (T, N_FAKE))
        t = r.mean(0) / (r.std(0, ddof=1) / np.sqrt(T))
        c2.append((np.abs(t) > 2).sum())
        c3.append((np.abs(t) > 3).sum())
        mx.append(np.abs(t).max())
    c2, c3, mx = map(np.array, (c2, c3, mx))
    log(f"  개월 수 {T}. |t|>2 인 가짜 팩터: 평균 {c2.mean():.1f}개 (5~95% {np.percentile(c2, 5):.0f}~{np.percentile(c2, 95):.0f}개)")
    log(f"  |t|>3 인 가짜 팩터: 평균 {c3.mean():.2f}개, 하나 이상 나온 비율 {(c3 > 0).mean():.1%}")
    log(f"  가장 큰 |t|: 중앙값 {np.median(mx):.2f}, 5~95% {np.percentile(mx, 5):.2f}~{np.percentile(mx, 95):.2f}")
    log(f"  정규분포 기준 기대값: |t|>2 {N_FAKE * p_two_sided(2):.1f}개, |t|>3 {N_FAKE * p_two_sided(3):.2f}개 "
        f"(t=3의 양쪽 p값 {p_two_sided(3):.2%})")
    one = np.random.default_rng(SEED + 1).normal(0, 0.03, (T, N_FAKE))
    one_t = one.mean(0) / (one.std(0, ddof=1) / np.sqrt(T))
    log(f"  차트에 그린 한 번의 예: |t|>2 {(np.abs(one_t) > 2).sum()}개, |t|>3 {(np.abs(one_t) > 3).sum()}개, "
        f"가장 큰 |t| {np.abs(one_t).max():.2f}")

    log("\n[3] 검정 횟수별 본페로니 기준 (양쪽 5%)")
    for n in [1, 5, 10, 32, 100, 316]:
        log(f"  {n:>4}번: t > {bonferroni_cut(n):.2f}")

    log(f"\n[4] 25개 포트폴리오 CAPM 절편, 1963-07 ~ {END}")
    pp, ff = p.loc["1963-07":], f.loc["1963-07":]
    ex = pp.sub(ff["RF"], axis=0)
    X = np.column_stack([np.ones(len(ff)), ff["Mkt-RF"].values])
    ts = []
    for c in pp.columns:
        y = ex[c].values
        coef, *_ = np.linalg.lstsq(X, y, rcond=None)
        res = y - X @ coef
        se = np.sqrt(res @ res / (len(y) - 2) * np.linalg.inv(X.T @ X)[0, 0])
        ts.append(coef[0] / se)
    ts = np.array(ts)
    pv = np.array([p_two_sided(t) for t in ts])
    log(f"  |t|>2: {(np.abs(ts) > 2).sum()}개, |t|>3: {(np.abs(ts) > 3).sum()}개, "
        f"본페로니(t>{bonferroni_cut(25):.2f}): {(np.abs(ts) > bonferroni_cut(25)).sum()}개, "
        f"홀름: {holm(pv).sum()}개, BHY: {bhy(pv).sum()}개")
    log("  절편 t값(작은 규모부터, 각 줄은 성장주 → 가치주): ")
    for i in range(5):
        log("   " + " ".join(f"{v:6.2f}" for v in ts[i * 5:(i + 1) * 5]))

    with open(os.path.join(OUT, "factor_hurdle.txt"), "w", encoding="utf-8") as fh:
        fh.write("\n".join(lines) + "\n")
    plot(one_t, t1)


def plot(one_t, t1):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.4), gridspec_kw={"width_ratios": [1.2, 1]})
    order = np.argsort(one_t)
    x = np.arange(len(one_t))
    colors = np.where(np.abs(one_t[order]) > 3, "#b22222", np.where(np.abs(one_t[order]) > 2, "#1f5fbf", "#8c8c8c"))
    ax1.scatter(x, one_t[order], s=9, c=colors)
    for lv, st in [(2, "--"), (3, "-")]:
        for sgn in (1, -1):
            ax1.axhline(sgn * lv, color="#333333", linestyle=st, linewidth=0.8)
    ax1.text(len(x) * 0.02, 2.08, "|t| = 2", fontsize=8, color="#333333")
    ax1.text(len(x) * 0.02, 3.08, "|t| = 3", fontsize=8, color="#333333")
    n2, n3 = int((np.abs(one_t) > 2).sum()), int((np.abs(one_t) > 3).sum())
    ax1.set_title(f"{len(one_t)} fake factors with zero true return: {n2} pass |t| > 2, {n3} pass |t| > 3", fontsize=10)
    ax1.set_xlabel("Fake factors, sorted by t-stat")
    ax1.set_ylabel("t-stat of mean monthly return")
    ax1.set_ylim(-4, 4)

    names = list(t1.index)
    xs = np.arange(len(names))
    ax2.bar(xs - 0.19, t1["표본 t"].astype(float), width=0.36, color="#1f5fbf", label="Original paper's sample")
    ax2.bar(xs + 0.19, t1["그 뒤 t"].astype(float), width=0.36, color="#8c8c8c", label=f"After the sample, to {END}")
    for lv, st in [(2, "--"), (3, "-")]:
        ax2.axhline(lv, color="#333333", linestyle=st, linewidth=0.8)
    ax2.axhline(0, color="#333333", linewidth=0.6)
    ax2.set_xticks(xs, ["Market", "SMB", "HML", "Momentum"])
    ax2.set_ylabel("t-stat of mean monthly return")
    ax2.set_title("US factors: in the paper's sample vs after", fontsize=10)
    ax2.legend(frameon=False, fontsize=8, loc="upper right")
    for ax in (ax1, ax2):
        ax.spines[["top", "right"]].set_visible(False)
        ax.grid(axis="y", color="#e6e6e6", linewidth=0.6)
        ax.set_axisbelow(True)
    fig.tight_layout()
    fig.savefig(os.path.join(OUT, "64_factor_hurdle.png"), dpi=130)
    plt.close(fig)


if __name__ == "__main__":
    main()