program ivexact, eclass
	/* This program performs estimation and nonparametric exact inference
		on the quantiles of potential outcomes for an endogenous binary
		regressor and binary instrument.
		This version uses the ML estimates as a starting place and finds the confidence intervals "from the
		inside out." Hopefully this will lead to much faster estimation. It does, for large sample sizes
		This version also speeds up computation by saving the critical values from each combination of parameters
	*/
	version 12
	syntax varlist(min=3 max=3 numeric) [if] [in]  , [QUANTILES(numlist)] ///
	[LEVEL(real .9)] [CI] [NOAt] [NONt] [NOIsily]

	/*
	 syntax: ivexact depvar assignvar treatvar, [QUANTILES(list of quantile indices)] [CI] [level(confidence level for CI estimation)] [NOAt] [NONt]
	example: ivexact y z d, ci quantiles(.05(.05).95) level(.95) noats
	Comments:
	the treatment and assignment variable must be binary
	if quantile option is omitted, default is inference on the median
	The option ci should be specified if a ocnfidence interval is desired. Default is only to compute ml estimates
	The default level for the CI is .9
	the options noat and nont mean zero always-takers or zero never-takers will be imposed (unless that data contradict that)

	If a CI is computed, this is the interpretation: the true quantile is included in [a,b) with a probability of at least
	the coverage level. Note the asymmetry in the closure of the interval.

	*/
	marksample touse
	preserve
	qui keep if `touse'
	tempname n P Q qlist xs xsci I J crits results0 results1 alpha x fmest pr pr prp T nCminmax nATminmax kATminmax twcrits
	tempvar Y Z D YFM DFM

	tokenize `varlist'

	gen `Y'=`1'
	gen `Z'=`2'
	gen `D'=`3'
	* The following are for estimating conventional (Frolich and Melly) quantiles of compliers' potential outcomes
	qui gen `YFM'=.
	qui gen `DFM'=.


	/*
	The idea behind estimation is maximum likelihood. Let's say I want an
	estimate for Q(tau). I find the combination
	of parameters (nC,kC=ceil(nC*tau),nAT,kAT(x),x) to maximize the likelihood
	of (P,Q,I(x),J(x)). I could do this by starting at the conventional estimates,
	and working out in x until I find a local maximum. And for each x I naturally find
	the (nC,nAT,kAT(x)) that maximizes the likelihood for that x.

	Here is the idea behind inference. The confidence interval for Q(tau) will consist
	of all those values x such that the null hypothesis ceil(nC*tau)=kC(x) is
	not rejected, where kC(x) is the number of
	compliers with potential outcomes less than or equal to x. The algorithm will
	start by looping through all possible values of nC that are consistent with the data.
	For each value of nC, it will set kC(x)=ceil(nC*tau), and will loop through all possible
	values of nAT that are consistent with nC and kC(x) and the data. Within that it will
	loop through all possible values of kAT(x) that are consistent with nC,kC(x), nAT, and the data.
	For each value of nC,kC(x),nAT,kAT(x) it will compute critical values for (P,Q,I(x),J(x)) under
	the null of those parameter values, and compare to the observed values to accept or reject.
	If any parameter sets are not rejected, then x is included in the confidence interval. In practice,
	since the extrema of the confidence interval are really of interest, the algorithm will start with extreme
	values of x and work in and stop as soon as an x is included in the CI.

	Now how to compute critical values for P,Q,I(x),J(x). I think I will do my "around the box" algorithm
	where I start at the extremes, and ratchet each bound in one at a time and compute the new coverage (by
	substracting out the probability of the last step) and stop ratcheting when the coverage drops below alpha.

	*/

	* initialize quantile vector
	if "`quantiles'"=="" {
		if "`noisily'"!="" display "quantiles not specified; using default q=.5"
		matrix `qlist' =.5
		local qlength=1
	}
	else {
		local qlength=0
		foreach quant of numlist `quantiles' {
			local ++qlength
		}
		* Initialize quantile vector
		matrix `qlist'=J(`qlength',1,0)
		* Populate quantile vector
		local q=0
		foreach quant of numlist `quantiles' {
			local ++q
			matrix `qlist'[`q',1]=`quant'
		}
	}
	mat colnames `qlist'=quantile
	if "`ci'"!="" {
		mat `results0'=J(`qlength',3,0)
		mat colnames `results0'=mle lower upper
	}
	else {
		mat `results0'=J(`qlength',1,0)
		mat colnames `results0'=mle
	}
	mat `results0'=`qlist',`results0'
	mat `results1'=`results0'

	if `level'>=1 | `level'<=0 {
		local level=.9
		display "Confidence level out of bounds; default will be used"
	}

	scalar `alpha'=1-`level'

	scalar `n'=_N

	* Check whether imposing noat and nont contradicts the data:
	if "`nont'"!="" {
		qui count if `Z'==1 & `D'==0
		if r(N)!= 0 display "No never-takers imposed, but contradicted by data. Ignoring nont option."
		else local nNTimp1 `nont'
	}
	if "`noat'"!="" {
		qui count if `Z'==0 & `D'==1
		if r(N)!= 0 display "No always-takers imposed, but contradicted by data. Ignoring noat option."
		else local nATimp1 `noat'
	}
	local nATimp0 `nNTimp1'
	local nNTimp0 `nATimp1'
	local bound0 lower
	local bound1 upper


	* loops through treated and untreated
	forvalues d=0/1 {
		if "`noisily'"!="" display "Estimating quantiles of y(`d')"
		* generate Frolich and Melly outcome and treatment
		qui replace `YFM'=`Y'*(`D'==`d')
		qui replace `DFM'=(`D'==`d')
		* compute P and Q
		qui count if `Z'==`d' & `D'==(1-`d')
		scalar `P'=r(N)

		qui count if `Z'==(1-`d') & `D'==`d'
		scalar `Q'=r(N)

		qui count if `Z'==`d'
		local m=r(N)
		* define starting evaluation vector (all observed values of Y(d))
		qui tab `Y' if `D'==`d',matrow(`xs')
		local numxs=rowsof(`xs')
		display "values to loop over: "
		mat li `xs'
		forvalues qind=1/`qlength' {
			local tau=`qlist'[`qind',1]
			if "`noisily'"!="" display "On quantile index = " `tau'
			if "`noisily'"!="" display "Computing maximum likelihood estimate"
			* First maximum likelihood estimation
			* Find the conventional estimate using Frolich and Melly's endogenous QTE estimator
			qui ivqte `YFM' (`DFM'=`Z'),quantiles(`tau')
			mat `fmest'=e(b)

			local ltx=0
			local xind=0
			while `ltx'==0 {
				local ++xind
				if `xind'<`numxs' local ltx=`fmest'[1,1]<(`xs'[`xind',1]-.0000001)
				else local ltx=1
			}

			* find maximum likelihood at x0=xs[xind0,1]
			* compute I(x0) and J(x0)
			qui count if `Z'==(1-`d') & `D'==`d' & `Y'<=(`xs'[`xind',1]+.00001)
			scalar `I'=r(N)

			qui count if `Z'==`d' & `D'==`d' & `Y'<=(`xs'[`xind',1]+.00001)
			scalar `J'=r(N)

			mat `T'=(`P',`Q',`I',`J')
			mata st_numscalar("`pr'" ///
			, maxjointpr(st_matrix("`T'") ///
			,strtoreal("`tau'") ///
			,st_numscalar("`n'") ///
			,`m' ///
			,"`nNTimp`d''", "`nATimp`d''"))

			* find maximum likelihood at xp=xs[xind0-1,1]
			* compute I(xp) and J(xp)
			qui count if `Z'==(1-`d') & `D'==`d' & `Y'<=(`xs'[`xind'-1,1]+.00001)
			scalar `I'=r(N)
			qui count if `Z'==`d' & `D'==`d' & `Y'<=(`xs'[`xind'-1,1]+.00001)
			scalar `J'=r(N)
			mat `T'=(`P',`Q',`I',`J')
			mata st_numscalar("`prp'" ///
			, maxjointpr(st_matrix("`T'") ///
			,strtoreal("`tau'") ///
			,st_numscalar("`n'") ///
			,`m' ///
			,"`nNTimp`d''", "`nATimp`d''"))

			* if prp is greater than pr0, keep looking in the downward direction until pr decreases. If prp is less than pr0, look
			* upwards until pr decreases
			if `prp'>=`pr' {
				local step=-1
				scalar `pr'=`prp'
				local xind=`xind'-1
			}
			else {
				local step=1
				scalar `prp'=`pr'
			}

			while `prp'>=`pr' & `xind'>0 & `xind'<=`numxs' {
				scalar `pr'=`prp'
				local xindp=`xind'+`step'
				* compute I(xp) and J(xp)
				qui count if `Z'==(1-`d') & `D'==`d' & `Y'<=(`xs'[`xindp',1]+.00001)
				scalar `I'=r(N)
				qui count if `Z'==`d' & `D'==`d' & `Y'<=(`xs'[`xindp',1]+.00001)
				scalar `J'=r(N)
				mat `T'=(`P',`Q',`I',`J')
				mata st_numscalar("`prp'" ///
				, maxjointpr(st_matrix("`T'") ///
				,strtoreal("`tau'") ///
				,st_numscalar("`n'") ///
				,`m' ///
				,"`nNTimp`d''", "`nATimp`d''"))
				local xind=`xindp'
			}
			* backtrack one to find the optimal x:
			local xindmax=`xind'-`step'
			matrix `results`d''[`qind',2]=`xs'[`xindmax',1]

			if "`ci'"!="" {
				if "`noisily'"!="" display "Computing exact confidence intervals"
				* concatenate a low value to the beginning of my xs matrix, because it it possible to test a value
				* strictly less than all observed outcome values (I=0 and J=0), but of course it is redundant
				* to test a value greater than all observed outcome values (because I=Q and J=m-P-Q
				* no matter how high x is
				matrix `xsci'=(`xs'[1,1]-1 \ `xs')
				local numxsci=rowsof(`xsci')

				* find lower bound, then upper bound, starting at the ML estimate
				forvalues upper=0/1 {
					* initialize xind so that the first x tried will be the MLE (to account for concatenation of xs)
					local xind=`xindmax'+1+(-1)^`upper'
					if "`noisily'"!="" display "Finding `bound`upper'' bound"
					local reject=0

					while `reject'==0 & `xind'>1 & `xind'<`numxsci' {
						* I initialize reject = 1 because below I have statements that kick me out of the loop as
						* soon as reject = 0 (because once I have failed to reject some combination of parameters,
						* I fail to reject the x value (this is what's behind the incidental parameters problem)
						local reject = 1

						* increment/decrement xind depending on if I'm working on upper or lower bound:
						local xind=`xind'-(-1)^`upper'
						scalar `x'=`xsci'[`xind',1]
						display "xind = `xind', x = " `x'

						* compute I(x) and J(x)
						qui count if `Z'==(1-`d') & `D'==`d' & `Y'<=(`x'+.00001)
						scalar `I'=r(N)
						qui count if `Z'==`d' & `D'==`d' & `Y'<=(`x'+.00001)
						scalar `J'=r(N)

						* calculate bounds using mata subroutines:
						mata st_matrix("`nCminmax'" ///
						,nCbounds(st_numscalar("`P'") ///
						,st_numscalar("`Q'") ///
						,st_numscalar("`n'") ///
						,`m' ///
						,"`nNTimp`d''", "`nATimp`d''"))

						local nCmin=`nCminmax'[1,1]
						local nCmax=`nCminmax'[1,2]

						forvalues nC=`nCmin'/`nCmax' {
							local kC=ceil(`nC'*`tau')

							* calculate bounds using mata subroutines:
							mata st_matrix("`nATminmax'" ///
							,nATbounds(st_numscalar("`P'") ///
							,st_numscalar("`Q'") ///
							,st_numscalar("`n'") ///
							,`m' ///
							,strtoreal("`nC'") ///
							,"`nNTimp`d''", "`nATimp`d''"))
							local nATmin=`nATminmax'[1,1]
							local nATmax=`nATminmax'[1,2]

							forvalues nAT=`nATmin'/`nATmax' {
								* calculate bounds using mata subroutines:
								mata st_matrix("`kATminmax'" ///
								,kATbounds(st_numscalar("`P'") ///
								,st_numscalar("`Q'") ///
								,st_numscalar("`I'") ///
								,st_numscalar("`J'") ///
								,`m' ///
								,strtoreal("`nC'") ///
								,strtoreal("`kC'") ///
								,strtoreal("`nAT'")))
								local kATmin=`kATminmax'[1,1]
								local kATmax=`kATminmax'[1,2]

								forvalues kAT=`kATmin'/`kATmax' {
									capture mat `crits'=`crits`nC'`kC'`nAT'`kAT'`m''
									if _rc!=0 {
										mata crits=critvals(st_numscalar("`alpha'"),`nC',`nAT',`kC',`kAT',st_numscalar("`n'"),`m')
										mata st_matrix("`crits'",crits)
										tempname crits`nC'`kC'`nAT'`kAT'`m'
										mat `crits`nC'`kC'`nAT'`kAT'`m''=`crits'
									}
									display "nC = `nC',kC = `kC', nAT = `nAT', kAT = `kAT'"
									display "(T,crits):"
									mat `twcrits'=((`P' \ . \ `Q' \ . \ `I' \ . \ `J' \ .),`crits'')
									mat li `twcrits'
									local reject=	  `P'<`crits'[1,1] | `P'>`crits'[1,2] ///
												| `Q'<`crits'[1,3] | `Q'>`crits'[1,4] ///
												| `I'<`crits'[1,5] | `I'>`crits'[1,6] ///
												| `J'<`crits'[1,7] | `J'>`crits'[1,8]
									if `reject'==0 continue , break
								}
								if `reject'==0 continue , break
							}
							if `reject'==0 continue , break
						}
					}
					* If I couldn't reject the extremes of my ci evaluation vector x, then bounds are infinite
					if `reject'==0 & (`xind'==1 | `xind'==`numxsci') matrix `results`d''[`qind',`upper'+3]=.
					else {
						* I get to backtrack if lower to find the lowest not rejected:
						scalar `x'=`xsci'[`xind'+1-`upper',1]
						matrix `results`d''[`qind',`upper'+3]=`x'
					}
				}

			}
		}
	}
	ereturn post,esample(`touse')
	ereturn matrix results1=`results1'
	ereturn matrix results0=`results0'
	ereturn scalar level=`level'
	ereturn local cmd "ivexact"
	display in gr _n _n "Exact Instrumental Variables Results"
	if "`ci'"!="" display in gr _n _n "Maximum likelihood estimates and `level'-level confidence interval for quantiles of Y(0)"
	else display in gr _n _n "Maximum likelihood estimates for quantiles of Y(0)"
	matlist e(results0),names(column)
	if "`ci'"!="" display in gr _n _n "Maximum likelihood estimates and `level'-level confidence interval for quantiles of Y(1)"
	else display in gr _n _n "Maximum likelihood estimates for quantiles of Y(1)"
	matlist e(results1),names(column)

	restore
end