{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 논문 1 · 챕터 2 — 미네소타 vs I-O 프라이어 비교\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/ch2/code/ch2_bvar_compare.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/ch2/data/\"          # 서버의 이 챕터 데이터 폴더\n",
    "FILES = [\"dlog_rates_sa.csv\", \"io_linkage_W.csv\", \"io_mapping.csv\", \"io_prior_omega.csv\", \"io_production_inducement_long.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",
    "P, LAM, WIN, TESTN = 4, 0.1, 30.0, 20\n",
    "NB, NR, NG = \"#1f4e79\", \"#c0392b\", \"#2e7d32\""
   ]
  },
  {
   "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",
    "names=list(g.columns); n=len(names); gv=g.values; T=gv.shape[0]\n",
    "test_start=T-TESTN\n",
    "Om=pd.read_csv(os.path.join(DATA,\"io_prior_omega.csv\"),encoding=\"utf-8-sig\",index_col=0).values\n",
    "off=~np.eye(n,dtype=bool)\n",
    "Om[off]=Om[off]/Om[off].mean()          # 오프대각 평균 1로 재정규화(공정)\n",
    "np.fill_diagonal(Om,1.0)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 7 · 함수 정의: build_XY, ar_sigma =====\n",
    "def build_XY(m,p):\n",
    "    Y=m[p:]; TT=Y.shape[0]; X=np.ones((TT,n*p+1))\n",
    "    for l in range(1,p+1): X[:,(l-1)*n:l*n]=m[p-l:-l]\n",
    "    return Y,X\n",
    "def ar_sigma(m,p):\n",
    "    s=np.zeros(n)\n",
    "    for i in range(n):\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 · 함수 정의: penalty_base =====\n",
    "def penalty_base(sig):\n",
    "    \"\"\"릿지 벌점 (l,j) = l²σ_j²/λ²  (미네소타, 방정식 공통). const≈0.\"\"\"\n",
    "    d=np.zeros(n*P+1)\n",
    "    for l in range(1,P+1):\n",
    "        for j in range(n): d[(l-1)*n+j]=(l**2)*(sig[j]**2)/(LAM**2)\n",
    "    d[-1]=1e-6\n",
    "    return d"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 9 · 함수 정의: fit =====\n",
    "def fit(m,p,mode,sig):\n",
    "    Y,X=build_XY(m,p); XtX=X.T@X; XtY=X.T@Y\n",
    "    d0=penalty_base(sig)\n",
    "    if mode==\"MIN\":\n",
    "        B=np.linalg.solve(XtX+np.diag(d0),XtY)         # 방정식 공통\n",
    "    else:  # I-O: 방정식별 (교차항 벌점 /= ω_ij)\n",
    "        B=np.zeros((n*P+1,n))\n",
    "        for i in range(n):\n",
    "            d=d0.copy()\n",
    "            for l in range(1,p+1):\n",
    "                for j in range(n):\n",
    "                    if j!=i: d[(l-1)*n+j]/=Om[i,j]      # 연계 강하면(ω>1) 벌점↓\n",
    "            B[:,i]=np.linalg.solve(XtX+np.diag(d),XtY[:,i])\n",
    "    return B"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 10 · 함수 정의: forecast =====\n",
    "def forecast(B,m,p,h):\n",
    "    buf=[r.copy() for r in m[-p:]]\n",
    "    for _ in range(h):\n",
    "        x=np.ones(n*p+1)\n",
    "        for l in range(1,p+1): x[(l-1)*n:l*n]=buf[-l]\n",
    "        buf.append(x@B)\n",
    "    return buf[-1]"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 11 · 표본외 검증 (확장창) =====\n",
    "# --- 표본외 검증 (확장창) ---\n",
    "def evaluate(mode):\n",
    "    res={1:[],4:[]}\n",
    "    for k in range(test_start,T):\n",
    "        for h in (1,4):\n",
    "            o=k-h\n",
    "            if o<P+4: continue\n",
    "            tr=gv[:o+1]; sig=ar_sigma(tr,P); B=fit(tr,P,mode,sig)\n",
    "            res[h].append((gv[k], forecast(B,tr,P,h)))\n",
    "    return res\n",
    "EV={m:evaluate(m) for m in [\"MIN\",\"IO\"]}"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 12 · 함수 정의: rmse =====\n",
    "def rmse(mode,h):\n",
    "    a=np.array([x[0] for x in EV[mode][h]]); f=np.array([x[1] for x in EV[mode][h]])\n",
    "    return np.sqrt(((a-f)**2).mean()), np.sqrt(((a-f)**2).mean(0))  # pooled, per-industry\n",
    "RM={m:{h:rmse(m,h) for h in (1,4)} for m in [\"MIN\",\"IO\"]}"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 13 · 전표본 계수 구조 (fig3) =====\n",
    "# --- 전표본 계수 구조 (fig3) ---\n",
    "sig_full=ar_sigma(gv,P)\n",
    "Bmin=fit(gv,P,\"MIN\",sig_full); Bio=fit(gv,P,\"IO\",sig_full)\n",
    "A1min=Bmin[0:n,:].T; A1io=Bio[0:n,:].T           # [반응 i, 충격 j]"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 14 · 콘솔 =====\n",
    "# --- 콘솔 ---\n",
    "print(\"=\"*60)\n",
    "print(\"챕터2 — 미네소타 vs I-O prior  (동일 test, 방정식별 릿지)\")\n",
    "print(\"=\"*60)\n",
    "print(f\"{'':10}{'1단계 pooled':>14}{'1년 pooled':>12}\")\n",
    "for m,lb in [(\"MIN\",\"미네소타\"),(\"IO\",\"I-O prior\")]:\n",
    "    print(f\"{lb:10}{RM[m][1][0]:>14.3f}{RM[m][4][0]:>12.3f}\")\n",
    "win1=int((RM['IO'][1][1]<RM['MIN'][1][1]).sum()); win4=int((RM['IO'][4][1]<RM['MIN'][4][1]).sum())\n",
    "print(f\"\\nI-O가 미네소타보다 나은 산업: 1단계 {win1}/{n}, 1년 {win4}/{n}\")\n",
    "print(f\"교차 |A1| 평균: 미네소타 {np.abs(A1min[off]).mean():.4f}  I-O {np.abs(A1io[off]).mean():.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 결과와 그림"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 15 · fig2: test RMSE 비교 (표와 동일수치) =====\n",
    "# --- fig2: test RMSE 비교 (표와 동일수치) ---\n",
    "fig,ax=plt.subplots(1,2,figsize=(13,4.6))\n",
    "fig.suptitle(\"그림 2. 미네소타 vs I-O prior — 동일 test 구간 예측 RMSE\",fontweight=\"bold\")\n",
    "x=np.arange(2); w=.38\n",
    "for hi,h in enumerate((1,4)):\n",
    "    vals=[RM[\"MIN\"][h][0],RM[\"IO\"][h][0]]\n",
    "    b=ax[hi].bar([\"미네소타\",\"I-O prior\"],vals,color=[NR,NB])\n",
    "    ax[hi].bar_label(b,fmt=\"%.3f\")\n",
    "    ax[hi].set_title(f\"({'a' if h==1 else 'b'}) {'1단계' if h==1 else '1년앞(4Q)'} pooled RMSE\")\n",
    "    ax[hi].set_ylabel(\"RMSE (%p)\")\n",
    "    lo=min(vals); ax[hi].set_ylim(lo*0.985,max(vals)*1.01)\n",
    "fig.tight_layout(rect=[0,0,1,.93]); fig.savefig(os.path.join(FIG,\"fig2_rmse_compare.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 16 · fig3: 계수 구조 (미네소타 vs I-O |A1|) =====\n",
    "# --- fig3: 계수 구조 (미네소타 vs I-O |A1|) ---\n",
    "def short(s): return (s.replace(\" 제조업\",\"\").replace(\" 및 \",\"·\").replace(\"서비스업\",\"서비스\"))[:9]"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 17 · 그림 그리기 =====\n",
    "lab=[short(s) for s in names]\n",
    "v=np.percentile(np.abs(np.r_[A1min[off],A1io[off]]),98)\n",
    "fig,ax=plt.subplots(1,2,figsize=(15,7.2))\n",
    "fig.suptitle(\"그림 3. 1차 계수 |A₁| 구조 — I-O prior는 연계 약한 쌍을 더 0으로 (열 j→행 i)\",fontweight=\"bold\")\n",
    "for a_,Mx,tt in [(ax[0],np.abs(A1min),\"(a) 미네소타\"),(ax[1],np.abs(A1io),\"(b) I-O prior\")]:\n",
    "    im=a_.imshow(Mx,cmap=\"YlOrRd\",vmin=0,vmax=v,aspect=\"auto\"); a_.set_title(tt)\n",
    "    a_.set_xticks(range(n)); a_.set_xticklabels(lab,rotation=90,fontsize=5)\n",
    "    a_.set_yticks(range(n)); a_.set_yticklabels(lab,fontsize=5)\n",
    "    fig.colorbar(im,ax=a_,shrink=.8,pad=.02)\n",
    "fig.tight_layout(rect=[0,0,1,.95]); fig.savefig(os.path.join(FIG,\"fig3_coef_structure.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 18 · 결과 저장 =====\n",
    "# --- 결과 저장 ---\n",
    "json.dump({\"rmse\":{m:{str(h):round(float(RM[m][h][0]),3) for h in (1,4)} for m in [\"MIN\",\"IO\"]},\n",
    "           \"win_io\":{\"h1\":f\"{win1}/{n}\",\"h4\":f\"{win4}/{n}\"},\n",
    "           \"absA1_cross\":{\"MIN\":round(float(np.abs(A1min[off]).mean()),4),\"IO\":round(float(np.abs(A1io[off]).mean()),4)}},\n",
    "          open(os.path.join(DATA,\"ch2_results.json\"),\"w\",encoding=\"utf-8\"),ensure_ascii=False,indent=2)\n",
    "print(\"\\n저장: ch2/figures/(fig2,fig3), ch2/data/ch2_results.json\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "ch2_bvar_compare.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}