# -*- coding: utf-8 -*-
"""
논문2 · 챕터 4 — 시나리오 전망 (S1–S5)
 DSGE 외생경로(선행연구와 동일 엔진)를 월별 70산업 VARX에 조건화해 2035년까지 전망한다.
 노동공급(인구 중·저·고위) × 대외환경의 5개 시나리오. 몬테카를로 팬으로 불확실성 표시.
 재현 데모: 6산업 VARX(외생 결합)에 DSGE 월별 외생을 먹여 조건부 전망(전체 70과 동일 원리).
 전체 70산업 결과는 저장소 06_scenario.py 산출(scenario_*.csv)을 제시.
산출: ch4/figures/(fig1_fan.png, fig2_labor.png, fig3_demo.png), ch4/data/chP4_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)
QMIDX=[(i*3+2) for i in range(104)]
EXO=["oil_g","trd_g","d_rrate","rfx_g","lab_g"]
SCN=["S1","S2","S3","S4","S5"]; SCN_KO={"S1":"기준","S2":"제약+위기","S3":"제약완화","S4":"기술도약","S5":"지정학"}
CLR={"S1":"#333","S2":"#c0392b","S3":"#1f4e79","S4":"#2e7d32","S5":"#b8892f"}

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, X):
    Tn,k=X.shape; rho_grid=np.r_[np.linspace(0,0.95,20),0.97,0.99]; Q=len(q_level); A=np.zeros((Q,Tn))
    for i,me in enumerate(QMIDX):
        for d in range(3):
            j=me-d
            if 0<=j<Tn: A[i,j]=1/3
    yq=np.asarray(q_level,float); vv=~np.isnan(yq); A=A[vv]; yq=yq[vv]; Q=len(yq); best=None
    for rho in rho_grid:
        lag=np.abs(np.subtract.outer(np.arange(Tn),np.arange(Tn))); 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); _,ld=np.linalg.slogdet(Om*s2); ll=-0.5*(Q*np.log(2*np.pi)+ld+Q)
        if best is None or ll>best[0]: best=(ll,rho,X@beta+Vu@A.T@Oi@r)
    return best[2]
def recover(qva, proxy):
    X=np.column_stack([np.ones(T),np.linspace(0,1,T),
        pd.Series(zscore_log(proxy)).interpolate(limit_direction="both").bfill().ffill().values])
    return chowlin_mean(np.log(np.asarray(qva,float)), X)

# ---- 6산업 월별 복원 + 성장률 ----
qv=pd.read_csv(os.path.join(DATA,"q_va_demo6.csv"),encoding="utf-8-sig")
mp=pd.read_csv(os.path.join(DATA,"m_prod_demo6.csv"),encoding="utf-8-sig")
labs=[c for c in qv.columns if c!="quarter"]
Mlvl=np.column_stack([recover(qv[c].values, mp[c].values) for c in labs])
Yg=100*np.diff(Mlvl,axis=0)                        # (T-1,6) 월별 성장, 2000-02~2025-12
w=qv[labs].iloc[-1].values; w=w/w.sum()            # 2025 VA 비중(총량 가중)

# ---- 월별 외생 히스토리 정렬 ----
ex=pd.read_csv(os.path.join(DATA,"exog_hist_m.csv"),encoding="utf-8-sig")
ex["date"]=pd.to_datetime(ex["date"]); ex=ex.set_index("date").reindex(MON)
X=ex[EXO].apply(pd.to_numeric,errors="coerce").interpolate(limit_direction="both").fillna(0).values
Xg=X[1:]                                           # y와 정렬(2000-02~)

# ---- VARX(p=2, 외생 시차0..1) 릿지 ----
P,Qx,LAM=2,1,50.0
def design(Y,Xe):
    n=Y.shape[1]; rows=[]; tgt=[]
    for t in range(P,len(Y)):
        r=[Y[t-l] for l in range(1,P+1)]+[Xe[t-s] for s in range(Qx+1)]+[[1.0]]
        rows.append(np.concatenate(r)); tgt.append(Y[t])
    return np.array(rows),np.array(tgt)
Z,Yt=design(Yg,Xg); n=Yg.shape[1]
D=np.ones(Z.shape[1]); D[-1]=0
B=np.linalg.solve(Z.T@Z+LAM*np.diag(D),Z.T@Yt)     # (P*n+(Qx+1)*ne+1, n)

# ---- DSGE 시나리오 외생 (분기→월, /3 반복) ----
dsge=pd.read_csv(os.path.join(DATA,"dsge_exog_paths.csv"),encoding="utf-8-sig")
Hf=120
def scen_exog(s):
    d=dsge[dsge["scenario"]==s][EXO].values          # (40,5) 분기
    return np.repeat(d/3.0,3,axis=0)[:Hf]            # (120,5) 월
def scen_forecast(s):
    xs=scen_exog(s); Xfull=np.vstack([Xg, xs]); start=len(Xg)
    buf=[r.copy() for r in Yg[-P:]]; out=[]
    for h in range(Hf):
        ti=start+h; row=[buf[-l] for l in range(1,P+1)]+[Xfull[ti-sx] for sx in range(Qx+1)]+[[1.0]]
        y=np.concatenate(row)@B; buf.append(y); out.append(y)
    return np.array(out)                              # (120,6) 월 성장

gdp_idx={}
for s in SCN:
    gr=scen_forecast(s); agg=(gr*w[None,:]).sum(1)     # 월 총량 성장%
    gdp_idx[s]=100*np.exp(np.cumsum(agg/100))
# 연평균 증가율(전년대비, 연말지수)
def yoy(idx):
    yl=idx[11::12]; return [100*(yl[k]/yl[k-1]-1) for k in range(1,len(yl))]

# ---- 실제 70산업 결과 ----
names=pd.read_csv(os.path.join(DATA,"industry_names.csv"),encoding="utf-8-sig").set_index("ind_id")["ind_name"].to_dict()
g=pd.read_csv(os.path.join(DATA,"scenario_gdp_growth.csv"),encoding="utf-8-sig")
le=pd.read_csv(os.path.join(DATA,"scenario_labor_effect.csv"),encoding="utf-8-sig"); le["name"]=le.ind_id.map(names).fillna(le.ind_id)
le2=le[~le.name.str.contains("잔여")].sort_values("labor_effect_S3_S2",ascending=False)

print("논문2·챕터4 — 시나리오 전망")
print("  [데모 6산업] 2035 총량지수(2025.12=100):")
for s in SCN: print(f"    {s} {SCN_KO[s]:8} {gdp_idx[s][-1]:.1f}")
print("  [전체 70] 2026 증가율:", g[g.year==2026].set_index("scenario")["g_point"].round(2).to_dict())
print("  [전체 70] 노동효과(S3-S2) 상위:")
for _,r in le2.head(4).iterrows(): print(f"    +{r.labor_effect_S3_S2:.1f}  {r['name']}")

# ---- fig1: 실제 70산업 시나리오 증가율 팬 ----
fig,ax=plt.subplots(figsize=(10,5))
s1=g[g.scenario=="S1"].set_index("year")
ax.fill_between(s1.index, s1.g_p5, s1.g_p95, color="#8899aa", alpha=.18, label="S1 90% 구간(몬테카를로)")
for s in SCN:
    d=g[g.scenario==s].set_index("year")
    ax.plot(d.index, d.g_point, color=CLR[s], lw=2, marker="o", ms=3, label=f"{s} {SCN_KO[s]}")
ax.set_title("그림 1. 시나리오별 총 실질부가가치 전년대비 증가율 (전체 70산업)", fontweight="bold")
ax.set_ylabel("전년대비 증가율 (%)"); ax.set_xlabel("연도"); ax.legend(fontsize=8.5, ncol=2); ax.grid(alpha=.25)
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig1_fan.png"),bbox_inches="tight"); plt.close(fig)

# ---- fig2: 실제 노동효과(S3-S2) 상위 산업 ----
fig,ax=plt.subplots(figsize=(8.5,5.2)); top=le2.head(12).iloc[::-1]
ax.barh(range(len(top)), top.labor_effect_S3_S2.values, color="#1f4e79")
ax.set_yticks(range(len(top))); ax.set_yticklabels(top.name.values, fontsize=8.5)
ax.set_title("그림 2. 노동공급 효과 = 고위(S3) - 저위(S2), 2035 누적 상위 산업", fontweight="bold")
ax.set_xlabel("2035 누적 부가가치 격차 (%p)"); ax.grid(alpha=.2, axis="x")
for i,v in enumerate(top.labor_effect_S3_S2.values): ax.text(v+0.1,i,f"+{v:.1f}",va="center",fontsize=8)
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig2_labor.png"),bbox_inches="tight"); plt.close(fig)

# ---- fig3: 재현 데모 — 6산업 조건부 시나리오 총량 경로 ----
fig,ax=plt.subplots(figsize=(9.5,4.8))
tf=pd.date_range("2026-01-01",periods=Hf,freq="MS")
for s in SCN: ax.plot(tf, gdp_idx[s], color=CLR[s], lw=1.8, label=f"{s} {SCN_KO[s]} ({gdp_idx[s][-1]:.0f})")
ax.set_title("그림 3. 재현 데모 — 6산업 VARX 조건부 시나리오 총량 경로 (2025.12=100)", fontweight="bold")
ax.set_ylabel("지수 (2025.12=100)"); ax.legend(fontsize=8.5); ax.grid(alpha=.25); ax.margins(x=.02)
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig3_demo.png"),bbox_inches="tight"); plt.close(fig)

json.dump({"demo_gdp2035_idx":{s:round(float(gdp_idx[s][-1]),1) for s in SCN},
           "full70_growth_2026":{s:round(float(g[(g.scenario==s)&(g.year==2026)].g_point.iloc[0]),2) for s in SCN},
           "full70_growth_2035":{s:round(float(g[(g.scenario==s)&(g.year==2035)].g_point.iloc[0]),2) for s in SCN},
           "labor_top":[[r["name"],round(float(r.labor_effect_S3_S2),1)] for _,r in le2.head(5).iterrows()]},
          open(os.path.join(DATA,"chP4_results.json"),"w",encoding="utf-8"),ensure_ascii=False,indent=2)
print("저장: ch4/figures/(fig1,fig2,fig3), ch4/data/chP4_results.json")