capt program drop clr_p
program define clr_p, eclass

version 11.2

// extract first inequality from the syntax 

local x ` '
local v ` '

local numeq = 1
gettoken eqn`numeq' 0: 0, match(leftover)

tokenize `eqn`numeq'' 	
local y `1'
macro shift

local nindep`numeq' = wordcount("`*'")/2

// Note that the number of independent variable and corresponding range variable should be the same
// If the number of input is not even, the procedure returns error 

if int(`nindep`numeq'') != `nindep`numeq'' {
	display as error "syntax error" 
	error 198
}

// get independent variables x and corresponding range variables v 		

forvalues j = 1/`nindep`numeq'' {
		local x `x' `1'
		macro shift
}
	
forvalues j = 1/`nindep`numeq'' {
		local v `v' `1'
		macro shift
}

// "indep_vector" records the number of indep. variables for each inequality.  

mat indep_vector = `nindep`numeq''

while "`leftover'" == "(" {
		
		local check word("`0'",1)
		
		if `check' ==  "," | `check' == "if" |`check' == "in"|`check' == "" {
			continue, break
		}
		
		// extract `numeq'-th ineqaulity from the syntax 
		
		local numeq = `numeq' + 1
		gettoken eqn`numeq' 0: 0 , match(leftover)
		
		tokenize `eqn`numeq''
		
		local y `y' `1'
		macro shift
		
		local nindep`numeq' = wordcount("`*'")/2
		
		if int(`nindep`numeq'') != `nindep`numeq'' {
			display as error "syntax error" 
			error 198
		}
		
		// get independent variables x and corresponding range variables v 		
		
		forvalues j = 1/`nindep`numeq'' {
			local x `x' `1'
			macro shift
		}
	
		forvalues j = 1/`nindep`numeq'' {
			local v `v' `1'
			macro shift
		}
		
		// "indep_vector" records the number of indep. variables for each inequality.  
		
		mat indep_vector = indep_vector \ `nindep`numeq''
		
}

// adjusting options

syntax [if] [in] [,LOWer LEVel(numlist >0 <1 sort) noAIS RND(integer 10000) *]
marksample touse, nov

// set confidence level vector

if "`level'" == "" {
	local level 0.5 0.9 0.95 0.99
}

local nlevel = wordcount("`level'")

mat level_vector = J(`nlevel',1,0)

forval i = 1/`nlevel' {
	mat level_vector[`i',1] = real(word("`level'",`i'))
}


// AIS option

if "`ais'" != "noais" {
	local vsemtd 1
}
else { 
	local vsemtd 0
}

// give random seed 

// determine lower or upper bound 

if "`lower'" != "lower" {
	local upp 1
}
else {
	local upp 0
}

// obtain estimator 

mata: bsp("`y'","`x'", "`v'",`vsemtd',`rnd',`upp',"`touse'")

// make ereturn list and display output 

tempname theta se V omega ai_selection bounds cl grid

mat `se' = r(se)
mat `theta' = r(theta)
mat `omega' = r(omega)
mat `bounds' = r(bounds)
mat `cl' = r(cvls)
mat `grid' = r(grid)

ereturn clear

local N = r(N)

ereturn local cmd = "clr_p"

if "`lower'" != "lower" {
	ereturn local title  = "CLR Intersection Upper Bounds (Parametric)"
}
else {
	ereturn local title  = "CLR Intersection Lower Bounds (Parametric)"	
}

if "`ais'" != "noais" {
	mat `ai_selection' = r(ais) 
}

ereturn local level = "`level'"
ereturn local depvar = "`y'"

ereturn scalar N = r(N)
ereturn scalar n_ineq = `numeq'
ereturn matrix omega = `omega'

display as text _newline e(title) _col(59) "Number of obs : " as result e(N)


// display grid and independent variables for each inequality 

local grid_count = 0
tokenize `x'

forval i = 1/`numeq' {
	
	tempname ais`i' se`i' theta`i'
	ereturn scalar grid`i' = `grid'[`i',1]
	
	local num_indep`i' = indep_vector[`i',1]
	forval j = 1/`num_indep`i'' {
		local indeplist`i' `indeplist`i'' `1' 
		macro shift
	}
	
	ereturn local indep`i' "`indeplist`i''"
	
	if "`ais'" != "noais" {
		mat `ais`i'' = `ai_selection'[`grid_count'+1..`grid_count'+e(grid`i'),1]
		ereturn matrix ais`i' = `ais`i''
	}
	mat `se`i'' = `se'[`grid_count'+1..`grid_count'+e(grid`i'),1]
	mat `theta`i'' = `theta'[`grid_count'+1..`grid_count'+e(grid`i'),1]
	
	ereturn matrix se`i' = `se`i''
	ereturn matrix theta`i' = `theta`i''

	display as text "Inequality #`i' : " word(e(depvar),`i') " (# of Grid Points : " as result e(grid`i') as text ", Independent Variables : " e(indep`i') " )"
	local grid_count = `grid_count' + e(grid`i')
}

