# -*- coding: utf-8 -*-
"""
논문2 · 챕터 1 — 혼합주기 상태공간(Mixed-frequency State-Space)
 분기(3개월에 한 번) 실질부가가치를 월별 지표(생산지수)의 '모양'으로 월 단위로 복원한다.
 방법: 누적기 상태공간의 닫힌형 해 = Chow-Lin BLUE(월별 로그수준의 분기 평균 = 분기 로그부가가치).
 자기완결: db 불필요, ch1/data 의 CSV만 사용. 원 구현은 저장소 _mf_kalman.py.
산출: ch1/figures/(fig1_recover.png, fig2_fit.png, fig3_rho.png), ch1/data/chP1_results.json
"""
import os, json
import numpy as np, pandas as pd
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
plt.rcParams["font.family"]="Malgun Gothic"; plt.rcParams["axes.unicode_minus"]=False
plt.rcParams["savefig.dpi"]=130
HERE=os.path.dirname(os.path.abspath(__file__)); CH=os.path.dirname(HERE)
DATA=os.path.join(CH,"data"); FIG=os.path.join(CH,"figures"); os.makedirs(FIG,exist_ok=True)

MON=pd.date_range("2000-01-01","2025-12-01",freq="MS"); T=len(MON)
QS=[f"{2000+i//4}Q{i%4+1}" for i in range(104)]
QMIDX=[(i*3+2) for i in range(104)]              # 각 분기의 마지막 월 인덱스(0-based)

def zscore_log(v):
    v=np.asarray(v,float); out=np.full_like(v,np.nan); m=v>0
    lg=np.log(v[m]); out[m]=(lg-lg.mean())/(lg.std()+1e-9); return out

def chowlin_mean(q_level, qmidx, X, rho_grid=None):
    """월별 로그수준 m_t 복원. 제약: 분기 로그부가가치 = 3개월 m_t의 '평균'(정합 재구성오차≈0)."""
    T,k=X.shape
    if rho_grid is None: rho_grid=np.r_[np.linspace(0,0.95,20),0.97,0.99]
    Q=len(q_level); A=np.zeros((Q,T))
    for i,me in enumerate(qmidx):
        for d in range(3):
            j=me-d
            if 0<=j<T: A[i,j]=1.0/3.0                 # 평균 집계
    yq=np.asarray(q_level,float); valid=~np.isnan(yq)
    A=A[valid]; yq=yq[valid]; Q=len(yq)
    best=None
    for rho in rho_grid:
        lag=np.abs(np.subtract.outer(np.arange(T),np.arange(T)))
        Vu=rho**lag/(1-rho**2+1e-12)
        AX=A@X; Om=A@Vu@A.T+1e-8*np.eye(Q); Oi=np.linalg.pinv(Om)
        beta=np.linalg.solve(AX.T@Oi@AX+1e-8*np.eye(k), AX.T@Oi@yq)
        r=yq-AX@beta; s2=float(r@Oi@r)/max(Q-k,1)
        _,logdet=np.linalg.slogdet(Om*s2); ll=-0.5*(Q*np.log(2*np.pi)+logdet+Q)
        if best is None or ll>best[0]:
            u=Vu@A.T@Oi@r; best=(ll,rho,beta,X@beta+u,s2)
    ll,rho,beta,m,s2=best
    return dict(m=m,rho=rho,loglik=ll,s2=s2)

def build_X(prod):
    cols=[np.ones(T), np.linspace(0,1,T)]
    p=pd.Series(zscore_log(prod)).interpolate(limit_direction="both").fillna(0).values
    cols.append(p); return np.column_stack(cols)

def recover(qva, prod, mask_last=0):
    ylog=np.log(np.asarray(qva,float)); qy=ylog.copy(); masked=None
    if mask_last>0: masked=qy[-mask_last:].copy(); qy[-mask_last:]=np.nan
    res=chowlin_mean(qy, QMIDX, build_X(prod)); res["masked"]=masked; res["qlog"]=ylog
    return res

# --- 입력 ---
qv=pd.read_csv(os.path.join(DATA,"q_va_demo.csv"),encoding="utf-8-sig")
mp=pd.read_csv(os.path.join(DATA,"m_prod_demo.csv"),encoding="utf-8-sig")
IND={"elec":"컴퓨터·전자·광학","auto":"자동차·운송장비"}
res={}; recon={}; hold={}
for key in ["elec","auto"]:
    r=recover(qv[key].values, mp[key].values); res[key]=r
    qfit=pd.Series(r["m"]).groupby(np.repeat(np.arange(104),3)).mean().values   # 분기평균
    recon[key]=float(np.nanmax(np.abs(qfit-r["qlog"])))
    h=recover(qv[key].values, mp[key].values, mask_last=8)
    qf2=pd.Series(h["m"]).groupby(np.repeat(np.arange(104),3)).mean().values
    pred=qf2[-8:]; act=h["masked"]
    hold[key]=dict(rmse=float(np.sqrt(np.nanmean((pred-act)**2))),
                   corr=float(np.corrcoef(pred,act)[0,1]))

diag=pd.read_csv(os.path.join(DATA,"latent_diagnostics.csv"),encoding="utf-8-sig")

print("논문2·챕터1 — 혼합주기 상태공간 (재현 데모)")
for key in ["elec","auto"]:
    print(f"  {IND[key]:14} rho={res[key]['rho']:.3f} | 재구성오차={recon[key]:.2e} | holdout corr={hold[key]['corr']:.3f}")
print(f"  [70산업 전체] rho 평균={diag.rho.mean():.3f}[{diag.rho.min():.2f},{diag.rho.max():.2f}] | 재구성오차 최대={diag.max_recon_err.max():.2e}")

