# -*- coding: utf-8 -*-
"""
논문2 · 챕터 2 — MECE 70산업 분할 (62 기준분류 + 신산업 8, carve-out)
 신산업은 '더하기'가 아니라 모산업에서 '떼어내기'(carve-out)다: 모산업 = 신산업 + 잔여.
 신산업 크기 = 모산업 × 시변비중(기준연도 비중 × 신산업/모산업 프록시 상대성장).
 자기완결: db 불필요, ch2/data 의 CSV만 사용. 원 구현은 저장소 10_mece_build.py·_va_shares.py.
산출: ch2/figures/(fig1_carveout.png, fig2_new8.png, fig3_shares.png), ch2/data/chP2_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)]

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, rho_grid=None):
    Tn,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,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_parent(qva, proxy):
    cols=[np.ones(T), np.linspace(0,1,T)]
    if proxy is not None:
        cols.append(pd.Series(zscore_log(proxy)).interpolate(limit_direction="both").fillna(0).values)
    m=chowlin_mean(np.log(np.asarray(qva,float)), np.column_stack(cols))
    return np.exp(m)                      # 월별 부모 수준

# 신산업 → (신프록시, 모프록시, 2020비중, 소스)
CFG={
 "반도체":("C261","C26",0.500),"이차전지":("C282","C28",0.220),"바이오의약":("C212","C20",0.173),
 "신재생에너지":(None,None,0.080),"AI·SW":("J582","J",0.444),"디지털플랫폼":("G4791","G",0.120),
 "콘텐츠":("J591","J",0.451),"디지털헬스케어":("Q869","Q86",0.050)}
HOST_KO={"반도체":"전자","이차전지":"전기장비","바이오의약":"화학","신재생에너지":"전기가스",
         "AI·SW":"정보서비스","디지털플랫폼":"도소매","콘텐츠":"출판영상","디지털헬스케어":"보건"}

hv=pd.read_csv(os.path.join(DATA,"host_va_q.csv"),encoding="utf-8-sig")
px=pd.read_csv(os.path.join(DATA,"proxies_m.csv"),encoding="utf-8-sig")
yr=MON.year.values; is2020=(yr==2020)

def tv_share(newk, hostk, s0):
    if newk is None: return np.full(T, s0)
    nv=pd.Series(px[newk].values).interpolate(limit_direction="both").bfill().ffill().values
    hv_=pd.Series(px[hostk].values).interpolate(limit_direction="both").bfill().ffill().values
    rel=(nv/np.nanmean(nv[is2020]))/(hv_/np.nanmean(hv_[is2020]))
    return np.clip(s0*rel, 0.01, 0.95)

new_series={}; new_growth={}; parent_cache={}; share_ts={}
for nm,(nk,hk,s0) in CFG.items():
    proxy=px[hk].values if hk else None
    parent=recover_parent(hv[nm].values, proxy); parent_cache[nm]=parent
    s=tv_share(nk,hk,s0); share_ts[nm]=s
    new=parent*s; new_series[nm]=new
    new_growth[nm]=float(np.nanmean(100*np.diff(np.log(new))))

# --- carve-out 항등 검증 (반도체) ---
p=parent_cache["반도체"]; s=share_ts["반도체"]; semi=p*s; resid=p*(1-s)
ident_err=float(np.nanmax(np.abs((semi+resid)-p)/p))

print("논문2·챕터2 — MECE 70산업 분할 (carve-out 재현)")
print(f"  carve-out 항등 최대오차 (반도체+잔여 vs 전자부모): {ident_err:.2e}")
print(f"  구조: 62 기준분류 + 8 신산업 + 8 잔여 = 70 (합=경제 보존)")
print(f"  {'신산업':12}{'모산업':8}{'2020비중':>8}{'월평균성장%':>11}")
for nm in CFG: print(f"  {nm:12}{HOST_KO[nm]:8}{CFG[nm][2]*100:>7.1f}%{new_growth[nm]:>11.2f}")

def idx2020(v): return 100*v/np.nanmean(v[is2020])
CLR=["#185FA5","#9a3b2e","#2e7d4f","#b8860b","#6a4c93","#0e7c86","#c0392b","#555"]

# --- fig1: 반도체 carve-out (전자 부모 = 반도체 + 전자잔여) — 공통기준(부모 2020) ---
fig,ax=plt.subplots(figsize=(10,4.6)); sel=MON>=pd.Timestamp("2008-01-01")
base=np.nanmean(p[is2020])                      # 부모 2020 = 100 공통기준(항등식 보존)
def bidx(v): return 100*v/base
ax.stackplot(MON[sel],[bidx(semi)[sel],bidx(resid)[sel]],
             labels=["반도체 (carve-out)","전자·컴퓨터·통신 잔여"],colors=["#185FA5","#c9d6e5"])
ax.plot(MON[sel], bidx(p)[sel], color="#111", lw=1.4, ls="--", label="전자 부모 = 반도체 + 잔여 (합)")
ax.set_title("그림 1. carve-out 항등 — 전자 부모 = 반도체 + 잔여 (2020=100)", fontweight="bold")
ax.set_ylabel("지수 (2020=100)"); ax.legend(fontsize=9, loc="upper left"); ax.grid(alpha=.25); ax.margins(x=.02)
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig1_carveout.png"),bbox_inches="tight"); plt.close(fig)

# --- fig2: 신규 8개 산업 월별 부가가치 지수 (2020=100) ---
fig,ax=plt.subplots(figsize=(10,5)); sel=MON>=pd.Timestamp("2005-01-01")
dash=["-","--","-","--","-.","--","-",":"]
for i,nm in enumerate(CFG):
    ax.plot(MON[sel], idx2020(new_series[nm])[sel], color=CLR[i], ls=dash[i], lw=1.6, label=nm)
ax.set_title("그림 2. MECE 분할로 복원한 신규 8개 산업 월별 실질부가가치 (2020=100)", fontweight="bold")
ax.set_ylabel("지수 (2020=100)"); ax.legend(fontsize=8.5, ncol=2); ax.grid(alpha=.25); ax.margins(x=.02)
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig2_new8.png"),bbox_inches="tight"); plt.close(fig)

# --- fig3: 2020 분할비중 ---
fig,ax=plt.subplots(figsize=(8,4.2))
nms=list(CFG); sh=[CFG[n][2]*100 for n in nms]; order=np.argsort(sh)
ax.barh([nms[i]+f" ({HOST_KO[nms[i]]})" for i in order],[sh[i] for i in order],color="#185FA5")
ax.set_title("그림 3. 신규 8개 산업의 2020 모산업 내 분할비중", fontweight="bold")
ax.set_xlabel("모산업 내 비중 (%)"); ax.grid(alpha=.2, axis="x")
for i,idx in enumerate(order): ax.text(sh[idx]+0.5,i,f"{sh[idx]:.0f}%",va="center",fontsize=8.5)
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig3_shares.png"),bbox_inches="tight"); plt.close(fig)

json.dump({"ident_err_semi":ident_err,"structure":"54 base + 8 emerging + 8 residual = 70",
           "new8":[{"name":nm,"host":HOST_KO[nm],"share_2020_pct":round(CFG[nm][2]*100,1),
                    "mean_growth_pct":round(new_growth[nm],2)} for nm in CFG]},
          open(os.path.join(DATA,"chP2_results.json"),"w",encoding="utf-8"),ensure_ascii=False,indent=2)
print("저장: ch2/figures/(fig1,fig2,fig3), ch2/data/chP2_results.json")