# -*- coding: utf-8 -*-
"""
ENB(유효 베팅 수) 실측 + 견고성 검증

  pip install pandas numpy        # 그 외 의존성 없음. data/ 폴더가 이 파일 옆에 있어야 한다.
  python enb_analyze.py           # → results/enb_results.json
  · 유니버스 3종: KOSPI 상위30 / 업비트 상위20 / 멀티에셋 8
  · 연도별, 윈저화, 회전방식(공분산PCA vs 상관PCA) 3중 검증
"""
import os, sys, json
import numpy as np, pandas as pd
sys.stdout.reconfigure(encoding='utf-8', errors='replace')

HERE = os.path.dirname(os.path.abspath(__file__))
DATA, OUT = os.path.join(HERE, "data"), os.path.join(HERE, "results")
os.makedirs(OUT, exist_ok=True)


# ───────────────────────────────────────────────────────── 유틸
def load(fn, names_fn=None):
    px = pd.read_csv(os.path.join(DATA, fn), index_col=0, encoding="utf-8-sig")
    px.index = pd.to_datetime(px.index)
    if names_fn and os.path.exists(os.path.join(DATA, names_fn)):
        px = px.rename(columns=json.load(open(os.path.join(DATA, names_fn), encoding="utf-8")))
    return px.sort_index()


def logret(px, min_obs_ratio=0.9):
    r = np.log(px / px.shift(1))
    keep = r.columns[r.notna().sum() >= min_obs_ratio * len(r)]
    return r[keep].dropna()


def winsorize(r, lo=0.01, hi=0.99):
    return r.clip(r.quantile(lo), r.quantile(hi), axis=1)


def _entropy_enb(c):
    c = np.clip(c, 0, None)
    p = c / c.sum()
    nz = p[p > 1e-15]
    return float(np.exp(-(nz * np.log(nz)).sum())), p


def enb_cov(S, w):
    """Meucci ENB — 공분산 주성분 회전"""
    vals, vecs = np.linalg.eigh(S)
    o = np.argsort(vals)[::-1]
    vals, vecs = vals[o], vecs[:, o]
    wt = vecs.T @ w
    return _entropy_enb((wt ** 2) * vals)


def enb_corr(C, w, sig):
    """상관 주성분 회전 (변동성 차이를 요인 정의에서 제거)"""
    v = w * sig                        # 위험 스케일 가중
    vals, vecs = np.linalg.eigh(C)
    o = np.argsort(vals)[::-1]
    vals, vecs = vals[o], vecs[:, o]
    vt = vecs.T @ v
    return _entropy_enb((vt ** 2) * vals)


def metrics(r, ann):
    N, T = r.shape[1], r.shape[0]
    w = np.ones(N) / N
    S = np.cov(r.values.T) * ann
    C = np.corrcoef(r.values.T)
    sig = np.sqrt(np.diag(S))
    sp = float(np.sqrt(w @ S @ w))
    off = C[~np.eye(N, dtype=bool)]

    e_cov, p_cov = enb_cov(S, w)
    e_cor, _ = enb_corr(C, w, sig)

    cvals = np.sort(np.linalg.eigvalsh(C))[::-1]
    q = N / T
    lam_p = (1 + np.sqrt(q)) ** 2

    pr = r.values @ w
    m = pr <= np.quantile(pr, 0.10)
    rho_down = float(np.corrcoef(r.values[m].T)[~np.eye(N, dtype=bool)].mean()) if m.sum() > N else np.nan

    mrc = S @ w
    rc = w * mrc / sp
    return dict(N=N, T=T, rho_bar=float(off.mean()), rho_down=rho_down,
                sigma_mean=float(sig.mean()), sigma_port=sp,
                ENC=float(1 / (w ** 2).sum()), DR=float((w @ sig) / sp),
                ENB=e_cov, ENB_corr=e_cor,
                pc1_share=float(p_cov[0]), pc_max_share=float(p_cov.max()),
                mp_lambda_plus=float(lam_p),
                n_eig_above_mp=int((cvals > lam_p).sum()),
                top_rc=float((rc / rc.sum()).max()),
                eigvals=[float(x) for x in cvals])


