{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 논문 2 · 챕터 1 — 혼합주기 상태공간(SSM)\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/ch1/code/chP1_mf_ssm.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/ch1/data/\"          # 서버의 이 챕터 데이터 폴더\n",
    "FILES = [\"latent_diagnostics.csv\", \"m_prod_demo.csv\", \"q_va_demo.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"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 5 · 계산 =====\n",
    "MON=pd.date_range(\"2000-01-01\",\"2025-12-01\",freq=\"MS\"); T=len(MON)\n",
    "QS=[f\"{2000+i//4}Q{i%4+1}\" for i in range(104)]\n",
    "QMIDX=[(i*3+2) for i in range(104)]              # 각 분기의 마지막 월 인덱스(0-based)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 6 · 함수 정의: zscore_log =====\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"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 7 · 함수 정의: chowlin_mean =====\n",
    "def chowlin_mean(q_level, qmidx, X, rho_grid=None):\n",
    "    \"\"\"월별 로그수준 m_t 복원. 제약: 분기 로그부가가치 = 3개월 m_t의 '평균'(정합 재구성오차≈0).\"\"\"\n",
    "    T,k=X.shape\n",
    "    if rho_grid is None: rho_grid=np.r_[np.linspace(0,0.95,20),0.97,0.99]\n",
    "    Q=len(q_level); A=np.zeros((Q,T))\n",
    "    for i,me in enumerate(qmidx):\n",
    "        for d in range(3):\n",
    "            j=me-d\n",
    "            if 0<=j<T: A[i,j]=1.0/3.0                 # 평균 집계\n",
    "    yq=np.asarray(q_level,float); valid=~np.isnan(yq)\n",
    "    A=A[valid]; yq=yq[valid]; Q=len(yq)\n",
    "    best=None\n",
    "    for rho in rho_grid:\n",
    "        lag=np.abs(np.subtract.outer(np.arange(T),np.arange(T)))\n",
    "        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)\n",
    "        r=yq-AX@beta; s2=float(r@Oi@r)/max(Q-k,1)\n",
    "        _,logdet=np.linalg.slogdet(Om*s2); ll=-0.5*(Q*np.log(2*np.pi)+logdet+Q)\n",
    "        if best is None or ll>best[0]:\n",
    "            u=Vu@A.T@Oi@r; best=(ll,rho,beta,X@beta+u,s2)\n",
    "    ll,rho,beta,m,s2=best\n",
    "    return dict(m=m,rho=rho,loglik=ll,s2=s2)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 8 · 함수 정의: build_X =====\n",
    "def build_X(prod):\n",
    "    cols=[np.ones(T), np.linspace(0,1,T)]\n",
    "    p=pd.Series(zscore_log(prod)).interpolate(limit_direction=\"both\").fillna(0).values\n",
    "    cols.append(p); return np.column_stack(cols)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 9 · 함수 정의: recover =====\n",
    "def recover(qva, prod, mask_last=0):\n",
    "    ylog=np.log(np.asarray(qva,float)); qy=ylog.copy(); masked=None\n",
    "    if mask_last>0: masked=qy[-mask_last:].copy(); qy[-mask_last:]=np.nan\n",
    "    res=chowlin_mean(qy, QMIDX, build_X(prod)); res[\"masked\"]=masked; res[\"qlog\"]=ylog\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 10 · 입력 =====\n",
    "# --- 입력 ---\n",
    "qv=pd.read_csv(os.path.join(DATA,\"q_va_demo.csv\"),encoding=\"utf-8-sig\")\n",
    "mp=pd.read_csv(os.path.join(DATA,\"m_prod_demo.csv\"),encoding=\"utf-8-sig\")\n",
    "IND={\"elec\":\"컴퓨터·전자·광학\",\"auto\":\"자동차·운송장비\"}\n",
    "res={}; recon={}; hold={}\n",
    "for key in [\"elec\",\"auto\"]:\n",
    "    r=recover(qv[key].values, mp[key].values); res[key]=r\n",
    "    qfit=pd.Series(r[\"m\"]).groupby(np.repeat(np.arange(104),3)).mean().values   # 분기평균\n",
    "    recon[key]=float(np.nanmax(np.abs(qfit-r[\"qlog\"])))\n",
    "    h=recover(qv[key].values, mp[key].values, mask_last=8)\n",
    "    qf2=pd.Series(h[\"m\"]).groupby(np.repeat(np.arange(104),3)).mean().values\n",
    "    pred=qf2[-8:]; act=h[\"masked\"]\n",
    "    hold[key]=dict(rmse=float(np.sqrt(np.nanmean((pred-act)**2))),\n",
    "                   corr=float(np.corrcoef(pred,act)[0,1]))\n",
    "\n",
    "diag=pd.read_csv(os.path.join(DATA,\"latent_diagnostics.csv\"),encoding=\"utf-8-sig\")"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 11 · 결과 출력 =====\n",
    "print(\"논문2·챕터1 — 혼합주기 상태공간 (재현 데모)\")\n",
    "for key in [\"elec\",\"auto\"]:\n",
    "    print(f\"  {IND[key]:14} rho={res[key]['rho']:.3f} | 재구성오차={recon[key]:.2e} | holdout corr={hold[key]['corr']:.3f}\")\n",
    "print(f\"  [70산업 전체] rho 평균={diag.rho.mean():.3f}[{diag.rho.min():.2f},{diag.rho.max():.2f}] | 재구성오차 최대={diag.max_recon_err.max():.2e}\")"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 12 · 지수화(2000.01=100) 도우미 =====\n",
    "# 지수화(2000.01=100) 도우미\n",
    "def idx100(m): return 100*np.exp(m-m[0])\n",
    "tm=MON"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 결과와 그림"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 13 · fig1: 전자 월별 복원 (2015~2025 확대) =====\n",
    "# --- fig1: 전자 월별 복원 (2015~2025 확대) ---\n",
    "fig,ax=plt.subplots(figsize=(10,4.6))\n",
    "m=res[\"elec\"][\"m\"]; lvl=idx100(m)\n",
    "proxy=zscore_log(mp[\"elec\"].values); proxy_idx=100*np.exp((proxy-np.nanmean(proxy))*0.15+ (m.mean()-m[0]))  # 배경 모양만\n",
    "# 분기 앵커: 각 분기의 월평균 로그를 지수화한 값을 분기 3개월에 계단으로\n",
    "qlog=res[\"elec\"][\"qlog\"]; qanchor=100*np.exp(qlog-m[0])\n",
    "sel=(tm>=pd.Timestamp(\"2015-01-01\"))\n",
    "ax.plot(tm[sel], lvl[sel], color=\"#185FA5\", lw=2.0, label=\"복원된 월별 부가가치\", zorder=3)\n",
    "# 분기 앵커 점(각 분기 마지막 월 위치)\n",
    "qx=[tm[QMIDX[i]] for i in range(104)]\n",
    "qsel=[i for i in range(104) if qx[i]>=pd.Timestamp(\"2015-01-01\")]\n",
    "ax.scatter([qx[i] for i in qsel],[qanchor[i] for i in qsel],color=\"#2e7d4f\",s=22,zorder=4,label=\"분기 공식치(앵커·평균)\")\n",
    "# 생산지수 프록시(표준화 모양, 우축)\n",
    "ax2=ax.twinx(); ax2.plot(tm[sel], pd.Series(mp[\"elec\"].values)[sel], color=\"#9a3b2e\", lw=1.0, ls=\"--\", alpha=.6, label=\"생산지수 프록시(우축)\")\n",
    "ax2.set_ylabel(\"생산지수\", color=\"#9a3b2e\", fontsize=9); ax2.tick_params(labelcolor=\"#9a3b2e\")\n",
    "ax.set_title(\"그림 1. 분기 부가가치를 월별로 복원 — 컴퓨터·전자·광학 (2015~2025)\", fontweight=\"bold\")\n",
    "ax.set_ylabel(\"지수 (2000.01=100)\"); ax.grid(alpha=.25)\n",
    "l1,la1=ax.get_legend_handles_labels(); l2,la2=ax2.get_legend_handles_labels()\n",
    "ax.legend(l1+l2, la1+la2, fontsize=8.5, loc=\"upper left\")\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig1_recover.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 14 · fig2: 정합성 + holdout =====\n",
    "# --- fig2: 정합성 + holdout ---\n",
    "fig,ax=plt.subplots(1,2,figsize=(12,4.4))\n",
    "# (a) full-sample 전자 월별 vs 분기평균\n",
    "m=res[\"elec\"][\"m\"]; ax[0].plot(tm, idx100(m), color=\"#185FA5\", lw=1.3, label=\"월별 복원\")\n",
    "qx=[tm[QMIDX[i]] for i in range(104)]\n",
    "ax[0].scatter(qx, 100*np.exp(res[\"elec\"][\"qlog\"]-m[0]), color=\"#2e7d4f\", s=9, label=\"분기 공식치\")\n",
    "ax[0].set_title(\"(a) 전 구간 정합 — 전자 (2000~2025)\"); ax[0].set_ylabel(\"지수(2000.01=100)\")\n",
    "ax[0].legend(fontsize=8.5); ax[0].grid(alpha=.25)\n",
    "# (b) holdout: 마지막 8분기 가리고 복원 예측 vs 실제\n",
    "for key,cl in [(\"elec\",\"#185FA5\"),(\"auto\",\"#9a3b2e\")]:\n",
    "    h=recover(qv[key].values, mp[key].values, mask_last=8)\n",
    "    qf2=pd.Series(h[\"m\"]).groupby(np.repeat(np.arange(104),3)).mean().values\n",
    "    ax[1].plot(range(8), h[\"masked\"], color=cl, lw=1.6, marker=\"o\", ms=4, label=f\"{IND[key]} 실제\")\n",
    "    ax[1].plot(range(8), qf2[-8:], color=cl, lw=1.2, ls=\"--\", marker=\"x\", ms=5, label=f\"{IND[key]} 복원(corr={hold[key]['corr']:.2f})\")\n",
    "ax[1].set_title(\"(b) 보류표본(마지막 8분기) 검증 — 로그부가가치\")\n",
    "ax[1].set_xlabel(\"보류 분기 (0=2024Q1 … 7=2025Q4)\"); ax[1].legend(fontsize=8); ax[1].grid(alpha=.25)\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig2_fit.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 15 · fig3: 70산업 rho 분포 =====\n",
    "# --- fig3: 70산업 rho 분포 ---\n",
    "fig,ax=plt.subplots(figsize=(8,3.8))\n",
    "ax.hist(diag.rho.values, bins=15, color=\"#185FA5\", edgecolor=\"white\")\n",
    "ax.axvline(diag.rho.mean(), color=\"#9a3b2e\", lw=1.5, ls=\"--\", label=f\"평균 {diag.rho.mean():.2f}\")\n",
    "ax.set_title(\"그림 3. 70개 산업의 잔차 지속성 rho 분포 (재구성오차 최대 %.0e)\"%diag.max_recon_err.max(), fontweight=\"bold\")\n",
    "ax.set_xlabel(\"rho (월별 잠재 잔차 AR(1) 계수)\"); ax.set_ylabel(\"산업 수\"); ax.legend(fontsize=9); ax.grid(alpha=.2)\n",
    "fig.tight_layout(); fig.savefig(os.path.join(FIG,\"fig3_rho.png\"),bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ===== 셀 16 · 결과 출력 =====\n",
    "json.dump({\"demo\":{k:{\"rho\":round(res[k][\"rho\"],3),\"recon_err\":recon[k],\n",
    "                      \"holdout_corr\":round(hold[k][\"corr\"],3),\"holdout_rmse\":round(hold[k][\"rmse\"],4)} for k in [\"elec\",\"auto\"]},\n",
    "           \"all70\":{\"rho_mean\":round(float(diag.rho.mean()),3),\"rho_min\":round(float(diag.rho.min()),3),\n",
    "                    \"rho_max\":round(float(diag.rho.max()),3),\"recon_max\":float(diag.max_recon_err.max()),\"n\":int(len(diag))}},\n",
    "          open(os.path.join(DATA,\"chP1_results.json\"),\"w\",encoding=\"utf-8\"),ensure_ascii=False,indent=2)\n",
    "print(\"저장: ch1/figures/(fig1,fig2,fig3), ch1/data/chP1_results.json\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "chP1_mf_ssm.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}