**----------------------------------------------**
		* Started on: December 28, 2019 *
		* Last updated on: July 9, 2023 *
**----------------------------------------------**
clear
clear matrix
mata: mata clear
capture log close
set more off

**Making global directories**

global data "~\GACC Study\Kenya Survey\Data"

global dofiles "~\GACC Study\Kenya Survey\Dofiles"

global log "~\GACC Study\Kenya Survey\Logfiles"

global results "~\GACC Study\Kenya Survey\Results"


ssc install mfx2
ssc install margeff


use "$data\Kenya_rosterdata_anon_9July2023.dta", clear



**************************************************************************************************************
**        1. MLOGIT of fuel choice (Charcoal, LPG, Kerosene)  (w/ log per capita monthly expenditures)      **
**************************************************************************************************************
			 
mlogit primary_fuel lpc_monthlyexp_usd hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust) baseoutcome (0)
margins, dydx(*) post predict (outcome(Kerosene))
outreg2 using "$results\mlogit_pcmoexp_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace

mlogit primary_fuel lpc_monthlyexp_usd hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust) baseoutcome (0)
margins, dydx(*) post predict (outcome(LPG))
outreg2 using "$results\mlogit_pcmoexp_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append


**************************************************************************************************************
**        1A. MLOGIT of fuel choice (Charcoal, LPG, Kerosene)  (w/ wealth index)   							**
**************************************************************************************************************

mlogit primary_fuel Wealthindex hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust) baseoutcome (0)
margins, dydx(*) post predict (outcome(Kerosene))
outreg2 using "$results\mlogit_pcmoexp_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

mlogit primary_fuel Wealthindex hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust) baseoutcome (0)
margins, dydx(*) post predict (outcome(LPG))
outreg2 using "$results\mlogit_pcmoexp_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append



**************************************************************************
**    2. Stacking of fuels (w/ log per capita monthly expenditures) ** 
**************************************************************************

reg stackfuels_num lpc_monthlyexp_usd hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust)
outreg2 using "$results\stack_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace

reg stackindex lpc_monthlyexp_usd hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust)
outreg2 using "$results\stack_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

probit binary_stack lpc_monthlyexp_usd hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust)
margins, dydx(_all) post
outreg2 using "$results\stack_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg charstackperc_new lpc_monthlyexp_usd hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust)
outreg2 using "$results\stack_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg kerostackperc_new lpc_monthlyexp_usd hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust)
outreg2 using "$results\stack_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg LPGstackperc_new lpc_monthlyexp_usd hhsize head_age head_female i.head_education i.cook_dmaker_educ_HH depratio private_latrine savings_account g_8 a_3, vce(robust)
outreg2 using "$results\stack_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append



***************************
** 	3A. WTP for CHARCOAL **
***************************

probit resp_charcoal log_firstprice_charcoal if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_charcoal_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace

probit resp_charcoal log_firstprice_charcoal lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 char_weekly_colltime binary_stack if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_charcoal_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

probit resp_charcoal log_firstprice_charcoal lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 char_weekly_colltime charstackperc_new if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_charcoal_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append


**************************
** 3B. WTP for KEROSENE **
**************************

probit resp_kerosene log_firstprice_kerosene if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_kerosene_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace

probit resp_kerosene log_firstprice_kerosene lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 kero_weekly_colltime binary_stack if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_kerosene_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

probit resp_kerosene log_firstprice_kerosene lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 kero_weekly_colltime kerostackperc_new if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_kerosene_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append
 

*********************
** 3C. WTP for LPG **
*********************

probit resp_lpg log_firstprice_lpg if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_lpg_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace

probit resp_lpg log_firstprice_lpg lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 lpg_weekly_colltime binary_stack if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_lpg_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

probit resp_lpg log_firstprice_lpg lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 lpg_weekly_colltime LPGstackperc_new if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_lpg_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append


********************************************
** 3D. WTP for Cooking Fuels POOLED MODEL **
********************************************

probit response_cv log_firstprice_all i.fuel_type if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_pooled_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace 

probit response_cv log_firstprice_all i.fuel_type lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 fuel_coll_weekly binary_stack if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_pooled_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

