# -*- coding: utf-8 -*-
"""
가계동향조사(통계청) 원자료 — 소득분위별 소주 가격탄력성 (2024 기준판매비율 자연실험)
  · 목적: H2 미시 탄력성 2차 검증(재정패널과 교차) + 시뮬레이션 파라미터 확정
  · 식별: 2024Q1 소주 기준판매비율(가격 인하)을 사건으로 한 DiD/이벤트
  · 핵심 통제: 소비문화 변화(저도주화·혼술·건강·고령화)를 유연추세+가구특성으로 흡수
  · 실행 환경: 통계청 MDIS 승인환경(가계동향조사 마이크로데이터, 분기)
  ─ 변수코드/파일경로는 MDIS 배포본에 맞게 수정(주석 표시).
"""
import numpy as np, pandas as pd
import statsmodels.formula.api as smf

# ── 0. 설정 ────────────────────────────────────────────────
YEARS = range(2018, 2025)                 # 2018~2024(분기). 필요시 확장
FILE = lambda y,q: f"hies_{y}Q{q}.csv"     # ← MDIS 파일명 패턴에 맞게 수정
CPI_FILE = "cpi_liquor_quarterly.csv"      # ym(YYYYQ), 소주, 맥주, 주류 CPI(분기)
POLICY = (2024,1)                          # 기준판매비율 시행 시점(2024Q1)

# 가계동향 변수코드(예시) ── MDIS 배포본 코드로 반드시 확인/수정
V = dict(alc="주류지출",        # 주류 지출(COICOP 02.1). 소주 분리 가능하면 소주지출 사용
         soju="소주지출",        # (분리 제공 시) 소주 단독 지출
         beer="맥주지출",        # (분리 제공 시) 맥주 단독 — DiD 대조군
         cons="소비지출", inc="경상소득", decile="소득10분위",
         age="가구주연령", edu="가구주교육", hhsize="가구원수",
         region="지역", weight="가중값")

# ── 1. 적재(분기 반복횡단) ─────────────────────────────────
rows=[]
for y in YEARS:
    for q in (1,2,3,4):
        try: d=pd.read_csv(FILE(y,q))
        except Exception: continue
        d["year"],d["q"]=y,q
        rows.append(d)
hies=pd.concat(rows,ignore_index=True)

# ── 2. 변수 구성 ───────────────────────────────────────────
cpi=pd.read_csv(CPI_FILE)                                  # 분기 CPI
hies["yq"]=hies["year"]*10+hies["q"]
cpi["yq"]=cpi["year"]*10+cpi["q"]
hies=hies.merge(cpi[["yq","소주","맥주","주류"]],on="yq",how="left")
# 실질 주류/소주 지출(가격 디플레이트) — 지출>0만 로그, 미구입=결측 or 2단계(Tobit/Heckman)
hies["real_soju"]=hies[V.get("soju",V["alc"])]/(hies["소주"]/100)
hies["ln_soju"]=np.log(hies["real_soju"].replace(0,np.nan))
hies["ln_inc"]=np.log(hies[V["inc"]].clip(lower=1))
hies["post"]=((hies["year"]>POLICY[0])|((hies["year"]==POLICY[0])&(hies["q"]>=POLICY[1]))).astype(int)
hies["low"]=(hies[V["decile"]]<=3).astype(int)
hies["high"]=(hies[V["decile"]]>=8).astype(int)
# 소비문화 추세: 연속 시간 + 2·3차 + 계절(분기) + 가구특성
hies["t"]=hies["yq"].rank(method="dense").astype(int); hies["t2"]=hies["t"]**2; hies["t3"]=hies["t"]**3

base = ("ln_soju ~ ln_inc + {age} + I({age}**2) + C({edu}) + {hh} + C({reg}) "
        "+ C(q) + t + t2 + t3").format(age=V["age"],edu=V["edu"],hh=V["hhsize"],reg=V["region"])

# ── 3. 추정 ───────────────────────────────────────────────
# (A) 2024 인하 이벤트 효과(소주 실질소비 변화) — 가중 OLS, 지역·분기·문화추세 통제
mA=smf.wls(base+" + post", data=hies, weights=hies[V["weight"]]).fit(cov_type="HC1")
print("[A] 2024 인하 후 소주소비 변화(post):", round(mA.params["post"],3),"(t=",round(mA.tvalues["post"],2),")")

# (B) 소득계층 이질성 — post×저/고소득
mB=smf.wls(base+" + post + post:C(low) + post:C(high)", data=hies, weights=hies[V["weight"]]).fit(cov_type="HC1")
print("[B] 저소득 추가반응 post:C(low):", round(mB.params.get("post:C(low)[T.1]",float('nan')),3),
      "| 고소득 post:C(high):", round(mB.params.get("post:C(high)[T.1]",float('nan')),3))

# (C) 소주 vs 맥주 이중차분 (품목 분리 제공 시) — 가장 깨끗한 인과 탄력성
#     맥주는 2024 세제변화 없음 → 대조군. 소주 상대가격 하락효과 식별
if V.get("beer") in hies.columns and V.get("soju") in hies.columns:
    long=hies.melt(id_vars=["year","q","post","ln_inc",V["age"],V["edu"],V["hhsize"],V["region"],V["decile"],V["weight"],"t","t2"],
                   value_vars=[V["soju"],V["beer"]],var_name="item",value_name="exp")
    long["soju"]=(long["item"]==V["soju"]).astype(int)
    long["ln_exp"]=np.log(long["exp"].replace(0,np.nan))
    mC=smf.wls(f"ln_exp ~ post*soju + ln_inc + {V['age']} + C({V['edu']}) + {V['hhsize']} + C({V['region']}) + C(q) + t + t2",
               data=long, weights=long[V["weight"]]).fit(cov_type="HC1")
    print("[C-DiD] 소주×post(순수 소주 가격인하 효과):", round(mC.params.get("post:soju",float('nan')),3))
    # → 가격변화율(소매 -6.3%)로 나누면 탄력성: ε ≈ 계수 / (-0.063)

# ── 4. 탄력성 환산 ────────────────────────────────────────
dprice=-0.063   # 2024 소주 소매가 변화율(관측). 실제 추정 소매가로 대체
print("\n=== 탄력성 환산 (ε = Δln수요 / Δln가격) ===")
print("전체 ε ≈", round(mA.params['post']/dprice,3),
      "| 저소득 ε ≈", round((mA.params['post']+mB.params.get('post:C(low)[T.1]',0))/dprice,3))
print("※ 미구입(0) 처리: 2단계모형(Heckman/Tobit) 또는 예산몫(share) 회귀로 강건성 확인 권장")