# -*- coding: utf-8 -*-
"""
가계동향조사 원자료 — 2024 소주 기준판매비율의 가구 파급(미시)
  · 목적: 가격인하 → 가정 주류소비 + 외식 파급 + 채널 전환(가정↔외식)
  · 식별: 2024Q1 소주만 가격하락(맥주·외식 미반영) 자연실험, 반복횡단 DiD/이벤트
  · 한계: 외식 주류는 외식비에 합산 → #2 외식비·#3 주점(분리 시)로 근사
  · 실행 환경: 통계청 MDIS 승인환경(가계동향조사 마이크로데이터, 분기)
  ─ 변수코드/경로는 MDIS 배포본에 맞게 수정(주석).
"""
import numpy as np, pandas as pd, statsmodels.formula.api as smf

YEARS=range(2018,2025); FILE=lambda y,q:f"hies_{y}Q{q}.csv"
CPI="cpi_liquor_quarterly.csv"      # year,q,소주,맥주,주류,소주외식
POLICY=(2024,1)
V=dict(home_alc="주류지출",         # 02.1 가정용 주류(분리 가능)
       dining="외식비",             # 11.1 음식점·주점(외식)
       pub="주점지출",              # (분리 제공 시) 외식 주류 근사 — MDIS 확인
       beer="맥주지출",             # 위약대조(2024 미변경)
       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)
h=pd.concat(rows,ignore_index=True)

# 2) 변수
h["post"]=((h.year>POLICY[0])|((h.year==POLICY[0])&(h.q>=POLICY[1]))).astype(int)
h=h[(h[V["inc"]]>0)&(h[V["cons"]]>0)]
h["ln_inc"]=np.log(h[V["inc"]])
h["low"]=(h[V["decile"]]<=3).astype(int); h["high"]=(h[V["decile"]]>=8).astype(int)
h["t"]=(h.year*4+h.q); h["t"]=h["t"]-h["t"].min(); h["t2"]=h["t"]**2
h["covid"]=h.year.isin([2020,2021]).astype(int)
# (선택) 실질화: 각 지출을 해당 CPI로 디플레이트 후 로그

CTRL=f"ln_inc + {V['age']} + I({V['age']}**2) + C({V['edu']}) + {V['hhsize']} + C({V['region']}) + C(q) + t + t2 + covid"

def run(dep,label):
    d=h.copy(); d["y"]=np.log(d[dep].replace(0,np.nan))   # 0은 결측(2단계모형 권장) 또는 share 사용
    # (A) 평균 파급
    mA=smf.wls(f"y ~ post + {CTRL}", data=d, weights=d[V['weight']]).fit(cov_type="HC1")
    # (B) 분위 이질성
    mB=smf.wls(f"y ~ post + post:C(low) + post:C(high) + {CTRL}", data=d, weights=d[V['weight']]).fit(cov_type="HC1")
    print(f"[{label}] 평균 post={mA.params['post']:+.3f}(t={mA.tvalues['post']:.2f}) "
          f"| 저소득추가={mB.params.get('post:C(low)[T.1]',float('nan')):+.3f} 고소득={mB.params.get('post:C(high)[T.1]',float('nan')):+.3f}")
    return mA

print("=== 2024 소주 인하 가구 파급 (가중 WLS, 통제·코로나더미) ===")
mHome=run(V["home_alc"],"가정 주류지출")          # 가격인하 직접반응(가정 채널)
run(V["dining"],"외식비(음식점)")                 # 외식 전반 파급(간접)
if V["pub"] in h.columns: run(V["pub"],"주점지출")  # 외식 주류 근사(분리 시)
run(V["beer"],"맥주지출(위약대조)")               # 2024 미변경 → 유의하면 다른 요인 경보

# 탄력성 환산(가정 주류): post계수 ÷ 소주 소매가 변화율
DPRICE=-0.063
print(f"\n가정 주류 가격탄력성 ≈ {mHome.params['post']/DPRICE:+.2f} (소매 -6.3% 가정)")

# 강건성 골격
print("\n[강건성] ① 코로나 제외(2022~) ② 플라시보(post=2023Q1) ③ 이벤트스터디(연도×분기 lead/lag) "
      "④ 비음용가구 대조 ⑤ 0지출 2단계(Heckman)·예산몫(share) 회귀")