capture log close
clear
log using elasticity_bounds_growing_linear_demand, replace

********************************************************************************
*
* Program ELASTICITY_BOUNDS_GROWING_LINEAR_DEMAND.DO
*
* Version of main program used for estimating elasticity bounds in a growing
* market using methods from Mullin & Snyder AEJ Micro paper. This version
* is meant to be generic, working for a researcher's simple dataset that has
* already been cleaned. In addition, the program does not analyze several 
* functional forms but just linear demand. This allows for simpler graphing
* commands. The need to cram linear and logit together in one figure in the 
* original paper required "tricking" Stata into having several sets of paired
* columns with whiskers using some non-transparent commands. 
*
* Instructions on use: User needs to scroll to "READ IN DATA" section of
* program and provide title of own data file (replacing 
* Mullin_Snyder_testdata.xlsx) that has three columns: year, price, and 
* quantity. Program assumes the first row of the Excel file has labels "year",
* "p", and "q", which it takes as variable names. If user has different 
* variable names, will have to add commands to rename the variables to year, p,
* and q. If user has a Stata dataset, the command invoking importation of Excel
* will have to be modified. 
*
* Mullin & Snyder December 2019
*
********************************************************************************

*** SET STATA ENVIRONMENTAL PARAMETERS ***

set seed 1776
set more 1
set type double
set maxvar 15000

********************************************************************************

*** SET BOOTSTRAPPING PARAMETERS ***

local bs = "10000" /* Bootstrap repetitions */
local k = "4" /* Number added spacings */
local kplus1 = `k' + 1

********************************************************************************

*** READ IN DATA ***

* As written, program uses cleaned test data from Mullin & Snyder's low-density
* polyethylene market. Researcher should supply own Excel file or other data
* source with the three cleaned variables year, p, and q.

import excel using Mullin_Snyder_testdata.xlsx, firstrow
summ

********************************************************************************

* CHECK FOR AND DROP VIOLATORS OF RECTANGULAR EXPANSION *** 

sort  year
 forvalues i = 1/`=scalar(_N)' {
   gen not_rect`i' = `i' < _n & p[`i'] > p & q[`i'] > q
}
egen not_rect = rowmax(not_rect*)
drop if not_rect[_n] == 1
drop not_*
 
* Create data range
gen t = _n

********************************************************************************

*** METHOD INCORPORATING LOCAL INFORMATION, WITH BOOTSTRAP ***

* Record pairwise comparisons
forvalues i = 1/`=scalar(_N)' {
   quietly gen b_u`i' = .
   quietly replace b_u`i' = abs(q[`i'] - q) / abs(p[`i'] - p) if (`i' < t & q[`i'] < q & p[`i'] > p) | (`i' > t & q[`i'] > q & p[`i'] < p)
   * Work in re_u = 1/e_u space so can take max rather than min
   quietly gen re0_u`i' = 0
   quietly replace re0_u`i' = 1 / (b_u`i' * p / q) if b_u`i' != .
   }

* Generate estimate
egen re_max_u = rowmax(re0_u*) 
gen e_min_u = 1 / re_max_u
drop b_u* re_max_u

* Generate the spacings, denoted d, normalized by "i" multiplication
reshape long re0_u, i(t) j(i)
by t: egen rank = rank(-re0_u), unique
drop i
reshape wide re0_u, i(t) j(rank)

forvalues i = 1/`k' {
   local lead = `i' + 1
   quietly gen d`i' = `i' * (re0_u`i' - re0_u`lead')
}

gen re_new = 0

* Bootstrap ala Zelterman. Note undoes normalization by dividing by "j"
forvalues i = 1/`bs' {
   quietly replace re_new = re0_u`kplus1'
   forvalues j = 1/`k' {
      quietly gen rn = ceil(uniform() * `k')
	  forvalues jj = 1/`k' {
	     quietly replace re_new = re_new + (rn==`jj')*d`jj'/`j'
      }   
      drop rn 
   }
   quietly gen e_bs`i' = 1 / re_new
   }

* Compute various bootstrapped confidence thresholds
egen e_1 = rowpctile(e_bs*), p(1)
egen e_5 = rowpctile(e_bs*), p(5)
egen e_10 = rowpctile(e_bs*), p(10)
egen e_90 = rowpctile(e_bs*), p(90)
egen e_95 = rowpctile(e_bs*), p(95)
egen e_99 = rowpctile(e_bs*), p(99)

* Collect bound estimate and bootstrapped confidence intervalus in file
preserve
rename e_min_u e
gen method = "local"
keep year e e_1 e_5 e_10 e_90 e_95 e_99 method
save boot_results_linear, replace
restore

********************************************************************************

*** METHOD INCORPORATING LIMITING INFORMATION, WITH BOOSTRAP ***

keep year e_min_u re0* p q t
rename e_min_u e0_u
gen b0_u = q * e0_u / p

*** Replace first-stage bound if this second stage results in improvement
forvalues i = 1/`=scalar(_N)' {
   quietly gen b1_u`i' = (q[`i'] - q + b0_u[`i'] * p[`i']) / p
   quietly replace b1_u`i' = q / (p[`i'] - p + q[`i']/b0_u[`i']) if `i' < t & p[`i'] + q[`i'] / b0_u[`i'] > p
   * Work in re_u = 1/e_u space so can take max rather than min
   quietly gen re1_u`i' = 0
   quietly replace re1_u`i' = 1 / (b1_u`i' * p / q) if b1_u`i' != .
}
   
* Generate estimate
egen re_max_u = rowmax(re1_u*) 
gen e_min_u = 1 / re_max_u
drop b* re_max_u

* Generate the spacings, denoted d, normalized by "i" multiplication
reshape long re0_u re1_u, i(t) j(i)
rename re0_u re0
rename re1_u re1
egen group = group(t i)
reshape long re, i(group) j(stage)
bysort t: egen rank = rank(-re), unique
drop i group stage
reshape wide re, i(t) j(rank)

forvalues i = 1/`k' {
   local lead = `i' + 1
   quietly gen d`i' = `i' * (re`i' - re`lead')
}

