capt program drop cmi_interval
program define cmi_interval, rclass

version 11.1

// Programmed by Wooyoung Kim (University of Wisconsin-Madison, wkim68@wisc.edu)
// This is the command implementing "Inference Based on Conditional Moment Inequalities" (Donald W.K. Andrews and Xiaoxia Shi, 2013).

/*
<OUTPUT LIST>




<INPUT LIST>

*/



local lbvar ` '
local ubvar ` '

gettoken lbvar 0: 0, match(leftover)
gettoken ubvar 0: 0, match(leftover)


syntax anything [if] [in] [, LEVel(real 0.95) DECI(integer 3) *]


local max_deci `deci'
local reject = 1 - `level'
local bon_reject = 1 - (1 + `level')/2

local num_lb = wordcount("`lbvar'")
local num_ub = wordcount("`ubvar'")

if `num_lb' != 0 {
	foreach a of local lbvar {
		quietly sum `a'
		local lowmin = `r(min)'
		local lowmax = `r(max)'
	}

	foreach a of local lbvar {
		quietly sum `a'
		local lowmin = min(`lowmin',`r(min)')
		local lowmax = max(`lowmax',`r(max)')
	}

	local ldist = `lowmax' - `lowmin'
	local lowmin = floor(`lowmin')
	local lowmax = ceil(`lowmax')

}

if `num_ub' != 0 {
	foreach a of local ubvar {
		quietly sum `a'
		local upmin = `r(min)'
		local upmax = `r(max)'
	}
	foreach a of local ubvar {
		quietly sum `a'
		local upmin = min(`upmin',`r(min)')
		local upmax = max(`upmax',`r(max)')
	}

	local udist = `upmax' - `upmin'
	local upmin = floor(`upmin')
	local upmax = ceil(`upmax')


}