probit response_cv log_firstprice_all i.fuel_type lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 fuel_coll_weekly charstackperc_new kerostackperc_new LPGstackperc_new if a_1==1, ro
margins, dydx(_all) post
outreg2 using "$results\Kenya_pooled_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append



*******************************************************************
** 4A. Change in percentage cooking for non-switchers (Charcoal) **
*******************************************************************

reg perc_cooking_charcoal log_firstprice_charcoal if (a_1==1 & resp_charcoal==1), ro
outreg2 using "$results\Kenya_charcoal_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace

reg perc_cooking_charcoal log_firstprice_charcoal lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_c if (a_1==1 & resp_charcoal==1), ro
outreg2 using "$results\Kenya_charcoal_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg perc_cooking_charcoal log_firstprice_charcoal lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_c binary_stack if (a_1==1 & resp_charcoal==1), ro
outreg2 using "$results\Kenya_charcoal_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg perc_cooking_charcoal log_firstprice_charcoal lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_c charstackperc_new if (a_1==1 & resp_charcoal==1), ro
outreg2 using "$results\Kenya_charcoal_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append


*******************************************************************
** 4B. Change in percentage cooking for non-switchers (Kerosene) **
*******************************************************************

reg perc_cooking_kerosene log_firstprice_kerosene if (a_1==1 & resp_kerosene==1), ro
outreg2 using "$results\Kenya_kerosene_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace

reg perc_cooking_kerosene log_firstprice_kerosene lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_d if (a_1==1 & resp_kerosene==1), ro
outreg2 using "$results\Kenya_kerosene_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg perc_cooking_kerosene log_firstprice_kerosene lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_d binary_stack if (a_1==1 & resp_kerosene==1), ro
outreg2 using "$results\Kenya_kerosene_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg perc_cooking_kerosene log_firstprice_kerosene lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_d kerostackperc_new if (a_1==1 & resp_kerosene==1), ro
outreg2 using "$results\Kenya_kerosene_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append


*******************************************************************
** 4C. Change in percentage cooking for non-switchers (LPG) **
*******************************************************************

reg perc_cooking_lpg log_firstprice_lpg if (a_1==1 & resp_lpg==1), ro
outreg2 using "$results\Kenya_lpg_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace

reg perc_cooking_lpg log_firstprice_lpg lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_a if (a_1==1 & resp_lpg==1), ro
outreg2 using "$results\Kenya_lpg_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg perc_cooking_lpg log_firstprice_lpg lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_a binary_stack if (a_1==1 & resp_lpg==1), ro
outreg2 using "$results\Kenya_lpg_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg perc_cooking_lpg log_firstprice_lpg lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 e_5_a LPGstackperc_new if (a_1==1 & resp_lpg==1), ro
outreg2 using "$results\Kenya_lpg_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append


*******************************************************************
** 4D. Change in percentage cooking for non-switchers (Pooled) **
*******************************************************************

reg perc_cooking_all log_firstprice_all i.fuel_type if (a_1==1 & response_cv==1), ro
outreg2 using "$results\Kenya_pooled_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label replace 

reg perc_cooking_all log_firstprice_all i.fuel_type lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 fuel_coll_weekly binary_stack if (a_1==1 & response_cv==1), ro
outreg2 using "$results\Kenya_pooled_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append

reg perc_cooking_all log_firstprice_all i.fuel_type lpc_monthlyexp_usd hhsize head_age head_female i.head_education depratio savings_account g_8 fuel_coll_weekly charstackperc_new kerostackperc_new LPGstackperc_new if (a_1==1 & response_cv==1), ro
outreg2 using "$results\Kenya_pooled_cooking_5April22.xls", se dec (3) stats (coef se) e(r2_p) label append



*******************************************************************************
		* 5. Contingent Valuation using Double-Bounded Dichotomous Choice *
*******************************************************************************

ssc install doubleb
* https://ideas.repec.org/c/boc/bocode/s457168.html *

log using "$log/Kenya_WTP Analysis_doubleb_5April22.log", replace

*Creating scalar values for calculating WTP* 

sum lpc_monthlyexp_usd, meanonly
scalar lpc_monthlyexp_usd_m = r(mean)

