# -*- coding: utf-8 -*-
"""
산업별 심층 분석 (ch2–ch3 사이):
 (1) 36개 산업 BVAR(미네소타) 10년 수준 전망 (2025Q4=100)
 (2) 각 산업의 성장을 가장 크게 '움직이는' 요인을 SHAP으로 규명
     설명모형: GradientBoosting(타깃=산업 i 성장, 특성=36산업 시차1 성장 + 외생5종)
     → TreeExplainer SHAP 평균|기여|로 상위 동인, 상관부호로 방향(+/−)
 (3) SHAP 결과·전망을 엮어 산업별 리포트 텍스트 생성(정적 임베드)
산출: repo/industry/data/shap_data.json  (페이지가 이 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)   # industry/
DATA=os.path.join(ROOT,"data"); os.makedirs(DATA,exist_ok=True)
REPO=os.path.dirname(ROOT)
P, LAM, WIN, HH = 4, 0.1, 30.0, 40
EXO=["oil_g","trd_g","d_rrate","rfx_g","lab_g"]; EXOK=["유가","교역","실질금리","실질환율","노동"]

# ---------- 데이터 ----------
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)
# ※ 이 단계(챕터 2, BVARX 이전)에서는 외생변수 X를 쓰지 않는다 —
#   오직 36개 산업 부가가치(GDP)와 산업연관표(IO)만으로 전망·설명한다.
names=list(g.columns); n=len(names)
gv=g.values; T=gv.shape[0]                      # 36산업 부가가치 성장(전체 표본)
def short(s): return (s.replace(" 제조업","").replace(" 및 ","·").replace("서비스업","서비스").replace(" 공급업","")).strip()[:12]
SH=[short(s) for s in names]

# ---------- (1) BVAR 미네소타 10년 수준 전망 ----------
def build_XY(m,p):
    Y=m[p:]; TT=Y.shape[0]; X=np.ones((TT,n*p+1))
    for l in range(1,p+1): X[:,(l-1)*n:l*n]=m[p-l:-l]
    return Y,X
def ar_sigma(m,p):
    s=np.zeros(n)
    for i in range(n):
        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 fit_bvar(m,p):
    Y,X=build_XY(m,p); sig=ar_sigma(m,p); d=np.zeros(n*p+1)
    for l in range(1,p+1):
        for j in range(n): d[(l-1)*n+j]=(l**2)*(sig[j]**2)/(LAM**2)
    d[-1]=1e-6
    return np.linalg.solve(X.T@X+np.diag(d),X.T@Y)
Bb=fit_bvar(gv,P)
def fpath(Bc):
    buf=[r.copy() for r in gv[-P:]]; path=[]
    for _ in range(HH):
        x=np.ones(n*P+1)
        for l in range(1,P+1): x[(l-1)*n:l*n]=buf[-l]
        y=x@Bc; buf.append(y); path.append(y)
    return np.array(path)
proj=fpath(Bb)                                  # (40,n) 성장률
projL=100*np.exp(np.cumsum(proj/100,axis=0))    # 수준지수(2025Q4=100)
logcum=np.cumsum(gv/100,axis=0); histL=100*np.exp(logcum-logcum[-1])
HW=20
lev2035=projL[-1]; cagr=(lev2035/100)**(1/10)-1

# --- 산업연관표(I-O) prior BVAR (방정식별 릿지, 교차벌점 /= ω_ij) ---
Om=pd.read_csv(os.path.join(REPO,"ch2","data","io_prior_omega.csv"),encoding="utf-8-sig",index_col=0)
Om=Om.reindex(index=names,columns=names).values
off=~np.eye(n,dtype=bool); Om[off]=Om[off]/Om[off].mean(); np.fill_diagonal(Om,1.0)
def fit_bvar_io(m,p):
    Y,X=build_XY(m,p); XtX=X.T@X; XtY=X.T@Y; sig=ar_sigma(m,p)
    d0=np.zeros(n*p+1)
    for l in range(1,p+1):
        for j in range(n): d0[(l-1)*n+j]=(l**2)*(sig[j]**2)/(LAM**2)
    d0[-1]=1e-6
    B=np.zeros((n*p+1,n))
    for i in range(n):
        d=d0.copy()
        for l in range(1,p+1):
            for j in range(n):
                if j!=i: d[(l-1)*n+j]/=Om[i,j]     # 연계 강하면(ω>1) 벌점↓
        B[:,i]=np.linalg.solve(XtX+np.diag(d),XtY[:,i])
    return B
Bio=fit_bvar_io(gv,P)
proj_io=fpath(Bio); projL_io=100*np.exp(np.cumsum(proj_io/100,axis=0))
lev2035_io=projL_io[-1]; cagr_io=(lev2035_io/100)**(1/10)-1

# ---------- (2) 산업별 SHAP (설명모형) ----------
# 특성: 36산업 시차1 성장 + 외생5종(동시), 타깃: 산업 i 성장
Xf=gv[:-1]                                       # 36산업 시차1 성장만 (외생 X 없음)
feat_names=list(names)
feat_type=["ind"]*n
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)
    expl=shap.TreeExplainer(gbm); sv=expl.shap_values(Xf)     # (samples, feats)
    imp=np.abs(sv).mean(0)                                    # 평균 |기여|
    imp=imp/(imp.sum()+1e-12)*100                             # 기여 비중 %
    order=np.argsort(imp)[::-1]
    own=int(feat_names.index(nm)) if nm in feat_names else -1
    own_share=float(imp[own]) if own>=0 else 0.0
    # 상위 외부 동인(자기 제외) 5개
    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)>=5: break
    records[nm]={"short":SH[i],"level_2035":round(float(lev2035[i]),1),
                 "cagr":round(float(cagr[i]*100),2),
                 "level_2035_io":round(float(lev2035_io[i]),1),
                 "cagr_io":round(float(cagr_io[i]*100),2),
                 "own_share":round(own_share,1),"drivers":drivers,
                 "hist":[round(float(v),1) for v in histL[-HW:,i]],
                 "proj":[round(float(v),1) for v in projL[:,i]],
                 "proj_io":[round(float(v),1) for v in projL_io[:,i]]}

# ---------- (3) 리포트 생성 (Claude 작성 서술 로직) — 36산업 간 상호작용만 ----------
def make_report(nm, r):
    lev=r["level_2035"]; cg=r["cagr"]; own=r["own_share"]; ds=r["drivers"]
    lev_io=r["level_2035_io"]
    gap=abs(lev-lev_io); prior_note=("두 사전분포의 전망은 거의 일치한다" if gap<3 else f"산업연관표 prior가 {'높게' if lev_io>lev else '낮게'} 갈린다({gap:.0f} 차이)")
    trend=("완만한 성장" if cg<1 else "견조한 성장" if cg<3 else "빠른 성장")
    top=ds[:3]
    parts=[]
    for d in top:
        rel=("동반 상승" if d["dir"]=="+" else "역행")
        parts.append(f"<b>{d['name']}</b>({d['dir']}, 기여 {d['share']}%) — 이 산업과 {rel}하는 연관 산업")
    momentum=("자체 과거 성장(관성)의 설명력이 "
              +("높아 추세가 안정적이다" if own>=25 else "중간 수준이다" if own>=12 else "낮아 다른 산업에 더 민감하다")
              +f"(자기 기여 {own}%).")
    body=(f"<p><b>{nm}</b>의 향후 10년 실질부가가치는 <b>{lev:.0f}</b>(미네소타)·<b>{lev_io:.0f}</b>(산업연관표) 수준(2025Q4=100), "
          f"연평균 <b>{cg:.1f}%</b>의 {trend}으로 전망된다 — {prior_note}.</p>"
          f"<p><b>가장 크게 움직이는 요인(SHAP).</b> 외생 거시변수 없이 <b>36개 산업 간 상호작용</b>만으로 분해했다(거시·외생은 챕터 3 BVARX에서). "
          f"① {parts[0]}; ② {parts[1]}; ③ {parts[2]}.</p>"
          f"<p>{momentum} 전망 경로는 이들 연관 산업의 동학에 따라 위아래로 조정될 수 있다.</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],
     "base":"2025Q4=100","quarters_hist":[str(q) for q in g.index[-HW:]],
     "data":records}
json.dump(out,open(os.path.join(DATA,"shap_data.json"),"w",encoding="utf-8"),ensure_ascii=False,indent=1)
print("산업별 분석 완료:",n,"개 산업 | shap_data.json 저장")
print("예시(컴퓨터·전자):")
ex_nm=[x for x in names if "컴퓨터" in x][0]
print("  2035수준",records[ex_nm]["level_2035"],"| 상위동인",[d["name"] for d in records[ex_nm]["drivers"][:3]])