{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 논문 1 · 챕터 3 — BVARX (외생변수 결합)\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/sdkparkforbi/gva-scenario-2025-2035/blob/main/ch3/code/ch3_bvarx.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-scenario-2025-2035/ch3/data/\"          # 서버의 이 챕터 데이터 폴더\n",
    "FILES = [\"dlog_rates_sa.csv\", \"exog_multipliers.csv\", \"exog_quarterly_filled.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\n",
    "from matplotlib.colors import TwoSlopeNorm"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 5 · 계산 =====\n",
    "plt.rcParams[\"savefig.dpi\"]=130\n",
    "NB,NR,NG=\"#1f4e79\",\"#c0392b\",\"#2e7d32\"\n",
    "P,Q,LAM,LAMX,WIN,TESTN = 4,2,0.1,1.0,30.0,20\n",
    "EXO=[\"oil_g\",\"trd_g\",\"d_rrate\",\"rfx_g\",\"lab_g\"]\n",
    "EXOK=[\"유가\",\"교역\",\"실질금리\",\"실질환율\",\"노동\"]"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 6 · 공통표본 =====\n",
    "# --- 공통표본 ---\n",
    "rate=pd.read_csv(os.path.join(DATA,\"dlog_rates_sa.csv\"),encoding=\"utf-8-sig\",index_col=0)\n",
    "rate.index=pd.PeriodIndex(rate.index,freq=\"Q\")\n",
    "g=rate.drop(columns=[\"총량GDP\"]).dropna().clip(-WIN,WIN)\n",
    "ex=pd.read_csv(os.path.join(DATA,\"exog_quarterly_filled.csv\"),encoding=\"utf-8-sig\")\n",
    "ex.index=pd.PeriodIndex(ex[\"분기\"],freq=\"Q\"); ex=ex[EXO].astype(float)\n",
    "common=g.index.intersection(ex.index).sort_values()\n",
    "g=g.loc[common]; ex=ex.loc[common]\n",
    "names=list(g.columns); n=len(names); ne=len(EXO)\n",
    "gv=g.values; xv=ex.values; T=gv.shape[0]; test_start=T-TESTN\n",
    "sig_x_full=xv.std(0)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 7 · 함수 정의: ar_sigma =====\n",
    "def ar_sigma(m,p):\n",
    "    s=np.zeros(m.shape[1])\n",
    "    for i in range(m.shape[1]):\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))])\n",
    "        b,*_=np.linalg.lstsq(Z,Y,rcond=None); s[i]=(Y-Z@b).std(ddof=Z.shape[1])\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 8 · 함수 정의: design =====\n",
    "def design(gvv,xvv,p,q,use_x):\n",
    "    T0=max(p,q); TT=len(gvv)-T0; blocks=[]\n",
    "    for l in range(1,p+1): blocks.append(gvv[[t-l for t in range(T0,len(gvv))]])\n",
    "    if use_x:\n",
    "        for s in range(0,q+1): blocks.append(xvv[[t-s for t in range(T0,len(gvv))]])\n",
    "    X=np.column_stack(blocks+[np.ones(TT)]); Y=gvv[T0:]; return Y,X"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 9 · 함수 정의: penalty =====\n",
    "def penalty(sig_y,sig_x,p,q,use_x):\n",
    "    d=[]\n",
    "    for l in range(1,p+1):\n",
    "        for j in range(n): d.append((l**2)*(sig_y[j]**2)/LAM**2)\n",
    "    if use_x:\n",
    "        for s in range(0,q+1):\n",
    "            for m in range(ne): d.append(((s+1)**2)*(sig_x[m]**2)/LAMX**2)\n",
    "    d.append(1e-6); return np.array(d)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 10 · 함수 정의: fit =====\n",
    "def fit(gvv,xvv,p,q,use_x):\n",
    "    Y,X=design(gvv,xvv,p,q,use_x); sy=ar_sigma(gvv,p); sx=xvv.std(0)\n",
    "    d=penalty(sy,sx,p,q,use_x)\n",
    "    B=np.linalg.solve(X.T@X+np.diag(d),X.T@Y); return B"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 11 · 함수 정의: fc1 =====\n",
    "def fc1(B,gvv,xvv,tgt,p,q,use_x):    # 조건부 1단계(외생 실제값 사용)\n",
    "    x=[]\n",
    "    for l in range(1,p+1): x.append(gvv[tgt-l])\n",
    "    if use_x:\n",
    "        for s in range(0,q+1): x.append(xvv[tgt-s])\n",
    "    x=np.concatenate(x+[[1.0]]); return x@B"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 12 · 조건부 1단계 OOS: BVAR vs BVARX =====\n",
    "# --- 조건부 1단계 OOS: BVAR vs BVARX ---\n",
    "def oos(use_x):\n",
    "    err=[]\n",
    "    for k in range(test_start,T):\n",
    "        o=k-1\n",
    "        if o<max(P,Q)+2: continue\n",
    "        B=fit(gv[:k],xv[:k],P,Q,use_x)   # 학습은 target 직전까지(외생도 실제 관측)\n",
    "        err.append(gv[k]-fc1(B,gv,xv,k,P,Q,use_x))\n",
    "    E=np.array(err); return np.sqrt((E**2).mean()), np.sqrt((E**2).mean(0))\n",
    "rm_bvar=oos(False); rm_bvarx=oos(True)\n",
    "win=int((rm_bvarx[1]<rm_bvar[1]).sum())"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 13 · 전표본 BVARX 동태승수 =====\n",
    "# --- 전표본 BVARX 동태승수 ---\n",
    "B=fit(gv,xv,P,Q,True)\n",
    "A=[B[(l-1)*n:l*n,:].T for l in range(1,P+1)]\n",
    "G=[B[n*P+s*ne:n*P+(s+1)*ne,:].T for s in range(Q+1)]   # (n x ne)\n",
    "H=20; D=[]\n",
    "for h in range(H+1):\n",
    "    Dh=G[h].copy() if h<=Q else np.zeros((n,ne))\n",
    "    for l in range(1,P+1):\n",
    "        if h-l>=0: Dh=Dh+A[l-1]@D[h-l]\n",
    "    D.append(Dh)\n",
    "cum=np.sum(D,axis=0)*sig_x_full[None,:]     # 20분기 누적, 1σ  (n x ne)\n",
    "Cdf=pd.DataFrame(cum,index=names,columns=EXOK)\n",
    "Cdf.round(3).to_csv(os.path.join(DATA,\"exog_multipliers.csv\"),encoding=\"utf-8-sig\")"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 14 · 함수 정의: short =====\n",
    "def short(s): return (s.replace(\" 제조업\",\"\").replace(\" 및 \",\"·\").replace(\"서비스업\",\"서비스\"))[:9]"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 15 · 결과 출력 =====\n",
    "print(\"=\"*58); print(\"챕터3 — BVAR vs BVARX  (조건부 1단계, test 2021Q1~2025Q4)\"); print(\"=\"*58)\n",
    "print(f\"  BVAR  pooled RMSE = {rm_bvar[0]:.3f}\")\n",
    "print(f\"  BVARX pooled RMSE = {rm_bvarx[0]:.3f}   (개선 {(1-rm_bvarx[0]/rm_bvar[0])*100:+.1f}%, 우세 {win}/{n} 산업)\")\n",
    "print(\"\\n[유가 1σ↑ 20분기 누적반응 — 상위/하위]\")\n",
    "oil=Cdf[\"유가\"].sort_values()\n",
    "for nm,v in list(oil.items())[-3:][::-1]: print(f\"  수혜 {v:+.2f}  {short(nm)}\")\n",
    "for nm,v in list(oil.items())[:3]: print(f\"  피해 {v:+.2f}  {short(nm)}\")"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 16 · 표본내 적합 R² (H3 검증) =====\n",
    "# --- 표본내 적합 R² (H3 검증) ---\n",
    "def insample_r2(use_x):\n",
    "    Y,X=design(gv,xv,P,Q,use_x); Bf=fit(gv,xv,P,Q,use_x); r=Y-X@Bf\n",
    "    return float(np.mean(1-r.var(0)/Y.var(0)))\n",
    "r2b=insample_r2(False); r2x=insample_r2(True)\n",
    "print(f\"\\n표본내 R²: BVAR {r2b:.3f} -> BVARX {r2x:.3f}  (외생블록으로 +{(r2x-r2b)*100:.1f}%p, 논문 H3의 표본내 개선)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 결과와 그림"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 17 · fig1: 표본내 R²(상승) vs 표본외 RMSE(불변) =====\n",
    "# --- fig1: 표본내 R²(상승) vs 표본외 RMSE(불변) ---\n",
    "fig,ax=plt.subplots(1,2,figsize=(11,4.4))\n",
    "fig.suptitle(\"그림 1. 외생블록 — 표본내 적합은 오르나 표본외 예측은 개선 없음\",fontweight=\"bold\")\n",
    "b0=ax[0].bar([\"BVAR\",\"BVARX\"],[r2b,r2x],color=[NR,NB]); ax[0].bar_label(b0,fmt=\"%.3f\")\n",
    "ax[0].set_title(\"(a) 표본내 R² — 외생블록으로 상승 (H3)\"); ax[0].set_ylabel(\"평균 R²\")\n",
    "b1=ax[1].bar([\"BVAR\",\"BVARX\"],[rm_bvar[0],rm_bvarx[0]],color=[NR,NB]); ax[1].bar_label(b1,fmt=\"%.3f\")\n",
    "ax[1].set_title(f\"(b) 표본외 1단계 RMSE — 대등/불변 (우세 {win}/{n})\"); ax[1].set_ylabel(\"RMSE (%p)\")\n",
    "lo=min(rm_bvar[0],rm_bvarx[0]); ax[1].set_ylim(lo*0.97,max(rm_bvar[0],rm_bvarx[0])*1.02)\n",
    "fig.tight_layout(rect=[0,0,1,.93]); fig.savefig(os.path.join(FIG,\"fig1_bvar_vs_bvarx.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 18 · fig2: 동태승수 히트맵 (산업 x 외생충격) =====\n",
    "# --- fig2: 동태승수 히트맵 (산업 x 외생충격) ---\n",
    "fig,ax=plt.subplots(figsize=(7,11))\n",
    "v=np.percentile(np.abs(cum),95)\n",
    "im=ax.imshow(cum,cmap=\"RdBu_r\",norm=TwoSlopeNorm(0,-v,v),aspect=\"auto\")\n",
    "ax.set_xticks(range(ne)); ax.set_xticklabels(EXOK,fontsize=9)\n",
    "ax.set_yticks(range(n)); ax.set_yticklabels([short(s) for s in names],fontsize=7)\n",
    "ax.set_title(\"그림 2. 외생충격 1σ의 산업별 20분기 누적반응\\n(붉음=수혜, 청=피해)\",fontsize=11)\n",
    "fig.colorbar(im,ax=ax,shrink=.7,pad=.02,label=\"%p\")\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": [
    "# ===== 셀 19 · fig3: 유가충격 산업별 반응(비대칭) =====\n",
    "# --- fig3: 유가충격 산업별 반응(비대칭) ---\n",
    "fig,ax=plt.subplots(figsize=(8,9))\n",
    "oil_s=Cdf[\"유가\"].sort_values()\n",
    "cols=[NR if x<0 else NG for x in oil_s.values]\n",
    "ax.barh(range(n),oil_s.values,color=cols)\n",
    "ax.set_yticks(range(n)); ax.set_yticklabels([short(s) for s in oil_s.index],fontsize=7)\n",
    "ax.axvline(0,color=\"#333\",lw=.8)\n",
    "ax.set_title(\"그림 3. 유가 1σ↑ 충격의 산업별 반응 (비대칭)\\n초록=수혜, 빨강=피해\")\n",
    "ax.set_xlabel(\"20분기 누적 %p\")\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig3_oil_response.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 20 · 결과 출력 =====\n",
    "json.dump({\"rmse\":{\"BVAR\":round(float(rm_bvar[0]),3),\"BVARX\":round(float(rm_bvarx[0]),3)},\n",
    "           \"r2\":{\"BVAR\":round(r2b,3),\"BVARX\":round(r2x,3)},\n",
    "           \"improve_pct\":round((1-rm_bvarx[0]/rm_bvar[0])*100,1),\"win\":f\"{win}/{n}\",\n",
    "           \"oil_top\":[[short(k),round(float(v),2)] for k,v in oil.tail(3).items()],\n",
    "           \"oil_bottom\":[[short(k),round(float(v),2)] for k,v in oil.head(3).items()]},\n",
    "          open(os.path.join(DATA,\"ch3_results.json\"),\"w\",encoding=\"utf-8\"),ensure_ascii=False,indent=2)\n",
    "print(\"\\n저장: ch3/figures/(fig1,fig2,fig3), ch3/data/(exog_multipliers.csv, ch3_results.json)\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "ch3_bvarx.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}