# 지수화(2000.01=100) 도우미
def idx100(m): return 100*np.exp(m-m[0])
tm=MON

# --- fig1: 전자 월별 복원 (2015~2025 확대) ---
fig,ax=plt.subplots(figsize=(10,4.6))
m=res["elec"]["m"]; lvl=idx100(m)
proxy=zscore_log(mp["elec"].values); proxy_idx=100*np.exp((proxy-np.nanmean(proxy))*0.15+ (m.mean()-m[0]))  # 배경 모양만
# 분기 앵커: 각 분기의 월평균 로그를 지수화한 값을 분기 3개월에 계단으로
qlog=res["elec"]["qlog"]; qanchor=100*np.exp(qlog-m[0])
sel=(tm>=pd.Timestamp("2015-01-01"))
ax.plot(tm[sel], lvl[sel], color="#185FA5", lw=2.0, label="복원된 월별 부가가치", zorder=3)
# 분기 앵커 점(각 분기 마지막 월 위치)
qx=[tm[QMIDX[i]] for i in range(104)]
qsel=[i for i in range(104) if qx[i]>=pd.Timestamp("2015-01-01")]
ax.scatter([qx[i] for i in qsel],[qanchor[i] for i in qsel],color="#2e7d4f",s=22,zorder=4,label="분기 공식치(앵커·평균)")
# 생산지수 프록시(표준화 모양, 우축)
ax2=ax.twinx(); ax2.plot(tm[sel], pd.Series(mp["elec"].values)[sel], color="#9a3b2e", lw=1.0, ls="--", alpha=.6, label="생산지수 프록시(우축)")
ax2.set_ylabel("생산지수", color="#9a3b2e", fontsize=9); ax2.tick_params(labelcolor="#9a3b2e")
ax.set_title("그림 1. 분기 부가가치를 월별로 복원 — 컴퓨터·전자·광학 (2015~2025)", fontweight="bold")
ax.set_ylabel("지수 (2000.01=100)"); ax.grid(alpha=.25)
l1,la1=ax.get_legend_handles_labels(); l2,la2=ax2.get_legend_handles_labels()
ax.legend(l1+l2, la1+la2, fontsize=8.5, loc="upper left")
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig1_recover.png"),bbox_inches="tight"); plt.close(fig)

# --- fig2: 정합성 + holdout ---
fig,ax=plt.subplots(1,2,figsize=(12,4.4))
# (a) full-sample 전자 월별 vs 분기평균
m=res["elec"]["m"]; ax[0].plot(tm, idx100(m), color="#185FA5", lw=1.3, label="월별 복원")
qx=[tm[QMIDX[i]] for i in range(104)]
ax[0].scatter(qx, 100*np.exp(res["elec"]["qlog"]-m[0]), color="#2e7d4f", s=9, label="분기 공식치")
ax[0].set_title("(a) 전 구간 정합 — 전자 (2000~2025)"); ax[0].set_ylabel("지수(2000.01=100)")
ax[0].legend(fontsize=8.5); ax[0].grid(alpha=.25)
# (b) holdout: 마지막 8분기 가리고 복원 예측 vs 실제
for key,cl in [("elec","#185FA5"),("auto","#9a3b2e")]:
    h=recover(qv[key].values, mp[key].values, mask_last=8)
    qf2=pd.Series(h["m"]).groupby(np.repeat(np.arange(104),3)).mean().values
    ax[1].plot(range(8), h["masked"], color=cl, lw=1.6, marker="o", ms=4, label=f"{IND[key]} 실제")
    ax[1].plot(range(8), qf2[-8:], color=cl, lw=1.2, ls="--", marker="x", ms=5, label=f"{IND[key]} 복원(corr={hold[key]['corr']:.2f})")
ax[1].set_title("(b) 보류표본(마지막 8분기) 검증 — 로그부가가치")
ax[1].set_xlabel("보류 분기 (0=2024Q1 … 7=2025Q4)"); ax[1].legend(fontsize=8); ax[1].grid(alpha=.25)
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig2_fit.png"),bbox_inches="tight"); plt.close(fig)

# --- fig3: 70산업 rho 분포 ---
fig,ax=plt.subplots(figsize=(8,3.8))
ax.hist(diag.rho.values, bins=15, color="#185FA5", edgecolor="white")
ax.axvline(diag.rho.mean(), color="#9a3b2e", lw=1.5, ls="--", label=f"평균 {diag.rho.mean():.2f}")
ax.set_title("그림 3. 70개 산업의 잔차 지속성 rho 분포 (재구성오차 최대 %.0e)"%diag.max_recon_err.max(), fontweight="bold")
ax.set_xlabel("rho (월별 잠재 잔차 AR(1) 계수)"); ax.set_ylabel("산업 수"); ax.legend(fontsize=9); ax.grid(alpha=.2)
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig3_rho.png"),bbox_inches="tight"); plt.close(fig)

json.dump({"demo":{k:{"rho":round(res[k]["rho"],3),"recon_err":recon[k],
                      "holdout_corr":round(hold[k]["corr"],3),"holdout_rmse":round(hold[k]["rmse"],4)} for k in ["elec","auto"]},
           "all70":{"rho_mean":round(float(diag.rho.mean()),3),"rho_min":round(float(diag.rho.min()),3),
                    "rho_max":round(float(diag.rho.max()),3),"recon_max":float(diag.max_recon_err.max()),"n":int(len(diag))}},
          open(os.path.join(DATA,"chP1_results.json"),"w",encoding="utf-8"),ensure_ascii=False,indent=2)
print("저장: ch1/figures/(fig1,fig2,fig3), ch1/data/chP1_results.json")