sum hhsize, meanonly
scalar hhsize_m = r(mean)

sum head_age, meanonly
scalar head_age_m = r(mean)

sum head_female, meanonly
scalar head_female_m = r(mean)

sum depratio, meanonly
scalar depratio_m = r(mean)

sum savings_account, meanonly
scalar savings_account_m = r(mean)

sum g_8, meanonly
scalar g_8_m = r(mean)

sum char_weekly_colltime, meanonly
scalar char_weekly_colltime_m = r(mean)

sum kero_weekly_colltime, meanonly
scalar kero_weekly_colltime_m = r(mean)

sum lpg_weekly_colltime, meanonly
scalar lpg_weekly_colltime_m = r(mean)

sum fuel_coll_weekly, meanonly
scalar fuel_coll_weekly_m = r(mean)

sum binary_stack, meanonly
scalar binary_stack_m = r(mean)

sum charstackperc_new, meanonly
scalar charstackperc_new_m = r(mean)

sum kerostackperc_new, meanonly
scalar kerostackperc_new_m = r(mean)

sum LPGstackperc_new, meanonly
scalar LPGstackperc_new_m = r(mean)


******************************************************************
**	5A. CALCULATING WILLINGNESS TO PAY (Double bound): CHARCOAL **	
******************************************************************

doubleb log_firstprice_charcoal log_secondprice_charcoal resp_charcoal resp_charcoal_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 char_weekly_colltime if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+char_weekly_colltime_m*_b[char_weekly_colltime])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpchar=xbeta 
sum wtpchar
sum wtpchar if resp_charcoal!=.
*/

doubleb log_firstprice_charcoal log_secondprice_charcoal resp_charcoal resp_charcoal_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 char_weekly_colltime binary_stack if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+char_weekly_colltime_m*_b[char_weekly_colltime]+binary_stack_m*_b[binary_stack])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpchar_stack=xbeta 
sum wtpchar_stack
sum wtpchar_stack if resp_charcoal!=.
*/

doubleb log_firstprice_charcoal log_secondprice_charcoal resp_charcoal resp_charcoal_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 char_weekly_colltime charstackperc_new if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+char_weekly_colltime_m*_b[char_weekly_colltime]+charstackperc_new_m*_b[charstackperc_new])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpchar_stack2=xbeta 
sum wtpchar_stack2
sum wtpchar_stack2 if resp_charcoal!=.

sum wtpchar wtpchar_stack wtpchar_stack2 
sum wtpchar wtpchar_stack wtpchar_stack2 if resp_charcoal!=.
*/


******************************************************************
**	5B. CALCULATING WILLINGNESS TO PAY (Double bound): KEROSENE **	
******************************************************************

doubleb log_firstprice_kerosene log_secondprice_kerosene resp_kerosene resp_kerosene_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 kero_weekly_colltime if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+kero_weekly_colltime_m*_b[kero_weekly_colltime])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpkero=xbeta 
sum wtpkero
sum wtpkero if resp_kerosene!=.
*/

doubleb log_firstprice_kerosene log_secondprice_kerosene resp_kerosene resp_kerosene_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 kero_weekly_colltime binary_stack if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+kero_weekly_colltime_m*_b[kero_weekly_colltime]+binary_stack_m*_b[binary_stack])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpkero_stack=xbeta 
sum wtpkero_stack
sum wtpkero_stack if resp_kerosene!=.
*/

doubleb log_firstprice_kerosene log_secondprice_kerosene resp_kerosene resp_kerosene_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 kero_weekly_colltime kerostackperc_new if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+kero_weekly_colltime_m*_b[kero_weekly_colltime]+kerostackperc_new_m*_b[kerostackperc_new])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpkero_stack2=xbeta 
sum wtpkero_stack2
sum wtpkero_stack2 if resp_kerosene!=.

sum wtpkero wtpkero_stack wtpkero_stack2 
sum wtpkero wtpkero_stack wtpkero_stack2 if resp_kerosene!=.
*/


*****************************************************************
**	5C. CALCULATING WILLINGNESS TO PAY (Double bound): LPG **	
*****************************************************************	