SPECS = [("KOSPI 상위 30종목", "kospi30_px.csv", "kospi30_names.json", 252),
         ("암호화폐 상위 20 (업비트)", "upbit20_px.csv", None, 365),
         ("멀티에셋 8자산 (ETF)", "multi_px.csv", "multi_names.json", 252)]

ALL = {}
for name, fn, nfn, ann in SPECS:
    px = load(fn, nfn)
    r_all = logret(px)
    full = metrics(r_all, ann)
    full_w = metrics(winsorize(r_all), ann)

    # 연도별
    by_year = {}
    for y, g in r_all.groupby(r_all.index.year):
        gg = g.dropna(axis=1, how="any")
        if len(gg) < 100 or gg.shape[1] < 4:
            continue
        by_year[int(y)] = metrics(gg, ann)

    ALL[name] = dict(period="%s ~ %s" % (r_all.index.min().date(), r_all.index.max().date()),
                     full=full, winsorized=full_w, by_year=by_year,
                     assets=list(r_all.columns))

    print("\n" + "=" * 82)
    print("  %s    N=%d  T=%d  (%s)" % (name, full["N"], full["T"], ALL[name]["period"]))
    print("=" * 82)
    print("  평균 쌍상관 ρ̄        %6.3f    (하락 상위10%% 구간 %6.3f)" % (full["rho_bar"], full["rho_down"]))
    print("  개별 연변동성 평균    %6.1f%%   등가중 포트 %5.1f%%" % (full["sigma_mean"]*100, full["sigma_port"]*100))
    print("  ENC (비중만)         %6.1f" % full["ENC"])
    print("  DR  (분산화 비율)     %6.2f" % full["DR"])
    print("  " + "-" * 78)
    print("  ENB  공분산PCA       %6.2f   |  윈저화 후 %5.2f  |  상관PCA %5.2f"
          % (full["ENB"], full_w["ENB"], full["ENB_corr"]))
    print("  최대 PC의 위험비중    %5.1f%%" % (full["pc_max_share"]*100))
    print("  MP상한 초과 고유값    %d개 / %d개" % (full["n_eig_above_mp"], full["N"]))
    if by_year:
        ys = sorted(by_year)
        print("  " + "-" * 78)
        print("  연도별 ENB(공분산PCA):")
        print("    " + "  ".join("%4d" % y for y in ys))
        print("    " + "  ".join("%4.2f" % by_year[y]["ENB"] for y in ys))
        print("  연도별 ρ̄:")
        print("    " + "  ".join("%4.2f" % by_year[y]["rho_bar"] for y in ys))
        e = [by_year[y]["ENB"] for y in ys]
        print("  → 연도별 ENB 범위 %.2f ~ %.2f (중앙값 %.2f)" % (min(e), max(e), float(np.median(e))))

json.dump(ALL, open(os.path.join(OUT, "enb_results.json"), "w", encoding="utf-8"),
          ensure_ascii=False, indent=1)

print("\n\n" + "#" * 82)
print("  최종 요약 — 등가중 포트폴리오의 유효 베팅 수")
print("#" * 82)
print("  %-26s %4s %6s %7s %7s %8s %8s" % ("유니버스", "N", "ρ̄", "ENC", "ENB", "윈저화", "연도범위"))
for k, v in ALL.items():
    ys = sorted(v["by_year"])
    e = [v["by_year"][y]["ENB"] for y in ys]
    rng = "%.1f~%.1f" % (min(e), max(e)) if e else "-"
    print("  %-26s %4d %6.2f %7.1f %7.2f %8.2f %10s"
          % (k[:26], v["full"]["N"], v["full"]["rho_bar"], v["full"]["ENC"],
             v["full"]["ENB"], v["winsorized"]["ENB"], rng))
print("\n→ results/enb_results.json")