# -*- coding: utf-8 -*-
"""
논문2 · 챕터 3 — 빈도주의 VARX & Diebold-Yılmaz 연결성
 월별 70산업 부가가치 성장률의 상호작용을 정칙화(릿지) VARX로 추정하고,
 일반화 예측오차분산분해(Pesaran-Shin)로 '누가 충격을 발신/수신하는가'(연결성)를 측정한다.
 재현 데모: 6개 제조업을 SSM으로 월별 복원 → 릿지 VARX → 6x6 연결성(전체 70산업 방법과 동일).
 전체 70산업 결과는 저장소 05_varx_estimate.py 산출(connectedness.csv 등)을 제시.
산출: ch3/figures/(fig1_net.png, fig2_multipliers.png, fig3_demo.png), ch3/data/chP3_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)]

# ---------- SSM (챕터1과 동일, 평균집계 Chow-Lin) ----------
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)

# ---------- VARX(p) 릿지 + 동반행렬 고유값 + GFEVD 연결성 ----------
def var_ridge(Y, p=2, lam=1.0):
    Tn,n=Y.shape; Z=[];
    for t in range(p,Tn): Z.append(np.concatenate([Y[t-l] for l in range(1,p+1)]+[[1.0]]))
    Z=np.array(Z); Yt=Y[p:]
    D=np.ones(Z.shape[1]); D[-1]=0.0                      # 절편 비축소
    B=np.linalg.solve(Z.T@Z+lam*np.diag(D), Z.T@Yt)       # (np+1, n)
    resid=Yt-Z@B; Sig=np.cov(resid.T)
    Blags=[B[l*n:(l+1)*n].T for l in range(p)]            # 각 (n,n)
    return Blags, Sig, B
def companion_eig(Blags):
    n=Blags[0].shape[0]; p=len(Blags)
    C=np.zeros((n*p,n*p)); C[:n]=np.hstack(Blags)
    if p>1: C[n:]=np.eye(n*(p-1),n*p)
    return np.max(np.abs(np.linalg.eigvals(C)))
def gfevd_connectedness(Blags, Sig, H=12):
    n=Blags[0].shape[0]; p=len(Blags)
    Psi=[np.eye(n)]
    for h in range(1,H+1):
        s=np.zeros((n,n))
        for l in range(1,p+1):
            if h-l>=0: s+=Blags[l-1]@Psi[h-l]
        Psi.append(s)
    sig=np.diag(Sig); theta=np.zeros((n,n))
    for i in range(n):
        denom=sum((Psi[h][i,:]@Sig@Psi[h][i,:].T) for h in range(H+1))
        for j in range(n):
            num=sum((Psi[h][i,:]@Sig[:,j])**2 for h in range(H+1))/sig[j]
            theta[i,j]=num/denom
    theta=theta/theta.sum(1,keepdims=True)                # 행 정규화
    FROM=theta.sum(1)-np.diag(theta); TO=theta.sum(0)-np.diag(theta)
    NET=TO-FROM; TCI=100*FROM.mean()
    return theta*100, FROM*100, TO*100, NET*100, TCI

# ---------- 재현 데모: 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"]
M=np.column_stack([recover(qv[c].values, mp[c].values) for c in labs])   # (T,6) 로그수준
Yg=100*np.diff(M,axis=0)                                                  # 월별 성장률
Bl,Sig,_=var_ridge(Yg,p=2,lam=1.0)
eig_ridge=companion_eig(Bl)
Bl_ols,Sig_ols,_=var_ridge(Yg,p=2,lam=1e-9)
eig_ols=companion_eig(Bl_ols)
theta,FROM,TO,NET,TCI=gfevd_connectedness(Bl,Sig,H=12)

print("논문2·챕터3 — VARX & 연결성 (재현 데모 6산업)")
print(f"  동반행렬 max|eig|: 릿지={eig_ridge:.3f}  OLS={eig_ols:.3f}  (안정=1 미만)")
print(f"  데모 총연결성 TCI={TCI:.1f}%")
for i,l in enumerate(sorted(range(len(labs)),key=lambda k:-NET[k])):
    print(f"    NET {NET[l]:+6.1f}  {labs[l]}")

# 전체 70산업 실제 결과
summ=pd.read_csv(os.path.join(DATA,"varx_summary.csv"),encoding="utf-8-sig").iloc[0]
conn=pd.read_csv(os.path.join(DATA,"connectedness.csv"),encoding="utf-8-sig")
names=pd.read_csv(os.path.join(DATA,"industry_names.csv"),encoding="utf-8-sig").set_index("ind_id")["ind_name"].to_dict()
conn["name"]=conn.ind_id.map(names).fillna(conn.ind_id)
print(f"  [전체 70산업] T={int(summ['T'])} n={int(summ['n'])} p={int(summ['p'])} lam={summ['lambda']} max|eig|={summ['max_eig']} R2={summ['mean_R2']} TCI={summ['total_connectedness']}%")

# ---------- fig1: 70산업 순연결성(NET) 상·하위 ----------
cs=conn.sort_values("NET"); top=cs.tail(8); bot=cs.head(6); rows=pd.concat([bot,top])
fig,ax=plt.subplots(figsize=(8.5,6.5))
ax.barh(range(len(rows)),rows.NET.values,color=["#9a3b2e" if v<0 else "#185FA5" for v in rows.NET])
ax.axvline(0,color="#333",lw=.8); ax.set_yticks(range(len(rows))); ax.set_yticklabels(rows.name.values,fontsize=8)
ax.set_title(f"그림 1. 70산업 순 연결성(NET) 상·하위 — 총연결성 TCI {summ['total_connectedness']}%", fontweight="bold")
ax.set_xlabel("순 연결성 NET (%p) : 양(+)=발신 / 음(-)=수신")
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig1_net.png"),bbox_inches="tight"); plt.close(fig)

# ---------- fig2: 외생충격 동태승수 히트맵 (대표 산업) ----------
em=pd.read_csv(os.path.join(DATA,"varx_exog_multipliers.csv"),encoding="utf-8-sig").rename(columns={"Unnamed: 0":"ind_id"})
em["name"]=em.ind_id.map(names).fillna(em.ind_id)
pick=["MFG_COKE","MFG_CHEM","E_SEMI","MFG_AUTO","MFG_METAL","MFG_MACH","MFG_TEX","CON_CIVIL","WHOLE","FIN","TRANS","AGR"]
sub=em[em.ind_id.isin(pick)].set_index("name")[["oil_g","trd_g","d_rrate","rfx_g","gpr_l","lab_g"]]
sub=sub.reindex([names.get(p,p) for p in pick if names.get(p,p) in sub.index])
XL={"oil_g":"유가","trd_g":"교역","d_rrate":"실질금리","rfx_g":"실질환율","gpr_l":"지정학","lab_g":"고용"}
fig,ax=plt.subplots(figsize=(7.5,6.2))
V=np.clip(sub.values,-1,1)
im=ax.imshow(V,cmap="RdBu_r",vmin=-1,vmax=1,aspect="auto")
ax.set_xticks(range(6)); ax.set_xticklabels([XL[c] for c in sub.columns]); ax.set_yticks(range(len(sub))); ax.set_yticklabels(sub.index,fontsize=8)
for i in range(len(sub)):
    for j in range(6): ax.text(j,i,f"{sub.values[i,j]:+.2f}",ha="center",va="center",fontsize=7,color="#222")
ax.set_title("그림 2. 외생충격 누적 동태승수 (대표 산업 × 6 외생)", fontweight="bold")
fig.colorbar(im,ax=ax,fraction=.046,pad=.04,label="누적 반응 (색 clip ±1)")
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig2_multipliers.png"),bbox_inches="tight"); plt.close(fig)

# ---------- fig3: 재현 데모 6x6 연결성 히트맵 ----------
fig,ax=plt.subplots(figsize=(7,5.8))
im=ax.imshow(theta,cmap="Blues",aspect="auto")
ax.set_xticks(range(len(labs))); ax.set_xticklabels(labs,rotation=35,ha="right",fontsize=8.5)
ax.set_yticks(range(len(labs))); ax.set_yticklabels(labs,fontsize=8.5)
for i in range(len(labs)):
    for j in range(len(labs)):
        ax.text(j,i,f"{theta[i,j]:.0f}",ha="center",va="center",fontsize=7.5,
                color="white" if theta[i,j]>theta.max()*0.55 else "#222")
ax.set_xlabel("← 로부터(발신 산업)"); ax.set_ylabel("수신 산업")
ax.set_title(f"그림 3. 재현 데모 — 6산업 연결성 행렬(%) · TCI {TCI:.0f}% · max|eig| {eig_ridge:.2f}", fontweight="bold", fontsize=11)
fig.colorbar(im,ax=ax,fraction=.046,pad=.04,label="분산기여 (%)")
fig.tight_layout(); fig.savefig(os.path.join(FIG,"fig3_demo.png"),bbox_inches="tight"); plt.close(fig)

json.dump({"demo":{"tci":round(TCI,1),"eig_ridge":round(eig_ridge,3),"eig_ols":round(eig_ols,3),
                   "net":{labs[i]:round(float(NET[i]),1) for i in range(len(labs))}},
           "full70":{"T":int(summ['T']),"n":int(summ['n']),"p":int(summ['p']),"lambda":float(summ["lambda"]),
                     "max_eig":float(summ['max_eig']),"mean_R2":float(summ['mean_R2']),"TCI":float(summ['total_connectedness'])},
           "senders":cs.tail(4)[["name","NET"]].to_dict("records"),
           "receivers":cs.head(4)[["name","NET"]].to_dict("records")},
          open(os.path.join(DATA,"chP3_results.json"),"w",encoding="utf-8"),ensure_ascii=False,indent=2,default=float)
print("저장: ch3/figures/(fig1,fig2,fig3), ch3/data/chP3_results.json")