/*******************************************************************************
    THIS DO FILE GENERATES THE MAIN RESULTS (GRAPHS AND TABLES) IN THE PAPER
*******************************************************************************/

clear all
set more off

* Directories
if c(username)=="aaguilares" {
	global panel "D:\Impuestos\Data\Standardized datasets"
	global graphs "D:\Impuestos\Tax_paper_AEJ\Graphs"
	global tables "D:\Impuestos\Tax_paper_AEJ\Tables"
	}
if c(username)=="aaeag" {
	global panel "C:\Users\aaeag\Dropbox (QFPD)\Kantar\Data"
	global graphs "C:\Users\aaeag\Dropbox (QFPD)\Kantar\AEJ_resubmission_AAE\Graphs"
	global tables "C:\Users\aaeag\Dropbox (QFPD)\Kantar\AEJ_resubmission_AAE\Tables"
	}

* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------
* SET OF RESULTS BASED ON THE HOUSEHOLD PANEL 2013-2014
* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------

use "$panel/Panel hogares/Panel_HH_week_FINAL_1314.dta", clear

* Time variables created
g year=year(dofw(fecha_sem))
g month=month(dofw(fecha_sem))
tab week, g(wk)

* Filter to select households that appear at least once in both 2013 and 2014
g yy0=(year==2012)
g yy1=(year==2013)
g yy2=(year==2014)
bysort iddomicilio: egen y2012=max(yy0)
bysort iddomicilio: egen y2013=max(yy1)
bysort iddomicilio: egen y2014=max(yy2)

* Define Outcomes
global Y_prices "index_P_index index_P_index_TaxBev index_P_index_NoTaxBev index_P_index_TaxFood index_P_index_NoTaxFood"
global Y_calories "kcal_tot SD_kcal nonSD_kcal HCF_kcal nonHCF_kcal"
global Y_nutrients "azucar_tot grasas_tot grasasat_tot carbohidratos_tot colesterol_tot sodio_tot proteinas_tot hierro_tot"
global Y_robustness "index_P_index kcal_tot index_P_index_TaxBev SD_kcal index_P_index_TaxFood HCF_kcal"
global Y_liters "SD_litros nonSD_litros"
global Y_kilos "HCF_kilos nonHCF_kilos"

* Running variable and treatment variable
g G=fecha_sem-wofd(mdy(1,1,2014))
g G2=G^2
g G3=G^3
g G4=G^4
g T=(G>=0)
g T_G=T*G
g T_G2=T*G2
g T_G3=T*G3
g T_G4=T*G4


* Keeping relevant observations
* Years
keep if year==2013 | year==2014
* Households
keep if y2013==1 & y2014==1

* Partial out week dummies and HH FE
foreach var of varlist $Y_prices {
	xtreg `var' wk2-wk52, fe
	predict `var'_res, resid
	}
foreach var of varlist $Y_calories $Y_nutrients $Y_liters $Y_kilos {
	g log_`var'=log(`var'+1)
	xtreg log_`var' wk2-wk52, fe
	predict `var'_res, resid
	}

* Create weekly averages by outcome
* Here we define how many weeks to aggregate in the conditional means for the scatter
g G_r=floor(G/4)
foreach var of varlist $Y_prices {
	bysort G: egen `var'_MNP=mean(`var')
	}
foreach var of varlist $Y_calories {
	bysort G: egen `var'_MNP=mean(`var')
	}
foreach var of varlist $Y_prices $Y_calories $Y_nutrients $Y_liters $Y_kilos {
	bysort G_r: egen `var'_MN=mean(`var'_res)
	}

* Rescale it to avoid many dots in the scatterplot
replace G_r=G_r*4

reg kcal_tot_res T G
g muestra=(e(sample)==1)

/* ------------------------------------------------------
                TABLE 1: DESCRIPTIVE STATS
   ------------------------------------------------------ */
preserve
keep if muestra==1
g gasto_ENIGH=4*gasto
g SD_ENIGH=4*SD_gasto
g Botana_ENIGH=4*HCF_gasto
cd "$tables\Paper"
sutex gasto_ENIGH SD_ENIGH Botana_ENIGH, file(Table1_Descriptive_stats_PartB.tex) replace
collapse (mean) edadjefefamilia hombre_jf prim sec terc imcjf obese overw perc_tv perc_comp perc_auto ses, by(iddomicilio year)
keep if year==2013
sutex edadjefefamilia hombre_jf prim sec terc imcjf obese overw perc_tv perc_comp perc_auto, file(Table1_Descriptive_stats_PartA.tex) replace


