# -*- coding: utf-8 -*-
"""
재정패널(NaSTaB) 가구 원자료 — 소득분위별 가격탄력성 추정
  · 목적: H2 가격탄력성 + 동태 세수 시뮬레이션 파라미터 확정
  · 식별: 2024 소주 기준판매비율(가격 인하) 자연실험 + 주종 CPI
  · 핵심 통제: 시간에 따른 '주류소비문화' 변화(저도주화·혼술·건강·고령화)를
             유연한 추세 + 코호트/연령/가구특성으로 흡수 → 가격효과 정화
  · 실행 환경: KIPF 승인환경(가구 원자료 wave 1~18, 조사 2008~2025)
  ─ 아래 변수코드/파일경로는 배포본에 맞게 수정하세요(주석 표시).
"""
import numpy as np, pandas as pd
# linearmodels 없으면: pip 불가 환경이면 statsmodels OLS+더미로 대체(하단 주석)
from linearmodels.panel import PanelOLS
import statsmodels.formula.api as smf

# ── 0. 설정 ────────────────────────────────────────────────
WAVES = range(1, 19)                 # 1~18차 (조사 2008~2025, 기준 2007~2024)
RAW = lambda w: f"nastab_{w:02d}.dta"   # ← KIPF 배포 파일명 패턴에 맞게 수정
CPI_FILE = "cpi_liquor_annual.csv"      # 연도(ref_year), 소주, 맥주, 주류종합 CPI
POLICY_YEAR = 2024                       # 기준판매비율 시행(소주 가격 인하)

# 변수코드(재정패널) ── findings 기준, 배포본에 맞게 확인/수정
V = dict(alc_home1="cd002",   # 가정용 주류비(맥주·소주·양주)
         alc_home2="cd004",   # 가정용 주류비(와인·막걸리·전통주 등)
         cons="ca002",        # 소비지출
         inc="ca004",         # 경상소득(개방형, 4차~)
         age="자녀가구주_연령",  # ← 가구주 연령 코드로 수정
         hhsize="가구원수",      # ← 가구원수 코드로 수정
         weight_tpl="h{w}cwt")   # 횡단가중치

# ── 1. 가구 패널 적재 ──────────────────────────────────────
frames=[]
for w in WAVES:
    try:
        d = pd.read_stata(RAW(w), convert_categoricals=False)
    except Exception as e:
        print(f"[wave {w}] 적재 실패: {e}"); continue
    ref_year = 2006 + w                      # wave w → 기준연도(ref) = 2006+w
    df = pd.DataFrame({
        "hid": d["hhid"],                     # ← 가구 고유ID 코드로 수정
        "wave": w, "year": ref_year,
        "alc": d[V["alc_home1"]].clip(lower=0).fillna(0) + d[V["alc_home2"]].clip(lower=0).fillna(0),
        "cons": d[V["cons"]], "inc": d[V["inc"]],
        "age": d.get(V["age"]), "hhsize": d.get(V["hhsize"]),
        "w": d.get(V["weight_tpl"].format(w=w), 1.0),
    })
    df = df[(df["inc"]>0) & (df["cons"]>0)]          # 결측/이상치 제거(-9 등)
    frames.append(df)
panel = pd.concat(frames, ignore_index=True)

# ── 2. 가격·문화 변수 결합 ─────────────────────────────────
cpi = pd.read_csv(CPI_FILE)                          # year, 소주, 맥주, 주류
panel = panel.merge(cpi, on="year", how="left")
panel["real_alc"] = panel["alc"] / (panel["주류"]/100)   # 실질 주류소비(가격 디플레이트)
panel["ln_realalc"] = np.log(panel["real_alc"].replace(0, np.nan))
panel["ln_P"]   = np.log(panel["소주"])              # 소주 가격(정책 대상). 종합이면 '주류'
panel["ln_inc"] = np.log(panel["inc"])
panel["post"]   = (panel["year"]>=POLICY_YEAR).astype(int)   # 2024 기준판매비율 이후