doubleb log_firstprice_lpg log_secondprice_lpg resp_lpg resp_lpg_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 lpg_weekly_colltime if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+lpg_weekly_colltime_m*_b[lpg_weekly_colltime])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtplpg=xbeta 
sum wtplpg
sum wtplpg if resp_lpg!=.
*/

doubleb log_firstprice_lpg log_secondprice_lpg resp_lpg resp_lpg_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 lpg_weekly_colltime binary_stack if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+lpg_weekly_colltime_m*_b[lpg_weekly_colltime]+binary_stack_m*_b[binary_stack])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtplpg_stack=xbeta 
sum wtplpg_stack
sum wtplpg_stack if resp_lpg!=.
*/

doubleb log_firstprice_lpg log_secondprice_lpg resp_lpg resp_lpg_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 lpg_weekly_colltime LPGstackperc_new if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+lpg_weekly_colltime_m*_b[lpg_weekly_colltime]+LPGstackperc_new_m*_b[LPGstackperc_new])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtplpg_stack2=xbeta 
sum wtplpg_stack2
sum wtplpg_stack2 if resp_lpg!=.

sum wtplpg wtplpg_stack wtplpg_stack2 
sum wtplpg wtplpg_stack wtplpg_stack2 if resp_lpg!=.
*/


*****************************************************************
**	5D. CALCULATING WILLINGNESS TO PAY (Double bound): Pooled **	
*****************************************************************

doubleb log_firstprice_all log_secondprice_all response_cv second_response_cv lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 fuel_coll_weekly if a_1==1 
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+fuel_coll_weekly_m*_b[fuel_coll_weekly])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtp_pooled=xbeta 
sum wtp_pooled
sum wtp_pooled if response_cv!=.
*/

doubleb log_firstprice_all log_secondprice_all response_cv second_response_cv lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 fuel_coll_weekly binary_stack if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+fuel_coll_weekly_m*_b[fuel_coll_weekly]+binary_stack_m*_b[binary_stack])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtp_pooled_stack=xbeta 
sum wtp_pooled_stack
sum wtp_pooled_stack if response_cv!=.

sum wtp_pooled wtp_pooled_stack 
sum wtp_pooled wtp_pooled_stack  if response_cv!=.
*/

doubleb log_firstprice_all log_secondprice_all response_cv second_response_cv lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 fuel_coll_weekly charstackperc_new kerostackperc_new LPGstackperc_new if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+fuel_coll_weekly_m*_b[fuel_coll_weekly]+charstackperc_new_m*_b[charstackperc_new]+kerostackperc_new_m*_b[kerostackperc_new]+LPGstackperc_new_m*_b[LPGstackperc_new])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtp_pooled_stack2=xbeta 
sum wtp_pooled_stack2
sum wtp_pooled_stack2 if response_cv!=.

sum wtp_pooled wtp_pooled_stack wtp_pooled_stack2
sum wtp_pooled wtp_pooled_stack wtp_pooled_stack2 if response_cv!=.
*/

log close 




log using "$log/Kenya_WTP Analysis_doubleb_NOLOG_5April22.log", replace

*Creating scalar values for calculating WTP* 

sum lpc_monthlyexp_usd, meanonly
scalar lpc_monthlyexp_usd_m = r(mean)

sum hhsize, meanonly
scalar hhsize_m = r(mean)

sum head_age, meanonly
scalar head_age_m = r(mean)

sum head_female, meanonly
scalar head_female_m = r(mean)

sum depratio, meanonly
scalar depratio_m = r(mean)

sum savings_account, meanonly
scalar savings_account_m = r(mean)

sum g_8, meanonly
scalar g_8_m = r(mean)

sum char_weekly_colltime, meanonly
scalar char_weekly_colltime_m = r(mean)

sum kero_weekly_colltime, meanonly
scalar kero_weekly_colltime_m = r(mean)

sum lpg_weekly_colltime, meanonly
scalar lpg_weekly_colltime_m = r(mean)

sum fuel_coll_weekly, meanonly
scalar fuel_coll_weekly_m = r(mean)

sum binary_stack, meanonly
scalar binary_stack_m = r(mean)

sum charstackperc_new, meanonly
scalar charstackperc_new_m = r(mean)

