# -*- coding: utf-8 -*-
"""3재화 LA-AIDS — 소주·국산맥주·수입맥주 (종가세기 2006~2019, 국세 물량·과세표준)
 동차성 부과(소주=numeraire), Marshallian 자체/교차/지출 탄력성. 소주 수입분은 무시(연 수십L).
 주의: n=14 연간·집계 → 정밀도 한계. 제품·용기·월별·기온통제 완전판은 스캐너 필요.
"""
import sys,io,json; sys.stdout=io.TextIOWrapper(sys.stdout.buffer,encoding='utf-8')
import numpy as np, pandas as pd, statsmodels.formula.api as smf
B='C:/Temp/주류세/공개데이터'
g=pd.read_csv(B+'/국세청_국세통계/주세신고_주종별_제조수입_2006_2025.csv')
g['값']=pd.to_numeric(g['값'].astype(str).str.replace(',',''),errors='coerce')
def P(주종,지표): return g[(g.주종==주종)&(g.지표==지표)].pivot_table(index='연도',columns='구분',values='값')
vB,xB=P('맥주','출고량'),P('맥주','과세표준_종가세')
vS,xS=P('희석식소주','출고량'),P('희석식소주','과세표준_종가세')
d=pd.DataFrame(index=vB.index)
d['q_impb']=vB['수입분']; d['q_domb']=vB['합계']-vB['수입분']; d['q_soju']=vS['합계']
d['x_impb']=xB['수입분']; d['x_domb']=xB['합계']-xB['수입분']; d['x_soju']=xS['합계']
d=d.dropna(); d=d[(d.index>=2006)&(d.index<=2019)]
for k in ['soju','domb','impb']:
    d['p_'+k]=d['x_'+k]/d['q_'+k]
d['X']=d.x_soju+d.x_domb+d.x_impb
for k in ['soju','domb','impb']: d['w_'+k]=d['x_'+k]/d.X
d['lnP']=sum(d['w_'+k]*np.log(d['p_'+k]) for k in ['soju','domb','impb'])
d['ln_XP']=np.log(d.X)-d.lnP
d['r_domb']=np.log(d.p_domb/d.p_soju); d['r_impb']=np.log(d.p_impb/d.p_soju)
print('===== 지출몫 2006~2019 (소주/국산맥주/수입맥주) =====')
for yr in [d.index.min(), d.index.max()]:
    r=d.loc[yr]; print(f'  {yr}: 소주 {r.w_soju:.3f} · 국산맥주 {r.w_domb:.3f} · 수입맥주 {r.w_impb:.3f}')
print(f'  평균: 소주 {d.w_soju.mean():.3f} · 국산맥주 {d.w_domb.mean():.3f} · 수입맥주 {d.w_impb.mean():.3f}')

goods=['soju','domb','impb']; names={'soju':'소주','domb':'국산맥주','impb':'수입맥주'}
G={}; Bc={}
for i in goods:
    m=smf.ols(f'w_{i} ~ r_domb + r_impb + ln_XP', data=d).fit()
    g_domb=m.params['r_domb']; g_impb=m.params['r_impb']; g_soju=-(g_domb+g_impb)
    G[i]={'soju':g_soju,'domb':g_domb,'impb':g_impb}; Bc[i]=m.params['ln_XP']
wbar={k:d['w_'+k].mean() for k in goods}
def elas(i,j): return (-1.0 if i==j else 0.0)+(G[i][j]-Bc[i]*wbar[j])/wbar[i]
print('\n===== Marshallian 가격탄력성 행렬 (행=수요, 열=가격) =====')
hdr='수요\\가격     '+''.join(f'{names[j]:>10}' for j in goods)
print(hdr)
for i in goods:
    row=f'{names[i]:<10}'+''.join(f'{elas(i,j):>10.2f}' for j in goods)
    print(row)
print('\n===== 자체가격·지출 탄력성 요약 =====')
res={}
for i in goods:
    own=elas(i,i); eta=1+Bc[i]/wbar[i]
    res[names[i]]=dict(own=round(own,2),expenditure=round(eta,2))
    print(f'  {names[i]}: 자체가격 {own:+.2f} · 지출 {eta:+.2f}')
print('\n[원논문 대조] 수입프리미엄 -4.8~-6.0 / 국산코어 -1.8~-1.9 / (소주는 원논문 대상 아님)')
print('[선행 소주탄력성 문헌범위] 대략 -0.3~-0.9 (증류주)')
json.dump({'own_expend':res,'elast_matrix':{names[i]:{names[j]:round(elas(i,j),3) for j in goods} for i in goods},
           'wbar':{names[k]:round(wbar[k],3) for k in goods},'n':len(d)},
          open(B+'/문헌/han2018_aids_3good.json','w',encoding='utf-8'),ensure_ascii=False,indent=2)
print('\n저장: 문헌/han2018_aids_3good.json')