/* ------------------------------------------------------
      TABLE 2: RDROBUST CALORIES INCLUDING 2013 Y 2014
   ------------------------------------------------------ */
restore
local contador=1
foreach var of varlist $Y_prices $Y_calories {
	rdrobust `var'_res G if muestra==1, c(0) p(2) kernel(triangular) h(52)
	if `contador'==1 {
	outreg2 using "$tables\Paper\Table2_RD_robust.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes replace
		}
	if `contador'>1 & `contador'<10 {
	outreg2 using "$tables\Paper\Table2_RD_robust.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	if `contador'==10 {
	outreg2 using "$tables\Paper\Table2_RD_robust.tex", tex ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	local contador=`contador'+1
	}


/* ------------------------------------------------------
      TABLE 3: RDROBUST NUTRIENTS INCLUDING 2013 Y 2014
   ------------------------------------------------------ */
	
local contador=1
foreach var of varlist kcal_tot $Y_nutrients {
	rdrobust `var'_res G, c(0) p(2) kernel(triangular) h(52)
	if `contador'==1 {
	outreg2 using "$tables\Paper\Table3_nutrients_RD_robust.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes replace
		}
	if `contador'>1 & `contador'<10 {
	outreg2 using "$tables\Paper\Table3_nutrients_RD_robust.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
	}
	if `contador'==8 {
	outreg2 using "$tables\Paper\Table3_nutrients_RD_robust.tex", tex ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	local contador=`contador'+1
	}

/* ------------------------------------------------------
        TABLE 4: HETEROGENEITY BY BMI AND SES
   ------------------------------------------------------ */

* Panel A: Baseline Parametric results

gl Y_heterog "index_P_index index_P_index_TaxBev index_P_index_TaxFood kcal_tot SD_kcal HCF_kcal"

local contador=1
foreach var of varlist $Y_heterog {
	rdrobust `var'_res G, c(0) p(2) kernel(triangular) h(52)
	if `contador'==1 {
	outreg2 using "$tables\Paper\Table4_heterog_base.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes replace
		}
	if `contador'>1 & `contador'<6 {
	outreg2 using "$tables\Paper\Table4_heterog_base.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	if `contador'==6 {
	outreg2 using "$tables\Paper\Table4_heterog_base.tex", tex ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	local contador=`contador'+1
	}

* Panel B: Heterogeneity by BMI

forvalues i=1/3 {
local contador=1
foreach var of varlist $Y_heterog {
	rdrobust `var'_res G if D_obesity==`i', c(0) p(2) kernel(triangular) h(52)
	if `contador'==1 {
	outreg2 using "$tables\Paper\TableA1_heterog_panelB`i'.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes replace
		}
	if `contador'>1 & `contador'<6 {
	outreg2 using "$tables\Paper\TableA1_heterog_panelB`i'.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	if `contador'==6 {
	outreg2 using "$tables\Paper\TableA1_heterog_panelB`i'.tex", tex ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	local contador=`contador'+1
	}
}

* Panel C: Heterogeneity by SES

forvalues i=1/3 {
local contador=1
foreach var of varlist $Y_heterog {
	rdrobust `var'_res G if D_ses==`i', c(0) p(2) kernel(triangular) h(52)
	if `contador'==1 {
	outreg2 using "$tables\Paper\TableA1_heterog_panelC`i'.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes replace
		}
	if `contador'>1 & `contador'<6 {
	outreg2 using "$tables\Paper\TableA1_heterog_panelC`i'.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	if `contador'==6 {
	outreg2 using "$tables\Paper\TableA1_heterog_panelC`i'.tex", tex ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	local contador=`contador'+1
	}
}


/* ------------------------------------------------------
              FIGURE 1: PRICE GRAPHS
   ------------------------------------------------------ */

foreach var of varlist $Y_prices {
twoway (qfitci `var'_res G if G<0, clcolor(navy) clwidth(medium)) ///
		(qfitci `var'_res G if G>=0, clcolor(navy) clwidth(medium)) ///
		(scatter `var'_MN G_r), ///
        yscale(range(10 10)) ylabel(-10(5)10) ytitle(Price index) ///
		xscale(range(-52 52)) xtitle(Weeks Before / After Tax) ///
		xline(0, lwidth(medthin) lpattern(solid) lcolor(black)) legend(off) ///
		graphregion(fcolor(white) ifcolor(white))
	graph export "$graphs\Paper\Figure1B_`var'.png", as(png) replace
	}

/* ------------------------------------------------------
        FIGURE 2: QFIT CALORIES 
   ------------------------------------------------------ */