gen re_new = 0

* Bootstrap ala Zelterman. Note undoes normalization by dividing by "j"
forvalues i = 1/`bs' {
   quietly replace re_new = re`kplus1'
   forvalues j = 1/`k' {
      quietly gen rn = ceil(uniform() * `k')
	  forvalues jj = 1/`k' {
	     quietly replace re_new = re_new + (rn==`jj')*d`jj'/`j'
      }   
      drop rn 
   }
   quietly gen e_bs`i' = 1 / re_new
   }

* Compute various bootstrapped confidence thresholds
egen e_1 = rowpctile(e_bs*), p(1)
egen e_5 = rowpctile(e_bs*), p(5)
egen e_10 = rowpctile(e_bs*), p(10)
egen e_90 = rowpctile(e_bs*), p(90)
egen e_95 = rowpctile(e_bs*), p(95)
egen e_99 = rowpctile(e_bs*), p(99)

* Collect bound estimate and bootstrapped confidence intervalus in file
rename e_min_u e
gen function = "linear"
gen method = "limiting"
keep year e e_1 e_5 e_10 e_90 e_95 e_99 method
append using boot_results_linear
save boot_results_linear, replace

********************************************************************************

*** DISPLAY TABLE OF RESULTS ***

* Can be used to construct table of numerical values and significance levels
* to complement visual results in figure.

gsort method year
list method year e e_90 e_95 e_99 

********************************************************************************

*** GRAPH OF RESULTS FROM METHOD INCORPORATING LOCAL INFORMATION ***

keep if method == "local"

* Compute common vertical axis cutoff (given by worst case across results) so 
* different graphs can be comparable
egen ymax1 = rowmax(e*)
egen ymax2 = max(ymax1)
gen ymax3 = ceil(ymax2)
local ymax = ymax3

egen xmin = min(year)
local xmin = xmin
egen xmax = max(year)
local xmax = xmax

#delimit;

twoway (bar e year, fcolor(gs8) lcolor(black) barwidth(.70)) 
	   (rcap e_5 e_95 year, lcolor(black)),
	   xscale(lwidth(medthin))
	   yscale(range(0 `ymax') lwidth(medthin))
	   xlabel(`xmin'(1)`xmax', noticks angle(90) labsize(small))
	   ylabel(0(1)`ymax', angle(0) labsize(small) glcolor(gs14))
   	   xtitle("Year", size(small))
	   ytitle("Elasticity", size(small))
	   title("Elasticity Bounds Assuming Linear Demand Using Method Incorporating Local Information", size(medsmall))
	   legend(cols(1) order(1 2) lab(1 "Elasticity bounds") lab(2 "Bootstrap c.i.") symxsize(4) size(small) ring(0) bplacement(neast))
	   plotregion(lcolor(black) lwidth(medthin) fcolor(white)) graphregion(color(white)) bgcolor(white);

#delimit cr
graph export "results_linear_local.eps", logo(off) fontface("Calibri") replace

********************************************************************************

*** GRAPH OF RESULTS FROM METHOD INCORPORATING LIMITING INFORMATION ***

clear
use boot_results_linear 
keep if method == "limiting"

#delimit;

twoway (bar e year, fcolor(gs8) lcolor(black) barwidth(.70)) 
	   (rcap e_5 e_95 year, lcolor(black)),
	   xscale(lwidth(medthin))
	   yscale(range(0 `ymax') lwidth(medthin))
	   xlabel(`xmin'(1)`xmax', noticks angle(90) labsize(small))
	   ylabel(0(1)`ymax', angle(0) labsize(small) glcolor(gs14))
   	   xtitle("Year", size(small))
	   ytitle("Elasticity", size(small))
	   title("Elasticity Bounds Assuming Linear Demand Using Method Incorporating Limiting Information", size(medsmall))
	   legend(cols(1) order(1 2) lab(1 "Elasticity bounds") lab(2 "Bootstrap c.i.") symxsize(4) size(small) ring(0) bplacement(neast))
	   plotregion(lcolor(black) lwidth(medthin) fcolor(white)) graphregion(color(white)) bgcolor(white);
    
#delimit cr
graph export "results_linear_limiting.eps", logo(off) fontface("Calibri") replace