# -*- coding: utf-8 -*-
"""
산업별 심층 분석 — 시나리오 편 (챕터 5 다음):
 이번엔 외생변수(유가·교역·금리·환율·노동)를 포함한다.
 (1) BVARX(내생 36산업 + 외생 5)에 DSGE 시나리오 외생경로(S1–S5)를 조건화 → 산업별 10년 수준 전망(시나리오 팬)
 (2) SHAP 설명모형에 외생 5종을 특성으로 추가 → 이제 거시요인도 동인으로 등장
 (3) 시나리오 스프레드 + SHAP 동인을 엮어 산업별 리포트 생성
산출: repo/industry-scenario/data/scenario_data.json
"""
import os, json
import numpy as np, pandas as pd
from sklearn.ensemble import GradientBoostingRegressor
import shap

HERE=os.path.dirname(os.path.abspath(__file__)); ROOT=os.path.dirname(HERE)
DATA=os.path.join(ROOT,"data"); os.makedirs(DATA,exist_ok=True); REPO=os.path.dirname(ROOT)
P,Q,LAM,LAMX,WIN,HH = 4,2,0.1,1.0,30.0,40
EXO=["oil_g","trd_g","d_rrate","rfx_g","lab_g"]; EXOK=["유가","교역","실질금리","실질환율","노동"]
SCN=["S1","S2","S3","S4","S5"]; SCN_KO={"S1":"기준","S2":"공급제약+충격","S3":"제약완화","S4":"기술도약","S5":"지정학"}

# ---------- 데이터 (내생∩외생 공통표본) ----------
rate=pd.read_csv(os.path.join(REPO,"ch1","data","dlog_rates_sa.csv"),encoding="utf-8-sig",index_col=0)
rate.index=pd.PeriodIndex(rate.index,freq="Q")
g=rate.drop(columns=["총량GDP"]).dropna().clip(-WIN,WIN)
ex=pd.read_csv(os.path.join(REPO,"ch3","data","exog_quarterly_filled.csv"),encoding="utf-8-sig")
ex.index=pd.PeriodIndex(ex["분기"],freq="Q"); ex=ex[EXO].astype(float)
common=g.index.intersection(ex.index).sort_values()
g=g.loc[common]; ex=ex.loc[common]
names=list(g.columns); n=len(names); ne=len(EXO)
gv=g.values; xv=ex.values; T=gv.shape[0]
def short(s): return (s.replace(" 제조업","").replace(" 및 ","·").replace("서비스업","서비스").replace(" 공급업","")).strip()[:12]
SH=[short(s) for s in names]

# ---------- BVARX 추정 (전표본) ----------
def ar_sigma(m,p):
    s=np.zeros(m.shape[1])
    for i in range(m.shape[1]):
        y=m[:,i]; Y=y[p:]; Z=np.column_stack([y[p-l:len(y)-l] for l in range(1,p+1)]+[np.ones(len(Y))])
        b,*_=np.linalg.lstsq(Z,Y,rcond=None); s[i]=(Y-Z@b).std(ddof=Z.shape[1])
    return s
def design(gvv,xvv,p,q):
    T0=max(p,q); TT=len(gvv)-T0; blocks=[]
    for l in range(1,p+1): blocks.append(gvv[[t-l for t in range(T0,len(gvv))]])
    for s in range(0,q+1): blocks.append(xvv[[t-s for t in range(T0,len(gvv))]])
    return gvv[T0:], np.column_stack(blocks+[np.ones(TT)])
sy=ar_sigma(gv,P); sx=xv.std(0); d=[]
for l in range(1,P+1):
    for j in range(n): d.append((l**2)*(sy[j]**2)/LAM**2)
for s in range(0,Q+1):
    for m in range(ne): d.append(((s+1)**2)*(sx[m]**2)/LAMX**2)
d.append(1e-6); d=np.array(d)
Y,X=design(gv,xv,P,Q); B=np.linalg.solve(X.T@X+np.diag(d),X.T@Y)

# ---------- DSGE 시나리오 외생경로 조건부 전망 ----------
dsge=pd.read_csv(os.path.join(REPO,"ch4","data","dsge_exog_paths.csv"),encoding="utf-8-sig")
def scen_forecast(scn):
    xs=dsge[dsge["시나리오"]==scn][EXO].values          # (40,5)
    exog_full=np.vstack([xv, xs]); start=len(xv)
    buf=[r.copy() for r in gv[-P:]]; out=[]
    for h in range(HH):
        ti=start+h; row=[buf[-l] for l in range(1,P+1)]+[exog_full[ti-sx2] for sx2 in range(Q+1)]
        y=np.concatenate(row+[[1.0]])@B; buf.append(y); out.append(y)
    return np.array(out)                                  # (40,n) 성장률
scen_growth={s:scen_forecast(s) for s in SCN}
def lvl(fc): return 100*np.exp(np.cumsum(fc/100,axis=0))
scen_lvl={s:lvl(scen_growth[s]) for s in SCN}             # (40,n) 수준지수
logcum=np.cumsum(gv/100,axis=0); histL=100*np.exp(logcum-logcum[-1]); HW=20