foreach var of varlist $Y_calories {
twoway (qfitci `var'_res G if G<0, clcolor(navy) clwidth(medium)) ///
		(qfitci `var'_res G if G>=0, clcolor(navy) clwidth(medium)) ///
		(scatter `var'_MN G_r), ///
        yscale(range(-0.2 0.2)) ylabel(-0.2(0.05)0.2) ytitle(Calories) ///
		xscale(range(-52 52)) xtitle(Weeks Before / After Tax) ///
		xline(0, lwidth(medthin) lpattern(solid) lcolor(black)) legend(off) ///
		graphregion(fcolor(white) ifcolor(white))
	graph export "$graphs\Paper\Figure2_`var'.png", as(png) replace
	}

/* ------------------------------------------------------
        FIGURE 3: QFIT NUTRIENTS 
   ------------------------------------------------------ */

foreach var of varlist $Y_nutrients {
twoway (qfitci `var'_res G if G<0, clcolor(navy) clwidth(medium)) ///
		(qfitci `var'_res G if G>=0, clcolor(navy) clwidth(medium)) ///
		(scatter `var'_MN G_r), ///
        yscale(range(-0.2 0.2)) ylabel(-0.2(0.05)0.2) ytitle(Nutrients) ///
		xscale(range(-52 52)) xtitle(Weeks Before / After Tax) ///
		xline(0, lwidth(medthin) lpattern(solid) lcolor(black)) legend(off) ///
		graphregion(fcolor(white) ifcolor(white))
	graph export "$graphs\Paper\Figure3_`var'.png", as(png) replace
	}

/* ------------------------------------------------------
              TABLE A.10: ANTICIPATORY EFFECTS
   ------------------------------------------------------ */

foreach var of varlist $Y_robustness {
	forvalues ll=0/3 {
	reg `var'_res T G G2 T_G T_G2 if year>=2013 & year<=2014 & abs(G)>`ll', vce(cluster iddomicilio)
		if `ll'==0 {
		outreg2 T using "$tables\Paper\TableA10_Anticipatory_`var'.txt", text bdec(3) sdec(4) ct(`var'_`ll') ///
			keep(T) nocons replace
		}
		if `ll'>0 & `ll'<3 {
		outreg2 T using "$tables\Paper\TableA10_Anticipatory_`var'.txt", text bdec(3) sdec(4) ct(`var'_`ll') ///
			keep(T) nocons append
		}
		if `ll'==3 {
		outreg2 T using "$tables\Paper\TableA10_Anticipatory_`var'.tex", tex bdec(3) sdec(4) ct(`var'_`ll') ///
			keep(T) nocons append
		}
	}
}

/* ------------------------------------------------------
    TABLE A.3: SENSITIVITY PARAMETRIC VS NON PARAMETRIC
   ------------------------------------------------------ */

* Parametric Degree 2
local contador=1
foreach var of varlist $Y_prices $Y_calories {
	reg `var'_res T G G2 T_G T_G2 if year>=2013 & year<=2014, vce(cluster iddomicilio)
	if `contador'==1 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_2.txt", text bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons replace
		}
	if `contador'>1 & `contador'<10 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_2.txt", text bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons append
		}
	if `contador'==10 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_2.tex", tex bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons append
		}
	local contador=`contador'+1
	}

* Parametric Degree 3
local contador=1
foreach var of varlist $Y_prices $Y_calories {
	reg `var'_res T G G2 G3 T_G T_G2 T_G3  if year>=2013 & year<=2014, vce(cluster iddomicilio)
	if `contador'==1 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_3.txt", text bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons replace
		}
	if `contador'>1 & `contador'<10 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_3.txt", text bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons append
		}
	if `contador'==10 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_3.tex", tex bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons append
		}
	local contador=`contador'+1
	}	

* Parametric Degree 4
local contador=1
foreach var of varlist $Y_prices $Y_calories {
	reg `var'_res T G G2 G3 G4 T_G T_G2 T_G3 T_G4 if year>=2013 & year<=2014, vce(cluster iddomicilio)
	if `contador'==1 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_4.txt", text bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons replace
		}
	if `contador'>1 & `contador'<10 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_4.txt", text bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons append
		}
	if `contador'==10 {
	outreg2 T using "$tables\Paper\TableA3_Parametric_4.tex", tex bdec(3) sdec(4) ct(`var') ///
		keep(T) nocons append
		}
	local contador=`contador'+1
	}	

	
/* -----------------------------------------------------------------
    FIGURE A.1: RELATION BETWEEN BMI AND CALORIES FROM SD AND HCF
   ----------------------------------------------------------------- */

xtile pct_bmi = imcjf, n(10)
bysort pct_bmi: egen MN_BMI=mean(imcjf)
bysort pct_bmi: egen MN_calSD=mean(SD_kcal)
bysort pct_bmi: egen MN_calHCF=mean(HCF_kcal)

twoway (lpolyci SD_kcal imcjf, fintensity(inten10)) ///
	   (scatter MN_calSD MN_BMI, msymbol(circle_hollow) fintensity(inten100)   ),  ///
		ytitle("SD calories", axis(1)) xtitle("BMI", axis(1) ) ///
		graphregion(color(white)) plotregion(lcolor(black) ) ///
		legend( label(3 "Weekly purchased calories from SDs") label(1 "Confidence Interval 95%") label(2 "Local Polynomial") ///
		cols(2) symxsize(10) symplacement(center) region(lcolor(white)) )

graph export "$graphs\Paper\FigureA1_BMI_SDcalories.png", as(png) replace	

twoway (lpolyci HCF_kcal imcjf, fintensity(inten10)) ///
	   (scatter MN_calHCF MN_BMI, msymbol(circle_hollow) fintensity(inten100)   ),  ///
		ytitle("HCF calories", axis(1)) xtitle("BMI", axis(1) ) ///
		graphregion(color(white)) plotregion(lcolor(black) ) ///
		legend( label(3 "Weekly purchased calories from HCFs") label(1 "Confidence Interval 95%") label(2 "Local Polynomial") ///
		cols(2) symxsize(10) symplacement(center) region(lcolor(white)) )

graph export "$graphs\Paper\FigureA1_BMI_HCFcalories.png", as(png) replace	




* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------
* ANALYSIS USING DATA AT THE BARCODE LEVEL
* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------

set more off
use "$panel/Panel productos/Panel_BC_week_FINAL_1314.dta", clear

* Time variables created
g year=year(dofw(fecha_sem))
g month=month(dofw(fecha_sem))
g week=week(dofw(fecha_sem))
tab week, g(wk)

* Keeping relevant observations
* Years
keep if year>=2013 & year<=2014

* Running variable and treatment variable
g G=fecha_sem-wofd(mdy(1,1,2014))
g G2=G^2
g G3=G^3
g G4=G^4
g T=(G>=0)
g T_G=T*G
g T_G2=T*G2
g T_G3=T*G3
g T_G4=T*G4

* Normalize calories per unit
sum CpU_kcal, d
drop if CpU_kcal>r(p99)
sum CpU_kcal, d
ren CpU_kcal CpU
g CpUc=CpU/100
g T_CpUc=T*CpUc
g T_CpU=T*CpU
g T_CpU2=T*CpU^2
g CpU2=CpU^2
g T_CpU3=T*CpU^3
g CpU3=CpU^3

g precio_PC=gasto/units_compra
drop _merge

* Merge product characteristics and create groups
merge m:1 idartigo using "$panel\base_estandarizada_barcodes_FINAL.dta", ///
      keepus(producto subproducto clas01 clas02 clas03 clas04 clas05 marca vol_unitario) keep(mas mat)

egen ID_prod=group(producto subproducto clas01 clas02 clas03 clas04 clas05), missing

bysort fecha_sem: egen agg_cal=sum(kcal_tot)
gen perc_cal=kcal_tot/agg_cal
bysort fecha_sem: egen agg_gasto=sum(gasto)
gen perc_gto=gasto/agg_gasto

xtile decSD=CpU if year>=2013 & year<=2014 & tax_beverage==1 [pw=perc_gto], nq(10)
xtile decHCF=CpU if year>=2013 & year<=2014 & tax_food==1 [pw=perc_gto], nq(10)

* Partial out week dummies and product FE
g log_preco=log(preco)
areg log_preco wk2-wk52 if year>=2013 & year<=2014 & tax_beverage==1 [pw=perc_gto], abs(decSD)
predict preco_SD_res, resid
areg log_preco wk2-wk52 if year>=2013 & year<=2014 & tax_food==1 [pw=perc_gto], abs(decHCF)
predict preco_HCF_res, resid
g log_units_compra=log(units_compra+1)
areg log_units_compra wk2-wk52 if year>=2013 & year<=2014 & tax_beverage==1 [pw=perc_gto], abs(decSD)
predict units_compra_SD_res, resid
areg log_units_compra wk2-wk52 if year>=2013 & year<=2014 & tax_food==1 [pw=perc_gto], abs(decHCF)
predict units_compra_HCF_res, resid


* Histogram of calories per unit
*histogram CpU_norm if tax_beverage==1 
*histogram CpU_norm if tax_food==1 
*histogram CpU_norm if tax_beverage==1 & CpU_kcal<2000
*histogram CpU_norm if tax_food==1 & CpU_kcal<2000

* Estimation
preserve
sum CpU if tax_beverage==1, d 
local CpU_SD_min=`r(p5)'
local CpU_SD_max=`r(p90)'
*drop if CpU<=`r(p5)' | CpU>=`r(p95)'
sort CpU 

/* -----------------------------------------------------------------
                FIGURE 5: MECHANISMS (SD)
   ----------------------------------------------------------------- */
 
foreach var of varlist preco units_compra {

	reg `var'_SD_res T T_CpU T_CpU2 G G2 T_G T_G2 CpU CpU2 if year>=2013 & year<=2014 & tax_beverage==1 [pw=perc_gto], vce(cluster idartigo)
	gen d2`var'_SD=_b[T] + _b[T_CpU]*CpU + _b[T_CpU2]*CpU^2 if e(sample)==1
	sort CpU
	twoway (histogram CpU, fcolor(gray) yaxis(2) yscale(alt range(0 0.004) axis(2)) fintensity(30)) (line d2`var'_SD CpU, yscale(alt) yaxis(1)) ///
			if CpU>=`CpU_SD_min' & CpU<=`CpU_SD_max', ///
		ytitle(`var'_SD) xtitle(Calories per Unit) legend(off) graphregion(fcolor(white) ifcolor(white))
	graph export "$graphs\Paper\Figure5_MechanismSD_`var'.png", as(png) replace

	}
	
restore 


sum CpU if tax_food==1, d 
local CpU_HCF_min=`r(p5)'
local CpU_HCF_max=`r(p90)'
*drop if CpU<=`r(p5)' | CpU>=`r(p90)'
sort CpU 

/* -----------------------------------------------------------------
                FIGURE 5: MECHANISMS (HCF)
   ----------------------------------------------------------------- */

foreach var of varlist preco units_compra {

	reg `var'_HCF_res T T_CpU T_CpU2 G G2 T_G T_G2 CpU CpU2 if year>=2013 & year<=2014 & tax_food==1 [pw=perc_gto], vce(cluster idartigo)
	gen d2`var'_HCF=_b[T] + _b[T_CpU]*CpU + _b[T_CpU2]*CpU^2 if e(sample)==1
	twoway (histogram CpU, fcolor(gray) yaxis(2) yscale(alt range(0 0.004) axis(2)) fintensity(30)) (line d2`var'_HCF CpU, yscale(alt) yaxis(1)) ///
		if CpU>=`CpU_HCF_min' & CpU<=`CpU_HCF_max', ///
		ytitle(`var'_HCF) xtitle(Calories per Unit) legend(off) graphregion(fcolor(white) ifcolor(white))
	graph export "$graphs\Paper\Figure5_MechanismHCF_`var'.png", as(png) replace

	}

/* -----------------------------------------------------------------
                TABLE A.4 AND A.5: TAXED PRODUCTS
   ----------------------------------------------------------------- */

collapse tax_beverage tax_food tax_ID, by(producto)

* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------
* RD USING 275 AS CUTOFF   kcal_100gr
* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------

set more off
use "$panel/Panel productos/Panel_BC_week_FINAL_1314.dta", clear

global Y_RD275 "D_preco"

* Time variables created
g year=year(dofw(fecha_sem))
g month=month(dofw(fecha_sem))
g week=week(dofw(fecha_sem))
tab week, g(wk)

* Keeping relevant observations
keep if tax_food!=.
*keep if year==2014

* Merge to recover the running variable
drop _m
merge m:1 idartigo using "$panel\base_estandarizada_barcodes_FINAL.dta", ///
      keepus(producto subproducto clas01 clas02 clas03 clas04 clas05 marca vol_unitario kcal_100gr) keep(mas mat)

* Variables for RD analysis
g G=kcal_100gr-275
g G2=G^2
g G3=G^3
g T=(G>=0)
g T_G=T*G
g T_G2=T*G2
g T_G3=T*G3

* Weights
bysort fecha_sem: egen agg_gasto=sum(gasto)
gen perc_gto=gasto/agg_gasto

* Dummies by type of product
egen ID_prod=group(producto subproducto clas01 clas02 clas03 clas04 clas05), missing
xtset idartigo fecha_sem

* Partial out week dummies and product FE
g log_preco=log(preco)
g D_preco=log_preco-L52.log_preco
reg D_preco wk2-wk52 
predict D_preco_res, resid
g log_units_compra=log(units_compra+1)
g D_units=log_units_compra-L52.log_units_compra
reg D_units wk2-wk52 
predict D_units_res, resid
g log_calories=log(1+kcal_tot)
g D_calories=log_calories-L52.log_calories
reg D_calories wk2-wk52 
predict D_calories_res, resid

* Create conditional averages by outcome
* Here we define how many weeks to aggregate in the conditional means for the scatter
g G_r=floor(G/5)
preserve
collapse (mean) D_preco D_units D_calories tax_food [pweight=perc_gto], by(G_r)

foreach var of varlist  D_preco D_units D_calories tax_food {
	ren `var' `var'_MN
	}
save "$panel\Means275.dta", replace
restore

* Rescale it to avoid many dots in the scatterplot
drop _merge
merge m:1 G_r using "$panel\Means275.dta", keepus(*_MN) keep(mas mat)
replace G_r=G_r*5


* Estimation
* TABLE
reg tax_food T G G2 T_G T_G2 [pweight=perc_gto] if abs(G)<=100, vce(robust)
*rdrobust tax_food G, c(0) p(2) kernel(triangular) h(100) 
	outreg2 using "$tables\Paper\TableXL_RD275.txt", text ct(`var') replace

local contador=1
foreach var of varlist $Y_RD275 {
	reg `var'_res T G G2 T_G T_G2 [pweight=perc_gto] if abs(G)<=100, vce(robust)
*	rdrobust `var'_res G, c(0) p(2) kernel(triangular) h(100)
	if `contador'<3 {
	outreg2 using "$tables\Paper\TableXL_RD275.txt", text ct(`var') append
		}
	if `contador'==3 {
	outreg2 using "$tables\Paper\TableXL_RD275.tex", tex ct(`var') append
		}
	local contador=`contador'+1
	}


/* -----------------------------------------------------------------
                FIGURE A.6: RD 	OF FOOD USING THE 275 RULE
   ----------------------------------------------------------------- */

foreach var of varlist $Y_RD275 {
twoway (qfitci `var'_res G [pweight=perc_gto] if G<0 & G>=-100, clcolor(navy) clwidth(medium)) ///
		(qfitci `var'_res G [pweight=perc_gto] if G>=0 & G<=100, clcolor(navy) clwidth(medium)) ///
		(scatter `var'_MN G_r if abs(G)<100), ///
        ytitle(`var') ///
		xscale(range(-100 100)) xtitle(Calories per 100gr (with respect to 275)) ///
		xline(0, lwidth(medthin) lpattern(solid) lcolor(black)) legend(off) ///
		graphregion(fcolor(white) ifcolor(white))
	graph export "$graphs\Paper\FigureL_RD275_`var'.png", as(png) replace
	}

* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------
* ANALYSIS USING 2012-2014 DATA
* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------

set more off	
use "$panel/Panel hogares/Panel_HH_week_FINAL_1214.dta", clear

* Time variables created
g year=year(dofw(fecha_sem))
g month=month(dofw(fecha_sem))
tab week, g(wk)
tab month, g(mth)

* Filter to select households that appear at least once in both 2013 and 2014
g yy0=(year==2012)
g yy1=(year==2013)
g yy2=(year==2014)
bysort iddomicilio: egen y2012=max(yy0)
bysort iddomicilio: egen y2013=max(yy1)
bysort iddomicilio: egen y2014=max(yy2)

* Running variable and treatment variable
g G=fecha_sem-wofd(mdy(1,1,2014))
g T=(G>=0)
g G2=G^2

* Keeping relevant observations
* Years
keep if year==2012 | year==2013 | year==2014
* Households
keep if  y2013==1 & y2014==1

* Partial out week dummies and HH FE
foreach var of varlist $Y_prices {
	xtreg `var' wk2-wk52 if year==2013 | year==2014, fe
	predict `var'_res, resid
	}
foreach var of varlist $Y_calories $Y_nutrients {
	g log_`var'=log(`var'+1)
	xtreg log_`var' wk2-wk52 if year==2013 | year==2014, fe
	predict `var'_res, resid
	}
foreach var of varlist $Y_liters $Y_kilos {
	g log_`var'=log(`var'+1)
	xtreg log_`var' mth1-mth11 if year==2013 | year==2014, fe
	predict `var'_res, resid
	}

* Create weekly averages by outcome
* Here we define how many weeks to aggregate in the conditional means for the scatter
g G_r=floor(G/4)
foreach var of varlist $Y_prices {
	bysort G: egen `var'_MNP=mean(`var')
	}
foreach var of varlist $Y_calories {
	bysort G: egen `var'_MNP=mean(`var')
	}
foreach var of varlist $Y_prices $Y_calories $Y_nutrients {
	bysort G_r: egen `var'_MN=mean(`var'_res)
	}
* Rescale it to avoid many dots in the scatterplot
replace G_r=G_r*4
	
/*
local contador=1
foreach var of varlist $Y_prices $Y_calories {
	rdrobust `var'_res G, c(0) p(2) kernel(triangular) h(52)
	if `contador'==1 {
	outreg2 using "$tables\Paper\Table2B_Check.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes replace
		}
	if `contador'>1 & `contador'<10 {
	outreg2 using "$tables\Paper\Table2B_Check.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	if `contador'==10 {
	outreg2 using "$tables\Paper\Table2B_Check.tex", tex ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	local contador=`contador'+1
	}

/* -----------------------------------------------------------------
                TABLE A.7: SENSITIVITY TO BANDWIDTHS
   ----------------------------------------------------------------- */

local contador=1

foreach bw of numlist 26 39 52 65 78 {

foreach var of varlist $Y_prices $Y_calories {

	rdrobust `var'_res G, c(0) p(2) kernel(triangular) h(`bw')
	
	* Tabla comparacion residuales
	if `contador'==1 {
	outreg2 using "$tables\Paper\TableA2_BW`bw'.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes replace
		}
	if `contador'>1 & `contador'<10 {
	outreg2 using "$tables\Paper\TableA2_BW`bw'.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	if `contador'==10 {
	outreg2 using "$tables\Paper\TableA2_BW`bw'.tex", tex ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	local contador=`contador'+1
	}
	local contador=1
	}

/* -----------------------------------------------------------------
                TABLE A.11: INVENTORY BEHAVIOR
   ----------------------------------------------------------------- */
*/	
g I1=(month==12 & year==2012)
g I2=(month==12 & year==2013)
local contador=1

foreach var of varlist $Y_liters $Y_kilos {
	reg `var'_res I1 I2 G G2 if year>=2012 & year<=2014, vce(cluster iddomicilio)
	if `contador'==1 {
		outreg2 I1 I2 using "$tables\Paper\TableA11_Inventory.txt", text bdec(3) sdec(4) ct(`var') ///
			keep(I1 I2) nocons replace
		}
	if `contador'>1 & `contador'<4 {
		outreg2 I1 I2 using "$tables\Paper\TableA11_Inventory.txt", text bdec(3) sdec(4) ct(`var') ///
			keep(I1 I2) nocons append
		}
	if `contador'==4 {
		outreg2 I1 I2 using "$tables\Paper\TableA11_Inventory.tex", tex bdec(3) sdec(4) ct(`var') ///
			keep(I1 I2) nocons append
		}
	local contador=`contador'+1
	}
	
* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------
* PERMUTATION ANALYSIS
* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------
set more off
use "$panel/Panel hogares/Panel_HH_week_FINAL_1214.dta", clear

* Time variables created
g year=year(dofw(fecha_sem))
g month=month(dofw(fecha_sem))
tab week, g(wk)

* Filter to select households that appear at least once in both 2013 and 2014
g yy0=(year==2012)
g yy1=(year==2013)
g yy2=(year==2014)
bysort iddomicilio: egen y2012=max(yy0)
bysort iddomicilio: egen y2013=max(yy1)
bysort iddomicilio: egen y2014=max(yy2)

* Keeping relevant observations
* Years
keep if year==2012 | year==2013 | year==2014
* Households
keep if y2013==1 & y2014==1

* Partial out week dummies and HH FE
foreach var of varlist $Y_prices {
	xtreg `var' wk2-wk52 if year==2013 | year==2014, fe
	predict `var'_res, resid
	}
foreach var of varlist $Y_calories $Y_nutrients $Y_liters $Y_kilos {
	g log_`var'=log(`var'+1)
	xtreg log_`var' wk2-wk52 if year==2013 | year==2014, fe
	predict `var'_res, resid
	}

*** Loop to generate the estimation by changing the date in which the tax got implemented
* Prices
foreach var of varlist $Y_prices {
local contador=1
	forvalues i=-78/26 {

	g G=fecha_sem-wofd(mdy(1,1,2014))-`i'
	g T=(G>=0)

/* -----------------------------------------------------------------
                FIGURE A.7: FISCHER EXACT TEST INPUT
   ----------------------------------------------------------------- */
	rdrobust `var'_res G, c(0) p(2) kernel(triangular) h(52) vce(cluster iddomicilio)
	if `contador'==1 {
	outreg2 using "$tables\Table_PermutationRD_`var'.xls", excel bdec(3) sdec(4) noaster ct(`var') ///
		nocons replace
		}
	if `contador'>1 {
	outreg2 using "$tables\Table_PermutationRD_`var'.xls", excel bdec(3) sdec(4) noaster ct(`var') ///
		nocons append
		}
	local contador=`contador'+1
	drop G T 
	}	
}

* Calorie consumption
foreach var of varlist $Y_calories {
local contador=1
	forvalues i=-78/26 {

	g G=fecha_sem-wofd(mdy(1,1,2014))-`i'
	g T=(G>=0)

	rdrobust `var'_res G, c(0) p(2) kernel(triangular) h(52) vce(cluster iddomicilio)
	if `contador'==1 {
	outreg2 using "$tables\Table_PermutationRD_`var'.xls", excel bdec(3) sdec(4) noaster ct(`var') ///
		nocons replace
		}
	if `contador'>1 {
	outreg2 using "$tables\Table_PermutationRD_`var'.xls", excel bdec(3) sdec(4) noaster ct(`var') ///
		nocons append
		}
	local contador=`contador'+1
	drop G T 
	}	
}
* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------
* SENSITIVITY ANALYSIS OF CHANGING THE TAX DEFINITION FOR PRODUCTS IN WHICH
* OUR CLASSIFICATION AND DELOITTE'S DOES NOT MATCH
* ---------------------------------------------------------------------------------
* ---------------------------------------------------------------------------------
set more off
use "$panel/Panel hogares/Panel_HH_week_FINAL_Sens_1314.dta", clear

* Variables de tiempo
g year=year(dofw(fecha_sem))
g month=month(dofw(fecha_sem))
*g week=week(dofw(fecha_sem))
tab week, g(wk)

* Filter to select households that appear at least once in both 2013 and 2014
g yy0=(year==2012)
g yy1=(year==2013)
g yy2=(year==2014)
bysort iddomicilio: egen y2012=max(yy0)
bysort iddomicilio: egen y2013=max(yy1)
bysort iddomicilio: egen y2014=max(yy2)

* Running variable and treatment variable
g G=fecha_sem-wofd(mdy(1,1,2014))
g G2=G^2
g G3=G^3
g G4=G^4
g T=(G>=0)
g T_G=T*G
g T_G2=T*G2
g T_G3=T*G3
g T_G4=T*G4


* Keeping relevant observations
* Years
keep if year==2013 | year==2014
* Households
keep if y2013==1 & y2014==1

* Partial out week dummies and HH FE
foreach var of varlist $Y_prices {
	xtreg `var' wk2-wk52, fe
	predict `var'_res, resid
	}
foreach var of varlist $Y_calories {
	g log_`var'=log(`var'+1)
	xtreg log_`var' wk2-wk52, fe
	predict `var'_res, resid
	}

* Create weekly averages by outcome
* Here we define how many weeks to aggregate in the conditional means for the scatter
g G_r=floor(G/4)
foreach var of varlist $Y_prices {
	bysort G: egen `var'_MNP=mean(`var')
	}
foreach var of varlist $Y_prices {
	bysort G_r: egen `var'_MN=mean(`var'_res)
	}
* Rescale it to avoid many dots in the scatterplot
replace G_r=G_r*4

reg kcal_tot_res T G
g muestra=(e(sample)==1)

/* -----------------------------------------------------------------
                TABLE A.9: SENSITIVITY TO TAX DEFINITION
   ----------------------------------------------------------------- */
local contador=1
foreach var of varlist $Y_prices $Y_calories {
	rdrobust `var'_res G, c(0) p(2) kernel(triangular) h(52)
	if `contador'==1 {
	outreg2 using "$tables\Paper\TableA9_RD_robust.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes replace
		}
	if `contador'>1 & `contador'<10 {
	outreg2 using "$tables\Paper\TableA9_RD_robust.txt", text ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	if `contador'==10 {
	outreg2 using "$tables\Paper\TableA9_RD_robust.tex", tex ct(`var') ///
		addstat(Bias corrected effect, e(tau_bc), Robust Std. Error, e(se_tau_rb), ///
		Polynomial order ,e(p), Bandwidth ,e(h_l), Num obs, `e(N_h_l)'+`e(N_h_r)') noobs nonotes append
		}
	local contador=`contador'+1
	}