*** CLEAR
set more off
capture clear *
capture log close
capture estimates clear
set memory 2000m
set matsize 2500
set linesize 200
//ssc install ivreg2
sysdir set PLUS "D:\Stata13 x64\ado\plus"

di as result "ESTIMATES OF IMPORT ELASTICITIES ARE PRODUCED USING THE METHOD DEVELOPED IN" _newline ///
				"	Soderbery, Anson, �Estimating Import Supply and Demand Elasticities: Analysis and Implications,� 2013." _newline ///
 				"PLEASE CITE ACCORDINGLY"
				
*** USER DEFINED VALUE, QUANTITY, AGGREGATION, COUNTRY, TIME PERIOD, DIRECTORY AND DATA
		global value   = "value_qy_y"
		global quan    = "value_qt_y"
		global good    = "hs_code"
		global country = "country"
		global time    = "year"
		global minyear = 2002
		global DIRECTORY = "E:\Research\GVCs\data\Es_Sigma\sigma"
		global trade_data = "2002-2007_cus_year_handled"

cd "$DIRECTORY" 


*** INITIALIZE NONLINEAR LIML PROGRAM
do "$DIRECTORY/mata_LIMLhybrid".do

*** LOOP OVER HS's *****
local counter = 1
use "$DIRECTORY/$trade_data.dta" if $value > 0 & $quan > 0, clear
levelsof $good, local(hs)
cap foreach good in `hs' {
	noisily di as txt "GOOD = `good'"
use "$DIRECTORY/$trade_data.dta" if $good=="`good'" & $value > 0 & $quan > 0, clear


****** LEVEL OF AGGREGATION, USER SPECIFIED AND STOP IF DATA UNAVAILABLE *******
rename $value  cusval             
rename $quan quantity
tostring $country $good, replace force
collapse (sum) quantity cusval, by($time $good $country) 
egen product = group($country $good)	/** PRODUCT IS YOUR VARIETY. IN F/BW VARIETY IS AN HS10 COUNTRY PAIR **/
bys product: egen periods = nvals($time)
egen maxperiods = max(periods)
egen products   = nvals(product)

if (products < 3 | maxperiods < 2) {
   noisily di as error "BAD APPLE: products < 3 | maxperiods < 2"
   continue
}

*** LOG UNIT VALUES AND SHARES 
drop if $time < $minyear
gen t = $time - $minyear + 1
gen meanprice = cusval/quantity

sort t $good
by t $good: egen totsum = sum(cusval)
gen s  = cusval/totsum
gen ls = ln(s)
gen lp = ln(meanprice)

*** FIRST DIFFERENCING
sort product t
by product: gen ls_dif = ls[_n]-ls[_n-1]
by product: gen lp_dif = lp[_n]-lp[_n-1]

sort t
save "$DIRECTORY/mfg_cusimport_hybrid.dta", replace


*** AUTOMATE CHOOSING REFERENCE VARIETY k (need to download xfill)
di "REFERENCE GOOD"
xtset product
bysort product: gen period = _N
sum period, detail
local maxt = r(max)

sum cusval, detail
local q_max = r(max)
local q_cutoff = r(p90)

gen ref =.
replace ref = product if period>=`maxt' & cusval>=`q_max'
xfill ref, i(product)
xfill ref, i(t)
local ref = ref

replace ref = product if ref == . & period>=`maxt' & cusval>=`q_cutoff'
xfill ref, i(product)
xfill ref, i(t)
local ref = ref

replace ref = product if ref == . & period>=`maxt' 
xfill ref, i(product)
xfill ref, i(t)
local ref = ref

gsort -ref
local ref = ref

gen ref_country = $country if ref == product
gsort -ref_country
local ref_country = ref_country

*** REFERENCE PRODUCT ***
keep if product==`ref'

rename ls_dif h_ls_dif
rename lp_dif h_lp_dif
keep t h_ls_dif h_lp_dif

sort t
save "$DIRECTORY/mfg_cref_hybrid.dta", replace
use "$DIRECTORY/mfg_cusimport_hybrid.dta", clear
merge m:1 t using "$DIRECTORY/mfg_cref_hybrid.dta"
tab _merge
drop _merge

/*********************** IV REGRESSION *************************************/

gen y  = (lp_dif-h_lp_dif)^2
gen x1 = (ls_dif-h_ls_dif)^2
gen x2 = (lp_dif-h_lp_dif)*(ls_dif-h_ls_dif)
drop if y == . | x1 == . | x2 == . 

save "$DIRECTORY/mfg_cregready_hybrid.dta", replace  /* THIS DATASET IS CREATED TO CHECK ANY ERRORS FOR THE CURRENT GOOD */