# 소득분위(연도 내 가중 10분위)
panel["decile"] = panel.groupby("year")["inc"].transform(
    lambda s: pd.qcut(s, 10, labels=False, duplicates="drop")+1)
panel["low"] = (panel["decile"]<=3).astype(int)     # 저소득(1~3)
panel["mid"] = panel["decile"].between(4,7).astype(int)
panel["high"]= (panel["decile"]>=8).astype(int)

# ★ 소비문화 추세 통제: 유연한 시간추세 + 코호트/연령/가구
panel["t"]   = panel["year"]-panel["year"].min()
panel["t2"]  = panel["t"]**2
panel["t3"]  = panel["t"]**3                         # 3차 추세(저도주화·혼술·건강트렌드 흡수)
panel["cohort"] = (panel["year"]-panel["age"]).round(-1)  # 출생코호트(10년) — 세대별 음주문화

# ── 3. 추정 ───────────────────────────────────────────────
panel = panel.dropna(subset=["ln_realalc","ln_P","ln_inc","age","hhsize"]).set_index(["hid","year"])

# (A) 평균 가격탄력성 — 가구FE + 유연추세(문화) 통제
mA = PanelOLS.from_formula(
    "ln_realalc ~ ln_P + ln_inc + age + I(age**2) + hhsize + t + t2 + t3 + EntityEffects",
    data=panel, weights=panel["w"]).fit(cov_type="clustered", cluster_entity=True)
print("[A] 평균 가격탄력성 =", round(mA.params["ln_P"],3), "(t=", round(mA.tstats["ln_P"],2), ")")

# (B) 소득계층별 가격탄력성 — ln_P × (저/중/고)
mB = PanelOLS.from_formula(
    "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",
    data=panel, weights=panel["w"]).fit(cov_type="clustered", cluster_entity=True)
print("[B] 계층별 가격탄력성:\n", mB.params.filter(like="ln_P"))

# (C) 2024 자연실험(DiD/IV) — 기준판매비율 가격인하의 인과적 탄력성
#     소주 상대가격 하락(post)×저소득 상호작용으로 이질적 반응 식별
mC = PanelOLS.from_formula(
    "ln_realalc ~ post + post:C(low) + post:C(high) + ln_inc + age + I(age**2) + hhsize + t + t2 + EntityEffects",
    data=panel, weights=panel["w"]).fit(cov_type="clustered", cluster_entity=True)
print("[C] 2024 인하효과(post):", round(mC.params["post"],3),
      "| 저소득 추가:", round(mC.params.get("C(low)[T.1]:post",float('nan')),3))

# (D) 단기 vs 장기 — 동태(부분조정): 전기 소비 lag
panel = panel.sort_index()
panel["ln_realalc_lag"] = panel.groupby(level=0)["ln_realalc"].shift(1)
mD = PanelOLS.from_formula(
    "ln_realalc ~ ln_P + ln_realalc_lag + ln_inc + age + I(age**2) + hhsize + t + t2 + EntityEffects",
    data=panel.dropna(subset=["ln_realalc_lag"]), weights=panel["w"]).fit(cov_type="clustered", cluster_entity=True)
sr = mD.params["ln_P"]; rho = mD.params["ln_realalc_lag"]
print(f"[D] 단기 탄력성={sr:.3f}, 장기 탄력성={sr/(1-rho):.3f} (조정속도 1-ρ={1-rho:.2f})")

# ── 4. 시뮬레이션 파라미터 출력 ────────────────────────────
print("\n=== 동태 세수 시뮬레이션 입력값 ===")
print("단기 ε:", round(sr,3), "| 장기 ε:", round(sr/(1-rho),3),
      "| 저소득 ε(계층별 B/C 참조)")
# statsmodels 대체(가구FE=hid 더미가 크면 within-변환 후 OLS 권장)