sum kerostackperc_new, meanonly
scalar kerostackperc_new_m = r(mean)

sum LPGstackperc_new, meanonly
scalar LPGstackperc_new_m = r(mean)


******************************************************************
**	5A. CALCULATING WILLINGNESS TO PAY (Double bound): CHARCOAL **	
******************************************************************

doubleb firstprice_charcoal_usd secondprice_charcoal_usd resp_charcoal resp_charcoal_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 char_weekly_colltime if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+char_weekly_colltime_m*_b[char_weekly_colltime])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpchar=xbeta 
sum wtpchar
sum wtpchar if resp_charcoal!=.
*/

doubleb firstprice_charcoal_usd secondprice_charcoal_usd resp_charcoal resp_charcoal_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 char_weekly_colltime binary_stack if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+char_weekly_colltime_m*_b[char_weekly_colltime]+binary_stack_m*_b[binary_stack])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpchar_stack=xbeta 
sum wtpchar_stack
sum wtpchar_stack if resp_charcoal!=.
*/

doubleb firstprice_charcoal_usd secondprice_charcoal_usd resp_charcoal resp_charcoal_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 char_weekly_colltime charstackperc_new if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+char_weekly_colltime_m*_b[char_weekly_colltime]+charstackperc_new_m*_b[charstackperc_new])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpchar_stack2=xbeta 
sum wtpchar_stack2
sum wtpchar_stack2 if resp_charcoal!=.

sum wtpchar wtpchar_stack wtpchar_stack2 
sum wtpchar wtpchar_stack wtpchar_stack2 if resp_charcoal!=.
*/


******************************************************************
**	5B. CALCULATING WILLINGNESS TO PAY (Double bound): KEROSENE **	
******************************************************************

doubleb firstprice_kerosene_usd secondprice_kerosene_usd resp_kerosene resp_kerosene_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 kero_weekly_colltime if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+kero_weekly_colltime_m*_b[kero_weekly_colltime])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpkero=xbeta 
sum wtpkero
sum wtpkero if resp_kerosene!=.
*/

doubleb firstprice_kerosene_usd secondprice_kerosene_usd resp_kerosene resp_kerosene_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 kero_weekly_colltime binary_stack if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+kero_weekly_colltime_m*_b[kero_weekly_colltime]+binary_stack_m*_b[binary_stack])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpkero_stack=xbeta 
sum wtpkero_stack
sum wtpkero_stack if resp_kerosene!=.
*/

doubleb firstprice_kerosene_usd secondprice_kerosene_usd resp_kerosene resp_kerosene_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 kero_weekly_colltime kerostackperc_new if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+kero_weekly_colltime_m*_b[kero_weekly_colltime]+kerostackperc_new_m*_b[kerostackperc_new])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtpkero_stack2=xbeta 
sum wtpkero_stack2
sum wtpkero_stack2 if resp_kerosene!=.

sum wtpkero wtpkero_stack wtpkero_stack2 
sum wtpkero wtpkero_stack wtpkero_stack2 if resp_kerosene!=.
*/


*****************************************************************
**	5C. CALCULATING WILLINGNESS TO PAY (Double bound): LPG **	
*****************************************************************	

doubleb firstprice_lpg_usd secondprice_lpg_usd resp_lpg resp_lpg_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 lpg_weekly_colltime if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+lpg_weekly_colltime_m*_b[lpg_weekly_colltime])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtplpg=xbeta 
sum wtplpg
sum wtplpg if resp_lpg!=.
*/

doubleb firstprice_lpg_usd secondprice_lpg_usd resp_lpg resp_lpg_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 lpg_weekly_colltime binary_stack if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+lpg_weekly_colltime_m*_b[lpg_weekly_colltime]+binary_stack_m*_b[binary_stack])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtplpg_stack=xbeta 
sum wtplpg_stack
sum wtplpg_stack if resp_lpg!=.
*/

doubleb firstprice_lpg_usd secondprice_lpg_usd resp_lpg resp_lpg_2 lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 lpg_weekly_colltime LPGstackperc_new if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+lpg_weekly_colltime_m*_b[lpg_weekly_colltime]+LPGstackperc_new_m*_b[LPGstackperc_new])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtplpg_stack2=xbeta 
sum wtplpg_stack2
sum wtplpg_stack2 if resp_lpg!=.