# ---------- SHAP (외생 포함) ----------
Xf=np.column_stack([gv[:-1], xv[1:]])                    # 36산업 시차1 + 외생5 동시
feat_names=list(names)+EXOK; feat_type=["ind"]*n+["exo"]*ne
Yf=gv[1:]
def direction(xcol,y):
    c=np.corrcoef(xcol,y)[0,1]; return "+" if c>=0 else "-"
records={}
for i,nm in enumerate(names):
    y=Yf[:,i]
    gbm=GradientBoostingRegressor(n_estimators=200,max_depth=3,learning_rate=0.05,subsample=0.9,random_state=0)
    gbm.fit(Xf,y)
    sv=shap.TreeExplainer(gbm).shap_values(Xf); imp=np.abs(sv).mean(0); imp=imp/(imp.sum()+1e-12)*100
    order=np.argsort(imp)[::-1]; own=feat_names.index(nm)
    drivers=[]
    for k in order:
        if k==own: continue
        drivers.append({"name":SH[k] if feat_type[k]=="ind" else feat_names[k],
                        "type":feat_type[k],"share":round(float(imp[k]),1),"dir":direction(Xf[:,k],y)})
        if len(drivers)>=6: break
    lev_by={s:round(float(scen_lvl[s][-1,i]),1) for s in SCN}
    records[nm]={"short":SH[i],"own_share":round(float(imp[own]),1),"drivers":drivers,
                 "lev2035":lev_by,"hist":[round(float(v),1) for v in histL[-HW:,i]],
                 "scen":{s:[round(float(v),1) for v in scen_lvl[s][:,i]] for s in SCN}}

# ---------- 리포트 ----------
EXO_CTX={"유가":{"+":"유가 상승이 순(順)으로 작용","-":"유가 상승이 비용 부담"},
 "교역":{"+":"세계 교역 확대에 수요가 민감","-":"교역 확대가 경쟁·대체로 작용"},
 "실질금리":{"+":"금리 상승기에 견조","-":"금리 상승에 취약(투자·내구재)"},
 "실질환율":{"+":"환율 상승(원화 약세)이 수출 경쟁력","-":"환율 상승이 수입비용 부담"},
 "노동":{"+":"취업자 증가에 민감한 노동집약","-":"고용보다 자본·기술이 주도"}}
def make_report(nm,r):
    lv=r["lev2035"]; vals=list(lv.values()); lo,hi=min(vals),max(vals)
    smin=min(lv,key=lv.get); smax=max(lv,key=lv.get); own=r["own_share"]; ds=r["drivers"]
    parts=[]
    for dd in ds[:3]:
        if dd["type"]=="exo":
            ctx=EXO_CTX.get(dd["name"],{}).get(dd["dir"],"반응"); parts.append(f"<b>{dd['name']}</b>({dd['dir']}, 기여 {dd['share']}%) — {ctx}")
        else:
            rel=("동반 상승" if dd["dir"]=="+" else "역행"); parts.append(f"<b>{dd['name']}</b>({dd['dir']}, 기여 {dd['share']}%) — 이 산업과 {rel}하는 연관 산업")
    exo_top=[dd for dd in ds if dd["type"]=="exo"][:1]
    macro=(f" 거시 요인 중에서는 <b>{exo_top[0]['name']}</b>({exo_top[0]['dir']})의 영향이 가장 크다." if exo_top else " 상위 동인은 대부분 연관 산업이다.")
    spread=(f"시나리오에 따라 <b>{lo:.0f}~{hi:.0f}</b> 사이(2025Q4=100)로, "
            f"<b>{smax} {SCN_KO[smax]}</b>에서 가장 높고 <b>{smin} {SCN_KO[smin]}</b>에서 가장 낮다")
    body=(f"<p><b>{nm}</b>의 10년 수준 전망은 {spread}. 격차 {hi-lo:.0f}p — 대외·인구 여건에 대한 이 산업의 민감도다.</p>"
          f"<p><b>가장 크게 움직이는 요인(SHAP, 외생 포함).</b> ① {parts[0]}; ② {parts[1]}; ③ {parts[2]}.{macro}</p>"
          f"<p>자기 관성 기여 {own}%. DSGE가 만든 외생 시나리오가 이 산업 경로를 어떻게 가르는지가 위 팬으로 나타난다.</p>")
    return body
for nm,r in records.items(): r["report"]=make_report(nm,r)

out={"industries":[{"name":nm,"short":records[nm]["short"]} for nm in names],
     "scenarios":[{"id":s,"ko":SCN_KO[s]} for s in SCN],"base":"2025Q4=100","data":records}
json.dump(out,open(os.path.join(DATA,"scenario_data.json"),"w",encoding="utf-8"),ensure_ascii=False,indent=1)
print("시나리오 산업분석 완료:",n,"산업 × 5 시나리오 | scenario_data.json")
exn=[x for x in names if "컴퓨터" in x][0]
print("예시 컴퓨터·전자 2035:",records[exn]["lev2035"],"| 상위동인",[d["name"] for d in records[exn]["drivers"][:3]])