{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 논문 2 · 챕터 4 — 시나리오 전망\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/ch4/code/chP4_scenario.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/ch4/data/\"          # 서버의 이 챕터 데이터 폴더\n",
    "FILES = [\"dsge_exog_paths.csv\", \"exog_hist_m.csv\", \"industry_names.csv\", \"m_prod_demo6.csv\", \"q_va_demo6.csv\", \"scenario_gdp_growth.csv\", \"scenario_industry_2035.csv\", \"scenario_labor_effect.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)]\n",
    "EXO=[\"oil_g\",\"trd_g\",\"d_rrate\",\"rfx_g\",\"lab_g\"]\n",
    "SCN=[\"S1\",\"S2\",\"S3\",\"S4\",\"S5\"]; SCN_KO={\"S1\":\"기준\",\"S2\":\"제약+위기\",\"S3\":\"제약완화\",\"S4\":\"기술도약\",\"S5\":\"지정학\"}\n",
    "CLR={\"S1\":\"#333\",\"S2\":\"#c0392b\",\"S3\":\"#1f4e79\",\"S4\":\"#2e7d32\",\"S5\":\"#b8892f\"}"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 5 · 함수 정의: zscore_log, chowlin_mean, recover =====\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]; 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); 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 · 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",
    "Mlvl=np.column_stack([recover(qv[c].values, mp[c].values) for c in labs])\n",
    "Yg=100*np.diff(Mlvl,axis=0)                        # (T-1,6) 월별 성장, 2000-02~2025-12\n",
    "w=qv[labs].iloc[-1].values; w=w/w.sum()            # 2025 VA 비중(총량 가중)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 7 · 월별 외생 히스토리 정렬 =====\n",
    "# ---- 월별 외생 히스토리 정렬 ----\n",
    "ex=pd.read_csv(os.path.join(DATA,\"exog_hist_m.csv\"),encoding=\"utf-8-sig\")\n",
    "ex[\"date\"]=pd.to_datetime(ex[\"date\"]); ex=ex.set_index(\"date\").reindex(MON)\n",
    "X=ex[EXO].apply(pd.to_numeric,errors=\"coerce\").interpolate(limit_direction=\"both\").fillna(0).values\n",
    "Xg=X[1:]                                           # y와 정렬(2000-02~)\n",
    "\n",
    "# ---- VARX(p=2, 외생 시차0..1) 릿지 ----\n",
    "P,Qx,LAM=2,1,50.0"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 8 · 함수 정의: design =====\n",
    "def design(Y,Xe):\n",
    "    n=Y.shape[1]; rows=[]; tgt=[]\n",
    "    for t in range(P,len(Y)):\n",
    "        r=[Y[t-l] for l in range(1,P+1)]+[Xe[t-s] for s in range(Qx+1)]+[[1.0]]\n",
    "        rows.append(np.concatenate(r)); tgt.append(Y[t])\n",
    "    return np.array(rows),np.array(tgt)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 9 · 데이터 로드 =====\n",
    "Z,Yt=design(Yg,Xg); n=Yg.shape[1]\n",
    "D=np.ones(Z.shape[1]); D[-1]=0\n",
    "B=np.linalg.solve(Z.T@Z+LAM*np.diag(D),Z.T@Yt)     # (P*n+(Qx+1)*ne+1, n)\n",
    "\n",
    "# ---- DSGE 시나리오 외생 (분기→월, /3 반복) ----\n",
    "dsge=pd.read_csv(os.path.join(DATA,\"dsge_exog_paths.csv\"),encoding=\"utf-8-sig\")\n",
    "Hf=120"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 10 · 함수 정의: scen_exog, scen_forecast =====\n",
    "def scen_exog(s):\n",
    "    d=dsge[dsge[\"scenario\"]==s][EXO].values          # (40,5) 분기\n",
    "    return np.repeat(d/3.0,3,axis=0)[:Hf]            # (120,5) 월\n",
    "def scen_forecast(s):\n",
    "    xs=scen_exog(s); Xfull=np.vstack([Xg, xs]); start=len(Xg)\n",
    "    buf=[r.copy() for r in Yg[-P:]]; out=[]\n",
    "    for h in range(Hf):\n",
    "        ti=start+h; row=[buf[-l] for l in range(1,P+1)]+[Xfull[ti-sx] for sx in range(Qx+1)]+[[1.0]]\n",
    "        y=np.concatenate(row)@B; buf.append(y); out.append(y)\n",
    "    return np.array(out)                              # (120,6) 월 성장"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 11 · 계산 =====\n",
    "gdp_idx={}\n",
    "for s in SCN:\n",
    "    gr=scen_forecast(s); agg=(gr*w[None,:]).sum(1)     # 월 총량 성장%\n",
    "    gdp_idx[s]=100*np.exp(np.cumsum(agg/100))"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 12 · 연평균 증가율(전년대비, 연말지수) =====\n",
    "# 연평균 증가율(전년대비, 연말지수)\n",
    "def yoy(idx):\n",
    "    yl=idx[11::12]; return [100*(yl[k]/yl[k-1]-1) for k in range(1,len(yl))]"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 13 · 실제 70산업 결과 =====\n",
    "# ---- 실제 70산업 결과 ----\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",
    "g=pd.read_csv(os.path.join(DATA,\"scenario_gdp_growth.csv\"),encoding=\"utf-8-sig\")\n",
    "le=pd.read_csv(os.path.join(DATA,\"scenario_labor_effect.csv\"),encoding=\"utf-8-sig\"); le[\"name\"]=le.ind_id.map(names).fillna(le.ind_id)\n",
    "le2=le[~le.name.str.contains(\"잔여\")].sort_values(\"labor_effect_S3_S2\",ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 14 · 결과 출력 =====\n",
    "print(\"논문2·챕터4 — 시나리오 전망\")\n",
    "print(\"  [데모 6산업] 2035 총량지수(2025.12=100):\")\n",
    "for s in SCN: print(f\"    {s} {SCN_KO[s]:8} {gdp_idx[s][-1]:.1f}\")\n",
    "print(\"  [전체 70] 2026 증가율:\", g[g.year==2026].set_index(\"scenario\")[\"g_point\"].round(2).to_dict())\n",
    "print(\"  [전체 70] 노동효과(S3-S2) 상위:\")\n",
    "for _,r in le2.head(4).iterrows(): print(f\"    +{r.labor_effect_S3_S2:.1f}  {r['name']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 결과와 그림"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 15 · fig1: 실제 70산업 시나리오 증가율 팬 =====\n",
    "# ---- fig1: 실제 70산업 시나리오 증가율 팬 ----\n",
    "fig,ax=plt.subplots(figsize=(10,5))\n",
    "s1=g[g.scenario==\"S1\"].set_index(\"year\")\n",
    "ax.fill_between(s1.index, s1.g_p5, s1.g_p95, color=\"#8899aa\", alpha=.18, label=\"S1 90% 구간(몬테카를로)\")\n",
    "for s in SCN:\n",
    "    d=g[g.scenario==s].set_index(\"year\")\n",
    "    ax.plot(d.index, d.g_point, color=CLR[s], lw=2, marker=\"o\", ms=3, label=f\"{s} {SCN_KO[s]}\")\n",
    "ax.set_title(\"그림 1. 시나리오별 총 실질부가가치 전년대비 증가율 (전체 70산업)\", fontweight=\"bold\")\n",
    "ax.set_ylabel(\"전년대비 증가율 (%)\"); ax.set_xlabel(\"연도\"); ax.legend(fontsize=8.5, ncol=2); ax.grid(alpha=.25)\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig1_fan.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 16 · fig2: 실제 노동효과(S3-S2) 상위 산업 =====\n",
    "# ---- fig2: 실제 노동효과(S3-S2) 상위 산업 ----\n",
    "fig,ax=plt.subplots(figsize=(8.5,5.2)); top=le2.head(12).iloc[::-1]\n",
    "ax.barh(range(len(top)), top.labor_effect_S3_S2.values, color=\"#1f4e79\")\n",
    "ax.set_yticks(range(len(top))); ax.set_yticklabels(top.name.values, fontsize=8.5)\n",
    "ax.set_title(\"그림 2. 노동공급 효과 = 고위(S3) - 저위(S2), 2035 누적 상위 산업\", fontweight=\"bold\")\n",
    "ax.set_xlabel(\"2035 누적 부가가치 격차 (%p)\"); ax.grid(alpha=.2, axis=\"x\")\n",
    "for i,v in enumerate(top.labor_effect_S3_S2.values): ax.text(v+0.1,i,f\"+{v:.1f}\",va=\"center\",fontsize=8)\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig2_labor.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 17 · fig3: 재현 데모 — 6산업 조건부 시나리오 총량 경로 =====\n",
    "# ---- fig3: 재현 데모 — 6산업 조건부 시나리오 총량 경로 ----\n",
    "fig,ax=plt.subplots(figsize=(9.5,4.8))\n",
    "tf=pd.date_range(\"2026-01-01\",periods=Hf,freq=\"MS\")\n",
    "for s in SCN: ax.plot(tf, gdp_idx[s], color=CLR[s], lw=1.8, label=f\"{s} {SCN_KO[s]} ({gdp_idx[s][-1]:.0f})\")\n",
    "ax.set_title(\"그림 3. 재현 데모 — 6산업 VARX 조건부 시나리오 총량 경로 (2025.12=100)\", fontweight=\"bold\")\n",
    "ax.set_ylabel(\"지수 (2025.12=100)\"); ax.legend(fontsize=8.5); ax.grid(alpha=.25); ax.margins(x=.02)\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": [
    "# ===== 셀 18 · 결과 출력 =====\n",
    "json.dump({\"demo_gdp2035_idx\":{s:round(float(gdp_idx[s][-1]),1) for s in SCN},\n",
    "           \"full70_growth_2026\":{s:round(float(g[(g.scenario==s)&(g.year==2026)].g_point.iloc[0]),2) for s in SCN},\n",
    "           \"full70_growth_2035\":{s:round(float(g[(g.scenario==s)&(g.year==2035)].g_point.iloc[0]),2) for s in SCN},\n",
    "           \"labor_top\":[[r[\"name\"],round(float(r.labor_effect_S3_S2),1)] for _,r in le2.head(5).iterrows()]},\n",
    "          open(os.path.join(DATA,\"chP4_results.json\"),\"w\",encoding=\"utf-8\"),ensure_ascii=False,indent=2)\n",
    "print(\"저장: ch4/figures/(fig1,fig2,fig3), ch4/data/chP4_results.json\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "chP4_scenario.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}