sum wtplpg wtplpg_stack wtplpg_stack2 
sum wtplpg wtplpg_stack wtplpg_stack2 if resp_lpg!=.
*/


*****************************************************************
**	5D. CALCULATING WILLINGNESS TO PAY (Double bound): Pooled **	
*****************************************************************

doubleb first_price_all_rescale second_price_all_rescale response_cv second_response_cv lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 fuel_coll_weekly if a_1==1 
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+fuel_coll_weekly_m*_b[fuel_coll_weekly])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtp_pooled=xbeta 
sum wtp_pooled
sum wtp_pooled if response_cv!=.
*/

doubleb first_price_all_rescale second_price_all_rescale response_cv second_response_cv lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 fuel_coll_weekly binary_stack if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+fuel_coll_weekly_m*_b[fuel_coll_weekly]+binary_stack_m*_b[binary_stack])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtp_pooled_stack=xbeta 
sum wtp_pooled_stack
sum wtp_pooled_stack if response_cv!=.

sum wtp_pooled wtp_pooled_stack 
sum wtp_pooled wtp_pooled_stack  if response_cv!=.
*/

doubleb first_price_all_rescale second_price_all_rescale response_cv second_response_cv lpc_monthlyexp_usd hhsize head_age head_female depratio savings_account g_8 fuel_coll_weekly charstackperc_new kerostackperc_new LPGstackperc_new if a_1==1
nlcom (WTP: (_b[_cons]+lpc_monthlyexp_usd_m*_b[lpc_monthlyexp_usd]+hhsize_m*_b[hhsize]+head_age_m*_b[head_age]+head_female_m*_b[head_female]+depratio_m*_b[depratio]+savings_account_m*_b[savings_account]+g_8_m*_b[g_8]+fuel_coll_weekly_m*_b[fuel_coll_weekly]+charstackperc_new_m*_b[charstackperc_new]+kerostackperc_new_m*_b[kerostackperc_new]+LPGstackperc_new_m*_b[LPGstackperc_new])), noheader
/*
matrix b=get(_b) 
matrix xbeta = b[.,1..12] 
matrix score wtp_pooled_stack2=xbeta 
sum wtp_pooled_stack2
sum wtp_pooled_stack2 if response_cv!=.

sum wtp_pooled wtp_pooled_stack wtp_pooled_stack2
sum wtp_pooled wtp_pooled_stack wtp_pooled_stack2 if response_cv!=.
*/

log close 



*-------------------------------------------------------------------------------
*Descriptive statistics - what fuels do people switch to?

*First we restrict the sample to kerosene users
preserve 

keep if primary_fuel==1

*Tabulate the responses of kerosene users to initial bid price increases
tab f_3b

*Categorize initial bid price increases into low and high price increases.
*We say low = 25-50% increases and high=100-200% increases.
 gen high_low =.
 replace high_low = 0 if percentage_label_a == 25 | percentage_label_a == 50
 replace high_low = 1 if percentage_label_a == 100 | percentage_label_a == 200
 label define hl 0"25-50% Increase" 1"100-200% Increase"
 label values high_low hl
 keep if high_low!=.
 
 *Divide the income distribution (OF kerosene users) into 2
 xtile bottomtop=lpc_monthlyexp, n(2) 

*What do Low income households do when faced with low price increases vs high price increases
tab f_3b if high_low==0 & bottomtop==1
 
tab f_3b if high_low==1 & bottomtop==1

*What do High income households do when faced with low price increases vs high price increases

tab f_3b if high_low==0 & bottomtop==2
 
tab f_3b if high_low==1 & bottomtop==2


*Plot the data. In the plot editing - I manually hide the 0 category which is required by slideplot
help plotplainblind

set scheme plotplain

set scheme plotplainblind

slideplot bar f_3b if bottomtop==1, by(high_low) pos(1 2 3 4 5) neg(0) percent ytitle("") title("Low Income Households")
slideplot bar f_3b if bottomtop==2, by(high_low) pos(1 2 3 4 5) neg(0) percent ytitle("") title("High Income Households")