clear
set mem 5g

*******************************************************************************
*Define globals for directories and files
*******************************************************************************
global path_main          	"H:\RCOptim\Results\~BLP delivery 16 CK\~BLP delivery 16 CK"
global path_temp          	"e:\crknittel"
global path_in            	"$path_main"
global path_out	  	  	"$path_main"

global file_in_optim            "$path_in\Optimization results\Optimization results.dta"
global file_optim_tmp		  "$path_temp\blp_optim_results_tmp.dta"

global file_in_mkt_power        "$path_in\market power results\blp_mkt_power_results_all.dta"
global file_in_mkt_power_tmp    "$path_temp\blp_mkt_power_results_all_tmp.dta"

global file_in_merger           "$path_in\merger results\blp_merger_results_all.dta"
global file_in_merger_tmp       "$path_temp\blp_merger_results_all_tmp.dta"

global file_BLP_analysis        "$path_temp\BLP_analysis.dta"
global file_BLP_analysis_final  "$path_temp\BLP_analysis_final.dta"

global path_figures              "$path_main\figures\"

********************************************************************************
*******************************************************************************
*Load market power results - BLP
*******************************************************************************
use "$file_in_mkt_power", clear
capture rename brand product
sort optmethod stvalue market product
save "$file_in_mkt_power_tmp", replace

*******************************************************************************
*Merge the optimization, convergence, market power and merger results - BLP
*Prepare a dataset to produce figures and tables
*******************************************************************************

use "$file_in_optim", clear

*******************************************************************************
*Convergence 
*******************************************************************************
capture gen conv2=0
replace conv2=1 if optmethod==1     & exitinfo==1
replace conv2=1 if optmethod==1     & exitinfo==2
replace conv2=1 if optmethod==2     & exitinfo==1
replace conv2=1 if optmethod==3     & exitinfo<4000
replace conv2=1 if optmethod==4     & exitinfo==2
replace conv2=1 if optmethod==5     & exitinfo==2
replace conv2=1 if optmethod==5     & exitinfo==4
replace conv2=1 if optmethod==7  	& fcnevals<4000
replace conv2=1 if optmethod==8     & exitinfo==1
replace conv2=1 if optmethod==9     & exitinfo==1
replace conv2=1 if optmethod==10    & exitinfo==1	& fcnevals==1020

save "$file_optim_tmp", replace
*******************************************************************************

sort  optmethod stvalue

merge optmethod stvalue using "$file_in_mkt_power_tmp"
tab   _merge
keep if _merge==3
drop _merge

sort  optmethod stvalue market product 
merge optmethod stvalue market product using "$file_in_merger_tmp"
tab _merge
keep if _merge==3
drop _merge
compress
save "$file_BLP_analysis", replace

use "$file_BLP_analysis", clear

gen     flag=0

*******************************************************************************
*Identify products with market shares
*corresponding to the four quartiles of the market share distribution
*******************************************************************************
sort optmethod stvalue share_obs market product, stable

global  market_25=market  in 554
global product_25=product in 554

global  market_50=market  in 1109
global product_50=product in 1109

global  market_75=market  in 1663
global product_75=product in 1663

global  market_100=market  in 2217
global product_100=product in 2217

gen flag_025  = (market==$market_25) *(product==$product_25)
gen flag_050  = (market==$market_50) *(product==$product_50)
gen flag_075  = (market==$market_75) *(product==$product_75)
gen flag_100  = (market==$market_100)*(product==$product_100)

table market product if flag_025 ==1, c(min share_obs) f(%8.4fc)
table market product if flag_050 ==1, c(min share_obs) f(%8.4fc)
table market product if flag_075 ==1, c(min share_obs) f(%8.4fc)
table market product if flag_100 ==1, c(min share_obs) f(%8.4fc)

*******************************************************************************
*Identify the results that correspond to the minimum value
*of the objective function for each of the 10 optimization routines
*******************************************************************************
gsort optmethod fval stvalue -conv2

egen fval_rank=group(optmethod fval stvalue conv2)
egen fval_rank_min=min(fval_rank), by(optmethod conv2)
replace fval_rank=fval_rank-fval_rank_min+1

gen      fval_optmethod_min_flag = 0
replace  fval_optmethod_min_flag = 1 if fval_rank==1&conv2==1

sum      fval if fval_optmethod_min_flag

egen hessians_min=rowmin(hessians_*)

cd "$path_figures"

*******************************************************************************
*Figure 24 - BLP
*******************************************************************************
global upp=0
qui sum elast_own, detail, if flag_025==1&conv2==1
qui return list
global low_025=-5

qui sum elast_own, detail, if flag_050==1&conv2==1
qui return list
global low_050=-6

qui sum elast_own, detail, if flag_075==1&conv2==1
qui return list
global low_075=-10.3

qui sum elast_own, detail, if flag_100==1&conv2==1
qui return list
global low_100=-10

global w_025=abs($low_025)/25
global w_025=abs($low_050)/25
global w_025=abs($low_075)/25
global w_025=abs($low_100)/25

#delimit ;
histogram elast_own if flag_025==1 & conv2==1 & elast_own>=$low_025  & elast_own<$upp, 
	w($w_025)
	title(25th Percentile) 
	xtitle(Own-Price Elasticity) 
	saving(BLP_25, replace) 
	xlabel(,format(%9.3f))
	freq;

histogram elast_own if flag_050==1 & conv2==1 & elast_own>=$low_050  & elast_own<$upp, 
	w($w_025)
	title(50th Percentile) 
	xtitle(Own-Price Elasticity) 
	saving(BLP_median, replace) 
	xlabel(,format(%9.3f))
	freq;
	
histogram elast_own if flag_075==1 & conv2==1 & elast_own>=$low_075  & elast_own<$upp, 
	w($w_025)
	title(75th Percentile) 
	xtitle(Own-Price Elasticity) 
	saving(BLP_75, replace) 
	xlabel(,format(%9.3f))
	freq;

histogram elast_own if flag_100==1 & conv2==1 & elast_own>=$low_100  & elast_own<$upp, 
	w($w_025)
	title(Maximum) 
	xtitle(Own-Price Elasticity) 
	saving(BLP_largest, replace) 
	xlabel(,format(%9.3f))
	freq;

graph combine BLP_25.gph BLP_median.gph BLP_75.gph BLP_largest.gph;
graph export figure24.wmf,replace;
#delimit cr

erase BLP_25.gph
erase BLP_median.gph 
erase BLP_75.gph
erase BLP_largest.gph

egen own_stddev=sd(elast_own), by(product market conv2)
sum own_stddev elast_own, detail, if conv2==1

egen own_stddev_best=sd(elast_own), by(product market conv2  fval_optmethod_min_flag)
sum own_stddev_best elast_own, detail, if conv2==1& fval_optmethod_min_flag==1

*******************************************************************************
*Figure 25 - BLP
*******************************************************************************
qui sum fval, detail, if conv2==1
qui return list
global fval_upp=r(p90)
#delimit ;
histogram fval if conv2==1 &fval<$fval_upp&product==1&market==16, 
	freq 
	w(5)
	xtitle(GMM Objective Values); 
graph export Figure25.wmf ,replace;