/****** IV 1 ****/
noisily di "  FIRST STEP"
sort product t
qui tab product if product != `ref', gen(c_I_)  
cap foreach var of varlist x1 x2 {
	regress `var' c_I_*, noc
	predict `var'hat
}

		if (_rc==2001 | _rc == 111) {
  			noisily di as error "ERROR: _rc = " _rc
  			continue
		}

cap ivreg2 y (x1 x2 = c_I_*), liml fuller(1)		/* CONSISTENT, BUT NOT EFFICIENT ESTIMATES */
if (_rc==504 | _rc==481 | _rc==3351) {
  noisily di as error "ERROR: _rc = " _rc
  continue
}


/**** OPTIMAL WEIGHTS ****/
cap predict uhat, resid
		if (_rc==2001 | _rc == 198) {
  			noisily di as error "_rc = " _rc
  			continue
		}
gen uhat2=uhat^2
drop if uhat2==.


qui regress uhat2 c_I_*, noc
predict uhat2hat
gen shat=sqrt(uhat2hat)

/**************************************************************************************
bysort product: gen period = _N							*BW 2006 WGHT OPTIONAL*
sort cname t
by cname: gen bw_wght = (period^1.5) * (1/quantity[_n] + 1/quantity[_n-1])^(-.5) 
****************************************************************************************/

********* WEIGHT DATA AND REESTIMATE ******
noisily di "     WEIGHT and REESTIMATE"
gen ones=1
foreach var of varlist y x1hat x2hat ones {
	gen `var'star = `var'/shat
}

cap ivreg2 y (x1 x2 = c_I_*) [aw = 1/uhat2hat], liml fuller(1)    
if (_rc==504 | _rc==481 | _rc == 3351 | _rc == 2000 | _rc == 198) {
   noisily di as error "_rc = " _rc
   continue
}

gen double theta1 = _b[x1]
gen double theta2 = _b[x2]
local cons = _b[_cons]
gen cons = `cons'
*** Calculate sigma and rho ***

gen double rho1=.5+sqrt(.25-1/(4+theta2^2/theta1)) if theta2>0
replace rho1=.5-sqrt(.25-1/(4+theta2^2/theta1)) if theta2<0
gen double sigma1=1+(2*rho-1)/((1-rho)*theta2)


foreach keep in theta1 theta2 rho1 sigma1 {
local `keep'_hat = `keep'
}
local omega_hat = `rho1_hat'/(`sigma1_hat' - 1 - (`sigma1_hat' * `rho1_hat'))
noisily di as result "  SIGMA = `sigma1_hat'    OMEGA = `omega_hat'   RHO = `rho1_hat'"
local constrained = 0
save "$DIRECTORY/mfg_cregready_hybrid.dta", replace

****** CONSTRAINED SEARCH IF NECESSARY ******
if (`theta1_hat' < 0 | `sigma1_hat' == . ) {
noisily di as txt "     INITIALIZE NONLINEAR LIML SEARCH"
local constrained = 1
local sigma1_hat = .
local rho1_hat = .

drop if c_I_1 == .


mata: Y = st_data(.,"y")
mata: X = st_data(.,("ones", "x1","x2"))
mata: Z = st_data(.,"c_I*")

mata: s = optimize_init()
mata:	optimize_init_evaluator(s,&LIML2())
mata:	optimize_init_evaluatortype(s, "v0")
mata:	d_start = (`cons', 1.5, .4)
mata:	optimize_init_params(s,d_start)
mata:	optimize_init_argument(s, 1, Y)
mata:	optimize_init_argument(s, 2, X)
mata:	optimize_init_argument(s, 3, Z)
mata:	optimize_init_which(s,"min")
mata:	optimize_init_conv_maxiter(s , 500)		
di "Optionally: If the constrained search is particularly time consuming the user may wish to truncate the maximum iterations (code line 234)"
cap mata: d = optimize(s)

if (_rc==430) {
   di as error "_rc = " _rc
   continue
}

mata: Qstar = optimize_result_value(s)
mata: st_matrix("d" , d)

	local cons           = d[1,1]
	local sigma_nl_hat   = d[1,2]
	local rho_nl_hat     = d[1,3]
	foreach var in sigma rho {
		capture replace `var'1 = ``var'_nl_hat' 
 		local `var'1_hat = ``var'_nl_hat' 
	}
	local omega_hat = `rho1_hat'/(`sigma1_hat' - 1 - (`sigma1_hat' * `rho1_hat'))
	if `omega_hat' < 0 {
		local omega_hat = 0
	}
	noisily di as result "  SIGMA = `sigma1_hat'    OMEGA = `omega_hat'   RHO = `rho1_hat'"

} 
	**** END ESTIMATION LOOP

** COLLECT ESTIMATES
drop _all
set obs 1
foreach keep in theta1 theta2 rho1 sigma1 {
capture gen double `keep'=``keep'_hat'
local `keep'_hat=.
}

contract theta* sigma* rho*
capture drop _freq
capture gen omega = rho/(sigma - 1 - (sigma * rho))
replace omega = 0 if omega < 0
capture gen  product = "`good'"
capture gen N = `N'
capture gen cons = `cons'
capture gen ref = "`ref_country'"
capture gen period = `maxt'
capture gen constrained = `constrained'

if (`counter'==1) {
  save "$DIRECTORY/elast_hybrid", replace
  local counter = `counter' + 1
} 
if (`counter'!=1) {
cap append using "$DIRECTORY/elast_hybrid"
duplicates drop *, force
cap save "$DIRECTORY/elast_hybrid", replace
}

}
	**** END GOODS LOOP

di as result "ESTIMATES OF IMPORT ELASTICITIES ARE PRODUCED USING THE METHOD DEVELOPED IN" _newline ///
				"	Soderbery, Anson, �Estimating Import Supply and Demand Elasticities: Analysis and Implications,� 2013." _newline ///
 				"PLEASE CITE ACCORDINGLY"

