from pathlib import Path
import pandas as pd, numpy as np, json
import statsmodels.formula.api as smf
P=Path(__file__).resolve().parent
src=P/'usu_individual_T425.txt'
if not src.exists(): src=Path('/mnt/data/usu_individual_T425.txt')
d=pd.read_csv(src,sep=';',decimal=',',low_memory=False)
assert (d.ANO4.eq(2025)&d.TRIMESTRE.eq(4)).all()
a=d.loc[d.ESTADO.eq(1)&d.CH06.between(18,65)&d.P21.gt(0)&d.PONDIIO.gt(0)&d.NIVEL_ED.between(1,6)&d.CH04.isin([1,2])].copy()
a['log_ingreso']=np.log(a.P21);a['edad']=a.CH06.astype(float);a['edad2']=a.edad**2;a['mujer']=(a.CH04==2).astype(int);a['educ']=a.NIVEL_ED.astype(int);a['region']=a.REGION.astype(str);a['peso']=a.PONDIIO.astype(float);a['hogar']=a.CODUSU.astype(str)+'_'+a.NRO_HOGAR.astype(str)
labels={1:'Primaria incompleta',2:'Primaria completa',3:'Secundaria incompleta',4:'Secundaria completa',5:'Superior incompleto',6:'Superior completo'}
regions={'1':'Gran Buenos Aires','40':'NOA','41':'NEA','42':'Cuyo','43':'Pampeana','44':'Patagonia'}
def wmedian(g):
 x=g.P21.to_numpy();w=g.peso.to_numpy();o=np.argsort(x);return float(x[o][np.searchsorted(np.cumsum(w[o]),w.sum()/2)])
def stats(g):return {'n':int(len(g)),'media':round(float(np.average(g.P21,weights=g.peso))),'mediana':round(wmedian(g))}
def fit(data,formula,weighted=True):
 model=(smf.wls(formula,data=data,weights=data.peso) if weighted else smf.ols(formula,data=data)).fit(cov_type='cluster',cov_kwds={'groups':data.hogar})
 terms=[]
 for k,v in model.params.items():
  lo,hi=model.conf_int().loc[k];terms.append({'variable':k,'coef':float(v),'li':float(lo),'ls':float(hi),'p':float(model.pvalues[k]),'efecto_pct':float(100*np.expm1(v)) if k!='Intercept' else None})
 return {'n':int(model.nobs),'r2':round(float(model.rsquared),4),'coeficientes':terms}
formula='log_ingreso ~ C(educ, Treatment(reference=4)) + edad + edad2 + mujer + C(region)'
main=fit(a,formula);unweighted=fit(a,formula,False)
# Hours worked in all occupations, positive and plausible; secondary model only
h=a.loc[a.PP3E_TOT.between(1,100)].copy();h['log_horas']=np.log(h.PP3E_TOT);h['categoria']=h.CAT_OCUP.astype(str)
hours=fit(h,formula+' + log_horas + C(categoria)')
byedu=[{'codigo':k,'educacion':v,**stats(a.loc[a.educ.eq(k)])} for k,v in labels.items()]
bysex=[{'grupo':v,**stats(a.loc[a.mujer.eq(k)])} for k,v in [(0,'Varones'),(1,'Mujeres')]]
byreg=[{'grupo':regions.get(k,k),**stats(g)} for k,g in a.groupby('region')]
# Education by sex cross-tab for visual comparisons
cross=[{'educacion':labels[k],'sexo':'Mujeres' if s else 'Varones',**stats(g)} for (k,s),g in a.groupby(['educ','mujer'])]
def find(model,token):
 for x in model['coeficientes']:
  if token in x['variable']:return x
 return None
ref=find(main,'[T.6]');female=find(main,'mujer')
res={'periodo':'2025-T4','total_personas':int(len(d)),'muestra':int(len(a)),'hogares_muestra':int(a.hogar.nunique()),'modelo':main,'sin_ponderar':unweighted,'con_horas_categoria':hours,'educacion':byedu,'sexo':bysex,'region':byreg,'educacion_sexo':cross,'principal_superior':ref,'principal_mujer':female,'filtros':'Ocupados 18–65 años, P21 > 0, PONDIIO > 0, NIVEL_ED 1–6, CH04 1–2','nota_horas':'PP3E_TOT: horas trabajadas en todas las ocupaciones; no son necesariamente horas de la ocupación principal. Submuestra 1–100 horas.'}
(P/'resultados.json').write_text(json.dumps(res,ensure_ascii=False,indent=2),encoding='utf-8')
pd.DataFrame(byedu).to_csv(P/'ingresos_educacion.csv',index=False)
pd.DataFrame(main['coeficientes']).to_csv(P/'coeficientes.csv',index=False)
print('N',len(a),'R2',main['r2'],'Robust N',hours['n'],'R2',hours['r2'])
print('Superior vs secundaria completa',ref)
print('Mujer',female)