{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# nastab_elasticity.py\n",
    "\n",
    "재정패널(NaSTaB) 가구 원자료 — 소득분위별 가격탄력성 추정\n",
    "  · 목적: H2 가격탄력성 + 동태 세수 시뮬레이션 파라미터 확정\n",
    "  · 식별: 2024 소주 기준판매비율(가격 인하) 자연실험 + 주종 CPI\n",
    "  · 핵심 통제: 시간에 따른 '주류소비문화' 변화(저도주화·혼술·건강·고령화)를\n",
    "             유연한 추세 + 코호트/연령/가구특성으로 흡수 → 가격효과 정화\n",
    "  · 실행 환경: KIPF 승인환경(가구 원자료 wave 1~18, 조사 2008~2025)\n",
    "  ─ 아래 변수코드/파일경로는 배포본에 맞게 수정하세요(주석 표시).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "재정패널(NaSTaB) 가구 원자료 — 소득분위별 가격탄력성 추정\n",
    "  · 목적: H2 가격탄력성 + 동태 세수 시뮬레이션 파라미터 확정\n",
    "  · 식별: 2024 소주 기준판매비율(가격 인하) 자연실험 + 주종 CPI\n",
    "  · 핵심 통제: 시간에 따른 '주류소비문화' 변화(저도주화·혼술·건강·고령화)를\n",
    "             유연한 추세 + 코호트/연령/가구특성으로 흡수 → 가격효과 정화\n",
    "  · 실행 환경: KIPF 승인환경(가구 원자료 wave 1~18, 조사 2008~2025)\n",
    "  ─ 아래 변수코드/파일경로는 배포본에 맞게 수정하세요(주석 표시).\n",
    "\"\"\"\n",
    "import numpy as np, pandas as pd\n",
    "# linearmodels 없으면: pip 불가 환경이면 statsmodels OLS+더미로 대체(하단 주석)\n",
    "from linearmodels.panel import PanelOLS\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# ── 0. 설정 ────────────────────────────────────────────────\n",
    "WAVES = range(1, 19)                 # 1~18차 (조사 2008~2025, 기준 2007~2024)\n",
    "RAW = lambda w: f\"nastab_{w:02d}.dta\"   # ← KIPF 배포 파일명 패턴에 맞게 수정\n",
    "CPI_FILE = \"cpi_liquor_annual.csv\"      # 연도(ref_year), 소주, 맥주, 주류종합 CPI\n",
    "POLICY_YEAR = 2024                       # 기준판매비율 시행(소주 가격 인하)\n",
    "\n",
    "# 변수코드(재정패널) ── findings 기준, 배포본에 맞게 확인/수정\n",
    "V = dict(alc_home1=\"cd002\",   # 가정용 주류비(맥주·소주·양주)\n",
    "         alc_home2=\"cd004\",   # 가정용 주류비(와인·막걸리·전통주 등)\n",
    "         cons=\"ca002\",        # 소비지출\n",
    "         inc=\"ca004\",         # 경상소득(개방형, 4차~)\n",
    "         age=\"자녀가구주_연령\",  # ← 가구주 연령 코드로 수정\n",
    "         hhsize=\"가구원수\",      # ← 가구원수 코드로 수정\n",
    "         weight_tpl=\"h{w}cwt\")   # 횡단가중치\n",
    "\n",
    "# ── 1. 가구 패널 적재 ──────────────────────────────────────\n",
    "frames=[]\n",
    "for w in WAVES:\n",
    "    try:\n",
    "        d = pd.read_stata(RAW(w), convert_categoricals=False)\n",
    "    except Exception as e:\n",
    "        print(f\"[wave {w}] 적재 실패: {e}\"); continue\n",
    "    ref_year = 2006 + w                      # wave w → 기준연도(ref) = 2006+w\n",
    "    df = pd.DataFrame({\n",
    "        \"hid\": d[\"hhid\"],                     # ← 가구 고유ID 코드로 수정\n",
    "        \"wave\": w, \"year\": ref_year,\n",
    "        \"alc\": d[V[\"alc_home1\"]].clip(lower=0).fillna(0) + d[V[\"alc_home2\"]].clip(lower=0).fillna(0),\n",
    "        \"cons\": d[V[\"cons\"]], \"inc\": d[V[\"inc\"]],\n",
    "        \"age\": d.get(V[\"age\"]), \"hhsize\": d.get(V[\"hhsize\"]),\n",
    "        \"w\": d.get(V[\"weight_tpl\"].format(w=w), 1.0),\n",
    "    })\n",
    "    df = df[(df[\"inc\"]>0) & (df[\"cons\"]>0)]          # 결측/이상치 제거(-9 등)\n",
    "    frames.append(df)\n",
    "panel = pd.concat(frames, ignore_index=True)\n",
    "\n",
    "# ── 2. 가격·문화 변수 결합 ─────────────────────────────────\n",
    "cpi = pd.read_csv(CPI_FILE)                          # year, 소주, 맥주, 주류\n",
    "panel = panel.merge(cpi, on=\"year\", how=\"left\")\n",
    "panel[\"real_alc\"] = panel[\"alc\"] / (panel[\"주류\"]/100)   # 실질 주류소비(가격 디플레이트)\n",
    "panel[\"ln_realalc\"] = np.log(panel[\"real_alc\"].replace(0, np.nan))\n",
    "panel[\"ln_P\"]   = np.log(panel[\"소주\"])              # 소주 가격(정책 대상). 종합이면 '주류'\n",
    "panel[\"ln_inc\"] = np.log(panel[\"inc\"])\n",
    "panel[\"post\"]   = (panel[\"year\"]>=POLICY_YEAR).astype(int)   # 2024 기준판매비율 이후\n",
    "\n",
    "# 소득분위(연도 내 가중 10분위)\n",
    "panel[\"decile\"] = panel.groupby(\"year\")[\"inc\"].transform(\n",
    "    lambda s: pd.qcut(s, 10, labels=False, duplicates=\"drop\")+1)\n",
    "panel[\"low\"] = (panel[\"decile\"]<=3).astype(int)     # 저소득(1~3)\n",
    "panel[\"mid\"] = panel[\"decile\"].between(4,7).astype(int)\n",
    "panel[\"high\"]= (panel[\"decile\"]>=8).astype(int)\n",
    "\n",
    "# ★ 소비문화 추세 통제: 유연한 시간추세 + 코호트/연령/가구\n",
    "panel[\"t\"]   = panel[\"year\"]-panel[\"year\"].min()\n",
    "panel[\"t2\"]  = panel[\"t\"]**2\n",
    "panel[\"t3\"]  = panel[\"t\"]**3                         # 3차 추세(저도주화·혼술·건강트렌드 흡수)\n",
    "panel[\"cohort\"] = (panel[\"year\"]-panel[\"age\"]).round(-1)  # 출생코호트(10년) — 세대별 음주문화\n",
    "\n",
    "# ── 3. 추정 ───────────────────────────────────────────────\n",
    "panel = panel.dropna(subset=[\"ln_realalc\",\"ln_P\",\"ln_inc\",\"age\",\"hhsize\"]).set_index([\"hid\",\"year\"])\n",
    "\n",
    "# (A) 평균 가격탄력성 — 가구FE + 유연추세(문화) 통제\n",
    "mA = PanelOLS.from_formula(\n",
    "    \"ln_realalc ~ ln_P + ln_inc + age + I(age**2) + hhsize + t + t2 + t3 + EntityEffects\",\n",
    "    data=panel, weights=panel[\"w\"]).fit(cov_type=\"clustered\", cluster_entity=True)\n",
    "print(\"[A] 평균 가격탄력성 =\", round(mA.params[\"ln_P\"],3), \"(t=\", round(mA.tstats[\"ln_P\"],2), \")\")\n",
    "\n",
    "# (B) 소득계층별 가격탄력성 — ln_P × (저/중/고)\n",
    "mB = PanelOLS.from_formula(\n",
    "    \"ln_realalc ~ ln_P:C(low) + ln_P:C(mid) + ln_P:C(high) + ln_inc + age + I(age**2) + hhsize + t + t2 + t3 + EntityEffects\",\n",
    "    data=panel, weights=panel[\"w\"]).fit(cov_type=\"clustered\", cluster_entity=True)\n",
    "print(\"[B] 계층별 가격탄력성:\\n\", mB.params.filter(like=\"ln_P\"))\n",
    "\n",
    "# (C) 2024 자연실험(DiD/IV) — 기준판매비율 가격인하의 인과적 탄력성\n",
    "#     소주 상대가격 하락(post)×저소득 상호작용으로 이질적 반응 식별\n",
    "mC = PanelOLS.from_formula(\n",
    "    \"ln_realalc ~ post + post:C(low) + post:C(high) + ln_inc + age + I(age**2) + hhsize + t + t2 + EntityEffects\",\n",
    "    data=panel, weights=panel[\"w\"]).fit(cov_type=\"clustered\", cluster_entity=True)\n",
    "print(\"[C] 2024 인하효과(post):\", round(mC.params[\"post\"],3),\n",
    "      \"| 저소득 추가:\", round(mC.params.get(\"C(low)[T.1]:post\",float('nan')),3))\n",
    "\n",
    "# (D) 단기 vs 장기 — 동태(부분조정): 전기 소비 lag\n",
    "panel = panel.sort_index()\n",
    "panel[\"ln_realalc_lag\"] = panel.groupby(level=0)[\"ln_realalc\"].shift(1)\n",
    "mD = PanelOLS.from_formula(\n",
    "    \"ln_realalc ~ ln_P + ln_realalc_lag + ln_inc + age + I(age**2) + hhsize + t + t2 + EntityEffects\",\n",
    "    data=panel.dropna(subset=[\"ln_realalc_lag\"]), weights=panel[\"w\"]).fit(cov_type=\"clustered\", cluster_entity=True)\n",
    "sr = mD.params[\"ln_P\"]; rho = mD.params[\"ln_realalc_lag\"]\n",
    "print(f\"[D] 단기 탄력성={sr:.3f}, 장기 탄력성={sr/(1-rho):.3f} (조정속도 1-ρ={1-rho:.2f})\")\n",
    "\n",
    "# ── 4. 시뮬레이션 파라미터 출력 ────────────────────────────\n",
    "print(\"\\n=== 동태 세수 시뮬레이션 입력값 ===\")\n",
    "print(\"단기 ε:\", round(sr,3), \"| 장기 ε:\", round(sr/(1-rho),3),\n",
    "      \"| 저소득 ε(계층별 B/C 참조)\")\n",
    "# statsmodels 대체(가구FE=hid 더미가 크면 within-변환 후 OLS 권장)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}