tokenize `v'
forval i = 1/`numeq' {
	forval j = 1/`num_indep`i'' {
		local range`i' `range`i'' `1' 
		macro shift
	}
	ereturn local range`i' "`range`i''"
}
	
if "`ais'" != "noais" {
	display as text _newline "AIS(adaptive inequality selection) is applied" 
}
else { 
	display as text _newline "AIS(adaptive inequality selection) is not applied" 
}


// display result 
	
display as text _newline _col(38) "{c |}" _col(55) "Value"
display as text "{hline 37}{c +}{hline 45}"

forval i = 1/`nlevel' {
	local bd_level = level_vector[`i',1]
	local bd_name = 100 * `bd_level'
	while (int(`bd_name') != `bd_name') {
		local bd_name = `bd_name' * 10 
	}
	
	ereturn scalar bd`bd_name' = `bounds'[`i',1]
	ereturn scalar cl`bd_name' = `cl'[`i',1]

	if `bd_level' == 0.5 {
		display as text "half-median-unbiased est." _col(38)"{c |}" _col(53) as result %9.7f e(bd`bd_name')
	}
	else {
		if "`lower'" == "lower" {
			display as text 100*`bd_level' "% one-sided confidence interval" _col(38)"{c |}" _col(48) "[ " as result %9.7f e(bd`bd_name') as text ", inf)"
		}
		else {
			display as text 100*`bd_level' "% one-sided confidence interval" _col(38)"{c |}" _col(48) "( -inf, " as result %9.7f e(bd`bd_name') as text " ]"
		}
	}

}

display as text "{hline 37}{c BT}{hline 45}"

return clear
mat drop level_vector
mat drop indep_vector
	
end


version 11.2
mata: mata clear
mata:
void bsp(string vector y, string vector x, string scalar v, real scalar vsemtd, real scalar rnd_num, real scalar upp, string scalar touse){
	
	
	// Import Data from Stata
	
	Y = X = V =.
	Y = st_data(.,y,touse)
	X = st_data(.,x,touse)
	V = st_data(.,v)
	
	
	level = st_matrix("level_vector")
	indep = st_matrix("indep_vector")
	
	if (upp == 1) {
		Y = -Y
	}
	
	// Check missing values and rearrange x,y matrices
	
	nonmissingrow = . 
	
	y_missing = rowmissing(Y)
	x_missing = rowmissing(X)
	
	for (i_count = 1; i_count <= rows(Y); i_count++){
		if (y_missing[i_count] + x_missing[i_count] == 0) {
			nonmissingrow = nonmissingrow\i_count 
		}
	}
	
	X = X[nonmissingrow[2::rows(nonmissingrow)],.]
	Y = Y[nonmissingrow[2::rows(nonmissingrow)],.]
	
	nn = rows(Y)
	j = cols(Y)
	nindep = cols(X) 
	st_numscalar("r(N)",nn)
	
	// Construct large y_adj and x_adj 
	
	y_adj = vec(Y) 
	x_adj = J(nn*j,nindep+j,0)
	inter = J(nn,1,1)
	
	
	// caculate the size of v 
	
	grid_vector = J(j,1,0)
	
	acc_indep = 1
	x_colcount = 0 
	
	for(i_indep = 1; i_indep <=j; i_indep++){
		c_indep = indep[i_indep]
		grid_mat = V[.,x_colcount+1::x_colcount+c_indep] 
		grid_missing = rowmissing(grid_mat) 
		grid_count = 0 
		for (i_count = 1; i_count <= rows(grid_mat); i_count++) {
			if (grid_missing[i_count] == 0) {
				grid_count++
			}
		}
		grid_vector[i_indep] = grid_count 
	}
	
	
	
	st_matrix("r(grid)",grid_vector)
	v_adj = J(sum(grid_vector),nindep+j,0)
	
	acc_indep = 1
	x_colcount = 0 
	v_index = 0 
	
	// maka pseudo kronecker product x_adj, v_adj   
	
	
	for (i_indep = 1; i_indep <= j; i_indep++){
		c_indep = indep[i_indep]
		x_adj[nn*(i_indep-1)+1::nn*i_indep,acc_indep] = inter
		x_adj[nn*(i_indep-1)+1::nn*i_indep,acc_indep+1::acc_indep+c_indep] = X[.,x_colcount+1::x_colcount+c_indep] 
		
		grid_mat = V[.,x_colcount+1::x_colcount+c_indep] 
		grid_missing = rowmissing(grid_mat) 
		grid_row = . 
		for (i_count = 1; i_count <= rows(grid_mat); i_count++) {
			if (grid_missing[i_count] == 0) {
				grid_row = grid_row\i_count
			}
		}
		
		grid_row = grid_row[2::rows(grid_row)]
		
		v_input = grid_mat[grid_row,.]
		v_inter = J(rows(v_input),1,1)
		
		v_adj[v_index+1::v_index+grid_vector[i_indep],acc_indep] = v_inter
		v_adj[v_index+1::v_index+grid_vector[i_indep],acc_indep+1::acc_indep+c_indep] = v_input 
		
		v_index = v_index + grid_vector[i_indep]
		x_colcount = x_colcount + c_indep 
		acc_indep = acc_indep + c_indep + 1 
	}
	
	st_matrix("r(grid)",grid_vector)
	
	// Random Number Generating Process
	
	
	argminset = v_adj
	
	rnd_MAT = .
	rnd_MAT = rnd_R(rows(argminset),rnd_num)
	
	if (vsemtd == 1) {
		gamma_n = 1-0.1/log(nn)
		
		para_est(y_adj,x_adj,argminset,gamma_n,rnd_MAT,j,ml,sel,cvl)
		
		cutoff = colmax(ml - sel:*cvl) :- 2*sel:*cvl
		setindex = (ml :>= cutoff)
		argindex = .
		
		for (i = 1; i <= rows(argminset); i++) {
			if (setindex[i] > 0) {
				argindex = argindex \ i 
			}
		}
		
		argminset = argminset[argindex[2::rows(argindex)],.]
		
		if (upp == 1) {
			ml = -ml
		}
		
		st_matrix("r(ais)", setindex)
		st_matrix("r(se)", sel)
		st_matrix("r(theta)",ml)
	
	}

	para_est(y_adj,x_adj,argminset,level,rnd_MAT,j,ml2,sel2,cvl2)

	l_bj = J(rows(level),1,0)
	
	for (i = 1; i <= rows(level) ; i++) {
		l_bj[i] = colmax(ml2 :- sel2 :* cvl2[i])
	}	
	
	if (upp == 1) {
		l_bj = -l_bj
		ml2 = -ml2 
	}
	
	st_matrix("r(bounds)",l_bj)
	st_matrix("r(cvls)",cvl2) 
	
	
	if (vsemtd == 0){
		
		st_matrix("r(se)", sel2)
		st_matrix("r(theta)",ml2)
	}	
}
end

// begin of rnd_R.mata
version 11.2
mata: 
// rnd_R 1.0.0 Wooyoung Kim 8July2012
real matrix rnd_R(real scalar vseq,real scalar rnd_num){
	rnd_MAT = rnormal(vseq,rnd_num,0,1)
	return(rnd_MAT)
}
end

// begin of sieve_est.mata
version 11.2 
mata: 
// sieve_est 1.0.0 Wooyoung Kim 12July2012
void para_est(real matrix y, real matrix x, real matrix z, real matrix level, real matrix rnd_MAT, real scalar j,fe,se,cv){ 
// x = data 
// z = points at which a function is estimated
// kt_x = numbers of approximating functions
// level = level of the confidence set
// rnd_MAT = randomly generated normal dist. 
// j = number of inequalities
// Output : function estimates, their standard errors, uniform critical values 

	tmp1 = pinv(x' * x)
	tmp2 = x'
	bhat = tmp1 * (tmp2 * y)
	uhat = y - x * bhat
	tmp3 = tmp2 * diag(uhat :^ 2) * tmp2'
	var = tmp1 * tmp3 * tmp1 
	
	ghat = z * bhat
	varhat = z * var * z'
	sehat = sqrt(diagonal(varhat))
	
	Sigmahat = diag(1:/sehat)*varhat*diag(1:/sehat)
	
	rnd_RR = cholesky(Sigmahat + (0.00000001)*I(rows(Sigmahat)))* rnd_MAT[(1::rows(z)),.]
	tmp_RR = colmax(rnd_RR)'
	
	
	acv = mm_quantile(tmp_RR,1,level)
	
	fe = ghat
	se = sehat
	cv = acv
	
	st_matrix("r(omega)",var)

}  
end