if `num_lb' == 0 {


	/********************************************************************/
	/******************* Upper Bound Only Case **************************/
	/********************************************************************/

	local dir = -1
	local ujump = 1
	local ubound `upmax'
	local inbound = 0
	local pinbound = 0
	if `udist' > 1 {
		while (`inbound' != 0 | `pinbound' != 1 | `ujump' != 1){


			local pinbound = `inbound'
			local ub_temp ` '

			foreach a of local ubvar {
				tempvar ubound_theta`i'
				gen `ubound_theta`i'' = `a' - `ubound'
				local dummy "`ubound_theta`i''"
				local ub_temp `ub_temp' `dummy'
				local i = `i' + 1
			}

			quietly cmi_test (`ub_temp') () `anything' `if' `in' , `options'

			local pval = `r(pval)'

			if `dir' == -1{
				if `pval' < `reject'{
					local ujump = `ujump' * 2
					local inbound = 0
				}
				else{
					local ujump = max(`ujump' /2, 1)
					local dir = -1 * `dir'
					local inbound = 1
				}
			}
			else{
				if `pval' >= `reject'{
					local ujump = `ujump' * 2
					local inbound = 1
				}
				else{
					local ujump = max(`ujump' /2, 1)
					local dir = -1 * `dir'
					local inbound = 0
				}
			}

			local ubound = `ubound' + `dir' * `ujump'

		}
	}

	forval deci = 1(1)`max_deci'{

		local ujump = 10^(-`deci')
		local inbound = 0

		while (`inbound' != 1){
			local ub_temp ` '

			foreach a of local ubvar {
				tempvar ubound_theta`i'
				gen `ubound_theta`i'' = `a' - `ubound'
				local dummy "`ubound_theta`i''"
				local ub_temp `ub_temp' `dummy'
				local i = `i' + 1
			}

			quietly cmi_test (`ub_temp') () `anything' `if' `in' , `options'
			local pval = `r(pval)'

			if `pval' > `reject' {
				local inbound = 1
			}
			else{
				local ubound = `ubound' - `ujump'
			}

		}

		local ubound = `ubound' + `ujump'
	}


}
else if `num_ub' == 0 {

	/********************************************************************/
	/******************* Lower Bound Only Case **************************/
	/********************************************************************/


	local dir = 1
	local ljump = 1
	local lbound `lowmin'
	local inbound = 0
	local pinbound = 0
	if `ldist' > 1 {
		while (`inbound' != 0 | `pinbound' != 1 | `ljump' != 1){

			local pinbound = `inbound'
			local lb_temp ` '

			foreach a of local lbvar {
				tempvar lbound_theta`i'
				gen `lbound_theta`i'' = `lbound' - `a'
				local dummy "`lbound_theta`i''"
				local lb_temp `lb_temp' `dummy'
				local i = `i' + 1
			}

			quietly cmi_test (`lb_temp') () `anything' `if' `in' , `options'

			local pval = `r(pval)'

			if `dir' == 1{
				if `pval' < `reject'{
					local ljump = `ljump' * 2
					local inbound = 0
				}
				else{
					local ljump = max(`ljump' /2, 1)
					local dir = -1 * `dir'
					local inbound = 1
				}
			}
			else{
				if `pval' >= `reject'{
					local ljump = `ljump' * 2
					local inbound = 1
				}
				else{
					local ljump = max(`ljump' /2, 1)
					local dir = -1 * `dir'
					local inbound = 0
				}
			}

			local lbound = `lbound' + `dir' * `ljump'

		}
	}

	forval deci = 1(1)`max_deci'{
		local ljump = 10^(-`deci')
		local inbound = 0

		while (`inbound' != 1){
			local lb_temp ` '

			foreach a of local lbvar {
				tempvar lbound_theta`i'
				gen `lbound_theta`i'' = `lbound' - `a'
				local dummy "`lbound_theta`i''"
				local lb_temp `lb_temp' `dummy'
				local i = `i' + 1
			}

			quietly cmi_test (`lb_temp') () `anything' `if' `in' , `options'
			local pval = `r(pval)'

			if `pval' > `reject' {
				local inbound = 1
			}
			else{
				local lbound = `lbound' + `ljump'
			}

		}

		local lbound = `lbound' - `ljump'

	}

}
else{
    /********************************************************************/
	/************** Both Upper and Lower Bound Case *********************/
	/********************************************************************/

	/* STEP1: Bonferroni Upper Bound */

	local dir = -1
	local ujump = 1
	local ubound `upmax'
	local inbound = 0
	local pinbound = 0

	if `udist' > 1{
		while (`inbound' != 0 | `pinbound' != 1 | `ujump' != 1){

			local pinbound = `inbound'
			local ub_temp ` '

			foreach a of local ubvar {
				tempvar ubound_theta`i'
				gen `ubound_theta`i'' = `a' - `ubound'
				local dummy "`ubound_theta`i''"
				local ub_temp `ub_temp' `dummy'
				local i = `i' + 1
			}

			quietly cmi_test (`ub_temp') () `anything' `if' `in' , `options'

			local pval = `r(pval)'

			if `dir' == -1{
				if (`pval' < `bon_reject'){
					local ujump = `ujump' * 2
					local inbound = 0
				}
				else{
					local ujump = max(`ujump' /2, 1)
					local dir = -1 * `dir'
					local inbound = 1
				}
			}
			else{
				if (`pval' >= `bon_reject'){
					local ujump = `ujump' * 2
					local inbound = 1
				}
				else{
					local ujump = max(`ujump' /2, 1)
					local dir = -1 * `dir'
					local inbound = 0
				}
			}

			local ubound = `ubound' + `dir' * `ujump'

		}
	}

	forval deci = 1(1)`max_deci'{
		local ujump = 10^(-`deci')
		local inbound = 0

		while (`inbound' != 1){
			local ub_temp ` '

			foreach a of local ubvar {
				tempvar ubound_theta`i'
				gen `ubound_theta`i'' = `a' - `ubound'
				local dummy "`ubound_theta`i''"
				local ub_temp `ub_temp' `dummy'
				local i = `i' + 1
			}

			quietly cmi_test (`ub_temp') () `anything' `if' `in' , `options'
			local pval = `r(pval)'

			if (`pval' > `bon_reject') {
				local inbound = 1
			}
			else{
				local ubound = `ubound' - `ujump'
			}

		}

		local ubound = `ubound' + `ujump'
	}


	/* Bonferroni Lower Bound */

	local dir = 1
	local ljump = 1
	local lbound `lowmin'
	local inbound = 0
	local pinbound = 0

	if `ldist' > 1 {
		while (`inbound' != 0 | `pinbound' != 1 | `ljump' != 1){

			local pinbound = `inbound'
			local lb_temp ` '

			foreach a of local lbvar {
				tempvar lbound_theta`i'
				gen `lbound_theta`i'' = `lbound' - `a'
				local dummy "`lbound_theta`i''"
				local lb_temp `lb_temp' `dummy'
				local i = `i' + 1
			}

			quietly cmi_test (`lb_temp') () `anything' `if' `in' , `options'

			local pval = `r(pval)'

			if `dir' == 1{
				if `pval' < `bon_reject'{
					local ljump = `ljump' * 2
					local inbound = 0
				}
				else{
					local ljump = max(`ljump' /2, 1)
					local dir = -1 * `dir'
					local inbound = 1
				}
			}
			else{
				if (`pval' >= `bon_reject'){
					local ljump = `ljump' * 2
					local inbound = 1
				}
				else{
					local ljump = max(`ljump' /2, 1)
					local dir = -1 * `dir'
					local inbound = 0
				}
			}

			local lbound = `lbound' + `dir' * `ljump'

		}
	}

	forval deci = 1(1)`max_deci'{
		local ljump = 10^(-`deci')
		local inbound = 0

		while (`inbound' != 1){
			local lb_temp ` '

			foreach a of local lbvar {
				tempvar lbound_theta`i'
				gen `lbound_theta`i'' = `lbound' - `a'
				local dummy "`lbound_theta`i''"
				local lb_temp `lb_temp' `dummy'
				local i = `i' + 1
			}

			quietly cmi_test (`lb_temp') () `anything' `if' `in' , `options'
			local pval = `r(pval)'

			if (`pval' > `bon_reject') {
				local inbound = 1
			}
			else{
				local lbound = `lbound' + `ljump'
			}

		}

		local lbound = `lbound' - `ljump'

	}


	/* STEP2: Calculate the strict bound */

	if (`lbound' > `ubound') {
		local empty = 1
	}
	else{
		local empty = 0
		local dist = `ubound' - `lbound'
		local jump = max(floor(`dist'/20),1)

		if (`dist' > 2){
			local inbound = 0
			while (`inbound' == 0 | `jump' > 1 ){
				local i = 1
				local lb_temp = ` '

				foreach a of local lbvar {
					tempvar lbound_theta`i'
					gen `lbound_theta`i'' = `lbound' - `a'
					local dummy "`lbound_theta`i''"
					local lb_temp `lb_temp' `dummy'
					local i = `i' + 1
				}

				local i = 1

				foreach a of local ubvar {
					tempvar ubound_theta`i'
					gen `ubound_theta`i'' = `a' - `lbound'
					local dummy "`ubound_theta`i''"
					local lb_temp `lb_temp' `dummy'
					local i = `i' + 1

				}

				quietly cmi_test (`lb_temp') () `anything' `if' `in' , `options'
				local pval = `r(pval)'
				if (`pval' > `reject') {
					local inbound = 1
					local lbound = `lbound' - `jump'
					local jump = max(floor(`jump'/2),1)

				}
				else if (`prepval' > `pval'){
					local inbound = 0.5
					local lbound = `lbound' - 2 * `jump'
					local jump = max(floor(`jump'/2),1)

				}
				else{
					local lbound = `lbound' + `jump'
					local inbound = 0
				}

			}

			if (`inbound' == 0.5) {
				local ubound = `lbound' + 2
			}
			else{
				local inbound = 0
				local dist = `ubound' - `lbound'
				local jump = max(floor(`dist'/20),1)
				while (`inbound' == 0 | `jump' > 1 ){
					local i = 1
					local ub_temp = ` '

					foreach a of local lbvar {
						tempvar lbound_theta`i'
						gen `lbound_theta`i'' = `ubound' - `a'
						local dummy "`lbound_theta`i''"
						local ub_temp `ub_temp' `dummy'
						local i = `i' + 1
					}

					local i = 1

					foreach a of local ubvar {
						tempvar ubound_theta`i'
						gen `ubound_theta`i'' = `a' - `ubound'
						local dummy "`ubound_theta`i''"
						local lb_temp `lb_temp' `dummy'
						local i = `i' + 1

					}

					quietly cmi_test (`ub_temp') () `anything' `if' `in' , `options'
					local pval = `r(pval)'
					if (`pval' > `reject') {
						local inbound = 1
						local ubound = `ubound' + `jump'
						local jump = max(floor(`jump'/2),1)

					}
					else{
						local ubound = `ubound' - `jump'
						local inbound = 0
					}

				}
			}


		}

		forval deci = 1(1)`max_deci'{

			local dist = `ubound' - `lbound'
			local jump = 10^(-`deci')

			if (`dist' > 20*`jump'){
				local inbound = 0
				while (`inbound' == 0){
					local i = 1
					local lb_temp ` '

					foreach a of local lbvar {
						tempvar lbound_theta`i'
						gen `lbound_theta`i'' = `lbound' - `a'
						local dummy "`lbound_theta`i''"
						local lb_temp `lb_temp' `dummy'
						local i = `i' + 1
					}

					local i = 1

					foreach a of local ubvar {
						tempvar ubound_theta`i'
						gen `ubound_theta`i'' = `a' - `lbound'
						local dummy "`ubound_theta`i''"
						local lb_temp `lb_temp' `dummy'
						local i = `i' + 1

					}

					quietly cmi_test (`lb_temp') () `anything' `if' `in' , `options'
					local pval = `r(pval)'

					if (`pval' > `reject') {
						local inbound = 1
					}
					else{
						local lbound = `lbound' + `jump'
					}

				}
				local lbound = `lbound' - `jump'

				local inbound = 0
				while (`inbound' == 0){
					local i = 1
					local ub_temp ` '

					foreach a of local lbvar {
						tempvar lbound_theta`i'
						gen `lbound_theta`i'' = `ubound' - `a'
						local dummy "`lbound_theta`i''"
						local ub_temp `ub_temp' `dummy'
						local i = `i' + 1
					}

					local i = 1

					foreach a of local ubvar {
						tempvar ubound_theta`i'
						gen `ubound_theta`i'' = `a' - `ubound'
						local dummy "`ubound_theta`i''"
						local ub_temp `ub_temp' `dummy'
						local i = `i' + 1

					}

					quietly cmi_test (`ub_temp') () `anything' `if' `in' , `options'

					local pval = `r(pval)'

					if (`pval' > `reject') {
						local inbound = 1
					}
					else{
						local ubound = `ubound' - `jump'
					}

				}
				local ubound = `ubound' + `jump'


			}
			else if `dist' > `jump'{

				local inbound = 0
				local pval = 0
				local prepval = 0

				while (`inbound' == 0 & `pval' >= `prepval'){

					local prepval = `pval'
					local i = 1
					local lb_temp ` '

					foreach a of local lbvar {
						tempvar lbound_theta`i'
						gen `lbound_theta`i'' = `lbound' - `a'
						local dummy "`lbound_theta`i''"
						local lb_temp `lb_temp' `dummy'
						local i = `i' + 1
					}

					local i = 1

					foreach a of local ubvar {
						tempvar ubound_theta`i'
						gen `ubound_theta`i'' = `a' - `lbound'
						local dummy "`ubound_theta`i''"
						local lb_temp `lb_temp' `dummy'
						local i = `i' + 1

					}

					quietly cmi_test (`lb_temp') () `anything' `if' `in' , `options'

					local pval = `r(pval)'

					/* A case in which we can find a point in the bound */
					if (`pval' > `reject') {
						local inbound = 1
					}
					else if (`pval' >= `prepval'){
						local lbound = `lbound' + `jump'
					}


				}
				local lbound = `lbound' - `jump'

				if (`pval' < `prepval'){
					/* we suspect a small interval here */
					local lbound = `lbound' - `jump'
					local ubound = `lbound' + `jump'
				}
				else{
					local inbound = 0

					while (`inbound' == 0){
						local i = 1
						local ub_temp ` '

						foreach a of local lbvar {
							tempvar lbound_theta`i'
							gen `lbound_theta`i'' = `ubound' - `a'
							local dummy "`lbound_theta`i''"
							local ub_temp `ub_temp' `dummy'
							local i = `i' + 1
						}

						local i = 1

						foreach a of local ubvar {
							tempvar ubound_theta`i'
							gen `ubound_theta`i'' = `a' - `ubound'
							local dummy "`ubound_theta`i''"
							local ub_temp `ub_temp' `dummy'
							local i = `i' + 1

						}

						quietly cmi_test (`ub_temp') () `anything' `if' `in' , `options'
						local pval = `r(pval)'


						if (`pval' > `reject') {
							local inbound = 1
							local ubound = `ubound' + `jump'
						}
						else{
							local inbound = 0
							local ubound = `ubound' - `jump'
						}

					}


				}


			}
		}

	}

}


local N = `r(N)'
local kappa = `r(kappa)'
local ncube = `r(ncube)'
local B = `r(B)'
local epsilon = `r(epsilon)'
local rep_cv = `r(rep_cv)'
local a_obs = `r(a_obs)'
local r_n = `r(r_n)'
local method = "`r(method)'"
local method_CV = "`r(method_CV)'"
local method_FUN = "`r(method_FUN)'"


return clear
return local method = "`method'"
return local method_CV = "`method_CV'"
return local method_FUN = "`method_FUN'"

return local x = "`anything'"
if `num_ub' != 0 {
	return local ubvar = "`ubvar'"
}
if `num_lb' != 0 {
	return local lbvar = "`lbvar'"
}
return local title = "Conditional Moment Inequalities Interval"
return local cmd = "cmi_interval"

return scalar N = `N'
if `num_lb' != 0 {
	return scalar lbound = `lbound'
}
if `num_ub' != 0 {
	return scalar ubound = `ubound'
}
return scalar level = `level'
return scalar kappa = `kappa'
return scalar ncube = `ncube'
return scalar B = `B'
return scalar epsilon = `epsilon'
return scalar rep_cv = `rep_cv'
return scalar a_obs = `a_obs'
return scalar r_n = `r_n'
return scalar deci = `max_deci'


display as text _newline "Conditional Moment Inequalities Interval" _col(59) "Number of obs : " as result r(N)
display as text "{hline 80}"
display as text "<Variables>"
if "`lbvar'" != "" {
	display as text "Variables for the Lower Bound : " as result "`lbvar'"
}
else{
	display as text "A Lower Bound is not computed"
}
if "`ubvar'" != "" {
	display as text "Variables for the Upper Bound : " as result "`ubvar'"
}
else{
	display as text "An Upper Bound is not computed"
}
display as text "Instruments : " as result "`anything'"
display as text "{hline 80}"

display as text "<Methods>"

display as text "`method'"
display as text "`method_CV'"
display as text "`method_FUN'"


display as text "{hline 80}"
display as text "<Results>"
display as result 100*r(level) as text "% confidence interval is:"
if `num_lb' == 0 {
	display as result as text "( " as result "-inf" as text " , " %5.`max_deci'f as result  `ubound' as text " )"
}
else if `num_ub' == 0 {
	display as result as text "( " %5.`max_deci'f as result `lbound' as text " , " as result "inf" as text " )"
}
else {
	if `empty' == 0{
		display as result as text "( " %5.`max_deci'f as result  `lbound' as text " , " as result %5.`max_deci'f `ubound' as text " )"
	}
	else{
		display as text "The confidence interval is empty under " `max_deci' " digits decimal"
	}
}



end