{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 논문 2 · 챕터 3 — VARX · 연결성\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/sdkparkforbi/gva-monthly-ssm-70/blob/main/ch3/code/chP3_varx.ipynb)\n",
    "\n",
    "> **Google Colab에서 바로 실행**됩니다. 배지를 누른 뒤 `런타임 → 모두 실행`.\n",
    "> 데이터는 아래 셀에서 자동으로 내려받고, 그림은 셀 실행 시 바로 표시됩니다.\n",
    "\n",
    "각 셀에는 **셀 번호 · 셀 제목**이 달려 있고, 코드에는 줄별 **주석**을 넣었습니다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. 준비 — Colab 환경과 데이터"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 1 · 한글 폰트 (koreanize-matplotlib) =====\n",
    "# 그림의 한글이 깨지지 않도록 koreanize-matplotlib 를 설치하고 불러옵니다.\n",
    "# import 하는 순간 matplotlib 한글 폰트가 자동 설정됩니다.\n",
    "_font_ok = False\n",
    "try:\n",
    "    import koreanize_matplotlib                       # 이미 설치돼 있으면 바로 사용\n",
    "    _font_ok = True\n",
    "except Exception:\n",
    "    try:\n",
    "        import subprocess, sys\n",
    "        subprocess.run([sys.executable,'-m','pip','install','-q','koreanize-matplotlib'], check=False)\n",
    "        import koreanize_matplotlib                   # 설치 후 다시 불러오기\n",
    "        _font_ok = True\n",
    "    except Exception:\n",
    "        _font_ok = False\n",
    "if not _font_ok:                                      # 폴백: 나눔폰트 직접 설치·등록\n",
    "    import subprocess\n",
    "    subprocess.run(['apt-get','-qq','-y','install','fonts-nanum'], check=False)\n",
    "    import matplotlib.font_manager as fm, matplotlib.pyplot as plt\n",
    "    _p = '/usr/share/fonts/truetype/nanum/NanumGothic.ttf'\n",
    "    try:\n",
    "        fm.fontManager.addfont(_p)\n",
    "        plt.rcParams['font.family'] = fm.FontProperties(fname=_p).get_name()  # 실제 등록 이름 사용\n",
    "    except Exception: pass\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['axes.unicode_minus'] = False           # 마이너스 기호 깨짐 방지\n",
    "print('한글 폰트 준비 완료')"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 2 · 데이터 다운로드 =====\n",
    "# 이 챕터의 입력 데이터를 서버에서 로컬 data/ 폴더로 내려받습니다.\n",
    "import os, urllib.request, urllib.parse\n",
    "BASE = \"https://aiforalab.com/gva-monthly-ssm-70/ch3/data/\"          # 서버의 이 챕터 데이터 폴더\n",
    "FILES = [\"connectedness.csv\", \"industry_names.csv\", \"m_prod_demo6.csv\", \"q_va_demo6.csv\", \"varx_exog_multipliers.csv\", \"varx_summary.csv\"]\n",
    "os.makedirs(\"data\", exist_ok=True)\n",
    "for f in FILES:                                      # 한글 파일명은 URL 인코딩 후 요청\n",
    "    urllib.request.urlretrieve(BASE + urllib.parse.quote(f), os.path.join('data', f))\n",
    "print('다운로드 완료:', len(FILES), '개 ->', sorted(os.listdir('data')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 방법 구현과 계산"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 3 · 라이브러리 임포트·설정 =====\n",
    "import os, json\n",
    "import numpy as np, pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 4 · Colab 경로: 데이터=data/, 그림=figures/ =====\n",
    "# ----- Colab 경로: 데이터=data/, 그림=figures/ -----\n",
    "import os\n",
    "DATA = \"data\"; PREP = \"data\"; FIG = \"figures\"; OUT = \"figures\"\n",
    "os.makedirs(\"data\", exist_ok=True); os.makedirs(\"figures\", exist_ok=True)\n",
    "plt.rcParams[\"savefig.dpi\"]=130\n",
    "MON=pd.date_range(\"2000-01-01\",\"2025-12-01\",freq=\"MS\"); T=len(MON)\n",
    "QMIDX=[(i*3+2) for i in range(104)]"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 5 · SSM (챕터1과 동일, 평균집계 Chow-Lin) =====\n",
    "# ---------- SSM (챕터1과 동일, 평균집계 Chow-Lin) ----------\n",
    "def zscore_log(v):\n",
    "    v=np.asarray(v,float); out=np.full_like(v,np.nan); m=v>0\n",
    "    lg=np.log(v[m]); out[m]=(lg-lg.mean())/(lg.std()+1e-9); return out\n",
    "def chowlin_mean(q_level, X):\n",
    "    Tn,k=X.shape; rho_grid=np.r_[np.linspace(0,0.95,20),0.97,0.99]\n",
    "    Q=len(q_level); A=np.zeros((Q,Tn))\n",
    "    for i,me in enumerate(QMIDX):\n",
    "        for d in range(3):\n",
    "            j=me-d\n",
    "            if 0<=j<Tn: A[i,j]=1/3\n",
    "    yq=np.asarray(q_level,float); vv=~np.isnan(yq); A=A[vv]; yq=yq[vv]; Q=len(yq)\n",
    "    best=None\n",
    "    for rho in rho_grid:\n",
    "        lag=np.abs(np.subtract.outer(np.arange(Tn),np.arange(Tn))); Vu=rho**lag/(1-rho**2+1e-12)\n",
    "        AX=A@X; Om=A@Vu@A.T+1e-8*np.eye(Q); Oi=np.linalg.pinv(Om)\n",
    "        beta=np.linalg.solve(AX.T@Oi@AX+1e-8*np.eye(k),AX.T@Oi@yq); r=yq-AX@beta\n",
    "        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)\n",
    "        if best is None or ll>best[0]: best=(ll,rho,X@beta+Vu@A.T@Oi@r)\n",
    "    return best[2]\n",
    "def recover(qva, proxy):\n",
    "    X=np.column_stack([np.ones(T),np.linspace(0,1,T),\n",
    "        pd.Series(zscore_log(proxy)).interpolate(limit_direction=\"both\").bfill().ffill().values])\n",
    "    return chowlin_mean(np.log(np.asarray(qva,float)), X)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 6 · VARX(p) 릿지 + 동반행렬 고유값 + GFEVD 연결성 =====\n",
    "# ---------- VARX(p) 릿지 + 동반행렬 고유값 + GFEVD 연결성 ----------\n",
    "def var_ridge(Y, p=2, lam=1.0):\n",
    "    Tn,n=Y.shape; Z=[];\n",
    "    for t in range(p,Tn): Z.append(np.concatenate([Y[t-l] for l in range(1,p+1)]+[[1.0]]))\n",
    "    Z=np.array(Z); Yt=Y[p:]\n",
    "    D=np.ones(Z.shape[1]); D[-1]=0.0                      # 절편 비축소\n",
    "    B=np.linalg.solve(Z.T@Z+lam*np.diag(D), Z.T@Yt)       # (np+1, n)\n",
    "    resid=Yt-Z@B; Sig=np.cov(resid.T)\n",
    "    Blags=[B[l*n:(l+1)*n].T for l in range(p)]            # 각 (n,n)\n",
    "    return Blags, Sig, B\n",
    "def companion_eig(Blags):\n",
    "    n=Blags[0].shape[0]; p=len(Blags)\n",
    "    C=np.zeros((n*p,n*p)); C[:n]=np.hstack(Blags)\n",
    "    if p>1: C[n:]=np.eye(n*(p-1),n*p)\n",
    "    return np.max(np.abs(np.linalg.eigvals(C)))\n",
    "def gfevd_connectedness(Blags, Sig, H=12):\n",
    "    n=Blags[0].shape[0]; p=len(Blags)\n",
    "    Psi=[np.eye(n)]\n",
    "    for h in range(1,H+1):\n",
    "        s=np.zeros((n,n))\n",
    "        for l in range(1,p+1):\n",
    "            if h-l>=0: s+=Blags[l-1]@Psi[h-l]\n",
    "        Psi.append(s)\n",
    "    sig=np.diag(Sig); theta=np.zeros((n,n))\n",
    "    for i in range(n):\n",
    "        denom=sum((Psi[h][i,:]@Sig@Psi[h][i,:].T) for h in range(H+1))\n",
    "        for j in range(n):\n",
    "            num=sum((Psi[h][i,:]@Sig[:,j])**2 for h in range(H+1))/sig[j]\n",
    "            theta[i,j]=num/denom\n",
    "    theta=theta/theta.sum(1,keepdims=True)                # 행 정규화\n",
    "    FROM=theta.sum(1)-np.diag(theta); TO=theta.sum(0)-np.diag(theta)\n",
    "    NET=TO-FROM; TCI=100*FROM.mean()\n",
    "    return theta*100, FROM*100, TO*100, NET*100, TCI"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 7 · 재현 데모: 6 제조업 =====\n",
    "# ---------- 재현 데모: 6 제조업 ----------\n",
    "qv=pd.read_csv(os.path.join(DATA,\"q_va_demo6.csv\"),encoding=\"utf-8-sig\")\n",
    "mp=pd.read_csv(os.path.join(DATA,\"m_prod_demo6.csv\"),encoding=\"utf-8-sig\")\n",
    "labs=[c for c in qv.columns if c!=\"quarter\"]\n",
    "M=np.column_stack([recover(qv[c].values, mp[c].values) for c in labs])   # (T,6) 로그수준\n",
    "Yg=100*np.diff(M,axis=0)                                                  # 월별 성장률\n",
    "Bl,Sig,_=var_ridge(Yg,p=2,lam=1.0)\n",
    "eig_ridge=companion_eig(Bl)\n",
    "Bl_ols,Sig_ols,_=var_ridge(Yg,p=2,lam=1e-9)\n",
    "eig_ols=companion_eig(Bl_ols)\n",
    "theta,FROM,TO,NET,TCI=gfevd_connectedness(Bl,Sig,H=12)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 8 · 결과 출력 =====\n",
    "print(\"논문2·챕터3 — VARX & 연결성 (재현 데모 6산업)\")\n",
    "print(f\"  동반행렬 max|eig|: 릿지={eig_ridge:.3f}  OLS={eig_ols:.3f}  (안정=1 미만)\")\n",
    "print(f\"  데모 총연결성 TCI={TCI:.1f}%\")\n",
    "for i,l in enumerate(sorted(range(len(labs)),key=lambda k:-NET[k])):\n",
    "    print(f\"    NET {NET[l]:+6.1f}  {labs[l]}\")"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 9 · 전체 70산업 실제 결과 =====\n",
    "# 전체 70산업 실제 결과\n",
    "summ=pd.read_csv(os.path.join(DATA,\"varx_summary.csv\"),encoding=\"utf-8-sig\").iloc[0]\n",
    "conn=pd.read_csv(os.path.join(DATA,\"connectedness.csv\"),encoding=\"utf-8-sig\")\n",
    "names=pd.read_csv(os.path.join(DATA,\"industry_names.csv\"),encoding=\"utf-8-sig\").set_index(\"ind_id\")[\"ind_name\"].to_dict()\n",
    "conn[\"name\"]=conn.ind_id.map(names).fillna(conn.ind_id)\n",
    "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']}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 결과와 그림"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 10 · fig1: 70산업 순연결성(NET) 상·하위 =====\n",
    "# ---------- fig1: 70산업 순연결성(NET) 상·하위 ----------\n",
    "cs=conn.sort_values(\"NET\"); top=cs.tail(8); bot=cs.head(6); rows=pd.concat([bot,top])\n",
    "fig,ax=plt.subplots(figsize=(8.5,6.5))\n",
    "ax.barh(range(len(rows)),rows.NET.values,color=[\"#9a3b2e\" if v<0 else \"#185FA5\" for v in rows.NET])\n",
    "ax.axvline(0,color=\"#333\",lw=.8); ax.set_yticks(range(len(rows))); ax.set_yticklabels(rows.name.values,fontsize=8)\n",
    "ax.set_title(f\"그림 1. 70산업 순 연결성(NET) 상·하위 — 총연결성 TCI {summ['total_connectedness']}%\", fontweight=\"bold\")\n",
    "ax.set_xlabel(\"순 연결성 NET (%p) : 양(+)=발신 / 음(-)=수신\")\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig1_net.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 11 · fig2: 외생충격 동태승수 히트맵 (대표 산업) =====\n",
    "# ---------- fig2: 외생충격 동태승수 히트맵 (대표 산업) ----------\n",
    "em=pd.read_csv(os.path.join(DATA,\"varx_exog_multipliers.csv\"),encoding=\"utf-8-sig\").rename(columns={\"Unnamed: 0\":\"ind_id\"})\n",
    "em[\"name\"]=em.ind_id.map(names).fillna(em.ind_id)\n",
    "pick=[\"MFG_COKE\",\"MFG_CHEM\",\"E_SEMI\",\"MFG_AUTO\",\"MFG_METAL\",\"MFG_MACH\",\"MFG_TEX\",\"CON_CIVIL\",\"WHOLE\",\"FIN\",\"TRANS\",\"AGR\"]\n",
    "sub=em[em.ind_id.isin(pick)].set_index(\"name\")[[\"oil_g\",\"trd_g\",\"d_rrate\",\"rfx_g\",\"gpr_l\",\"lab_g\"]]\n",
    "sub=sub.reindex([names.get(p,p) for p in pick if names.get(p,p) in sub.index])\n",
    "XL={\"oil_g\":\"유가\",\"trd_g\":\"교역\",\"d_rrate\":\"실질금리\",\"rfx_g\":\"실질환율\",\"gpr_l\":\"지정학\",\"lab_g\":\"고용\"}\n",
    "fig,ax=plt.subplots(figsize=(7.5,6.2))\n",
    "V=np.clip(sub.values,-1,1)\n",
    "im=ax.imshow(V,cmap=\"RdBu_r\",vmin=-1,vmax=1,aspect=\"auto\")\n",
    "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)\n",
    "for i in range(len(sub)):\n",
    "    for j in range(6): ax.text(j,i,f\"{sub.values[i,j]:+.2f}\",ha=\"center\",va=\"center\",fontsize=7,color=\"#222\")\n",
    "ax.set_title(\"그림 2. 외생충격 누적 동태승수 (대표 산업 × 6 외생)\", fontweight=\"bold\")\n",
    "fig.colorbar(im,ax=ax,fraction=.046,pad=.04,label=\"누적 반응 (색 clip ±1)\")\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig2_multipliers.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 12 · fig3: 재현 데모 6x6 연결성 히트맵 =====\n",
    "# ---------- fig3: 재현 데모 6x6 연결성 히트맵 ----------\n",
    "fig,ax=plt.subplots(figsize=(7,5.8))\n",
    "im=ax.imshow(theta,cmap=\"Blues\",aspect=\"auto\")\n",
    "ax.set_xticks(range(len(labs))); ax.set_xticklabels(labs,rotation=35,ha=\"right\",fontsize=8.5)\n",
    "ax.set_yticks(range(len(labs))); ax.set_yticklabels(labs,fontsize=8.5)\n",
    "for i in range(len(labs)):\n",
    "    for j in range(len(labs)):\n",
    "        ax.text(j,i,f\"{theta[i,j]:.0f}\",ha=\"center\",va=\"center\",fontsize=7.5,\n",
    "                color=\"white\" if theta[i,j]>theta.max()*0.55 else \"#222\")\n",
    "ax.set_xlabel(\"← 로부터(발신 산업)\"); ax.set_ylabel(\"수신 산업\")\n",
    "ax.set_title(f\"그림 3. 재현 데모 — 6산업 연결성 행렬(%) · TCI {TCI:.0f}% · max|eig| {eig_ridge:.2f}\", fontweight=\"bold\", fontsize=11)\n",
    "fig.colorbar(im,ax=ax,fraction=.046,pad=.04,label=\"분산기여 (%)\")\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig3_demo.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 13 · 결과 출력 =====\n",
    "json.dump({\"demo\":{\"tci\":round(TCI,1),\"eig_ridge\":round(eig_ridge,3),\"eig_ols\":round(eig_ols,3),\n",
    "                   \"net\":{labs[i]:round(float(NET[i]),1) for i in range(len(labs))}},\n",
    "           \"full70\":{\"T\":int(summ['T']),\"n\":int(summ['n']),\"p\":int(summ['p']),\"lambda\":float(summ[\"lambda\"]),\n",
    "                     \"max_eig\":float(summ['max_eig']),\"mean_R2\":float(summ['mean_R2']),\"TCI\":float(summ['total_connectedness'])},\n",
    "           \"senders\":cs.tail(4)[[\"name\",\"NET\"]].to_dict(\"records\"),\n",
    "           \"receivers\":cs.head(4)[[\"name\",\"NET\"]].to_dict(\"records\")},\n",
    "          open(os.path.join(DATA,\"chP3_results.json\"),\"w\",encoding=\"utf-8\"),ensure_ascii=False,indent=2,default=float)\n",
    "print(\"저장: ch3/figures/(fig1,fig2,fig3), ch3/data/chP3_results.json\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "chP3_varx.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}