*! version 1.0.0  11sep2015
*! version 1.0.1  5apr2016
*! version 1.0.2  8jul2016
*! version 1.0.3  29aug2016
*! version 1.0.4  15sep2016

program define teradialbc, eclass sortpreserve

	version 11
	local version : di "version " string(_caller()) ", missing :"

	syntax anything(equalok) [if] [in] [, Rts(string) Base(string) REFerence(string) ///
		Kappa(real 0.7) REPS(integer 999) Level(real 95) ///
		TEname(string) TEBC(string) BIASboot(string) VARboot(string) ///
		BIASSQVAR(string)	TELOWer(string) TEUPPer(string) ///
		noPRINT noDOTS SMOOTHED HETerogeneous SUBSampling]

	if "`print'" != "noprint" {
		local ls   = c(linesize)-7
		local lshl = c(linesize)-7-8
	}

	capture mata mata which kdens_bw()
	if _rc {
		display "{p 1 1 7 `ls'} {error}kdens_bw() not found; type -net install kdens.pkg- to obtain it{p_end}"
    exit 3499
	}
	capture mata mata which mm_quantile()
	if _rc {
		display "{p 1 1 7 `ls'} {error}mm_quantile() from -moremata- not found; type -ssc install moremata- to obtain it{p_end}"
   	exit 3499
	}

	* handle plugin
	local os "upper(substr("`c(os)'",1,3))"
	local vers = c(stata_version)
	if `vers' >= 12 {
		local bit = c(bit)
		*di `bit'
	}
	*di `vers'
	if `os' == "MAC"{
		capture program te_radial, plugin using ("te_radial.plugin")
	}
	else if `os' == "UNI"{
		local mach "upper(substr("`c(machine_type)'",1,3))"
		if `mach' == "MAC"{
			capture program te_radial, plugin using ("te_radial.plugin")
		}
		else {
			capture program te_radial, plugin using ("te_radial_ubuntu.plugin")
		}
	}
	else if `os' == "WIN"{
		if `vers' >= 12 {
			if `bit' == 64 {
				*di "64"
				capture program te_radial, plugin using ("te_radial_windows.plugin")
			}
			else {
				*di "32_1"
				capture program te_radial, plugin using ("te_radial_windows32.plugin")
			}
		}
		else {
			*di "32_2"
			capture program te_radial, plugin using ("te_radial_windows32.plugin")
		}
	}
	else {
		di "{p 1 1 7}{error}plugin to solve linear programming problem is not available for your system; contact developers{p_end}"
		exit 199
	}

	if `level' < 10 |  `level' > 99.99 {
		di "{p 1 1 7}{error}level() must be between 10 and 99.99 inclusive{p_end}"
		exit 198
	}

	if `kappa' < 0.5 |  `kappa' > 1 {
		di "{p 1 1 7}{error}kappa() must be between 0.5 and 1{p_end}"
		exit 198
	}

	scalar nrep = `reps'

	if `reps' < 100 {
		di "{p 1 1 7}{error}'reps' must be at least 100{p_end}"
		exit 498
	}

	if `reps' < 200 {
		di "{p 1 1 7}{error}Warning: Statistical inference may be unreliable for small number of bootstrap replications{p_end}"
	}

	if `reps' > 2000 {
		di "{p 1 1 7}{error}Warning: Unnecessary too many bootstrap replications; consider setting 'reps' smaller than 2000{p_end}"
		*exit 498
	}

	/*
		handle RTS
	*/

	local rts = substr("`rts'",1,1)
	if "`rts'" == "" | upper("`rts'") == "C" {
		scalar rt = 3
		local rt = "CRS"
	}
	else if upper("`rts'") == "N" {
			scalar rt = 2
			local rt = "NIRS"
		}
		else if upper("`rts'") == "V" {
				scalar rt = 1
				local rt = "VRS"
			}
			else {
				di "{p 1 1 7}{error}invalid returns to scale{p_end}"
				exit 198
			}

	/*
		handle base
	*/

	local base = substr("`base'",1,1)
	if "`base'" == "" | upper("`base'") == "O" {
		scalar ba = 2
		local ba = "output"
	}
	else if upper("`base'") == "I" {
			scalar ba = 1
			local ba = "input"
		}
		else {
			di "{p 1 1 7}{error}invalid base{p_end}"
			exit 198
		}

	/*
		check the names of variable(s) to be created
	*/

	if "`tename'" != ""{
		quietly generate `tename' = .
		local tevarname "Farrell `ba'-based TE measure under `rt'"
		label var `tename' "`tevarname'"
	}
	if "`tebc'" != ""{
		quietly generate `tebc' = .
		local tevarname "Bias-corrected Farrell `ba'-based TE measure under `rt'"
		label var `tebc' "`tevarname'"
	}
	if "`biasboot'" != ""{
		quietly generate `biasboot' = .
		local tevarname "Bootstrap bias estimate for Farrell `ba'-based TE measure under `rt'"
		label var `biasboot' "`tevarname'"
	}
	if "`varboot'" != ""{
		quietly generate `varboot' = .
		local tevarname "Bootstrap variance estimate for Farrell `ba'-based TE measure under `rt'"
		label var `varboot' "`tevarname'"
	}
	if "`biassqvar'" != ""{
		quietly generate `biassqvar' = .
		local tevarname "Bias-Variance ratio for Farrell `ba'-based TE measure under `rt'"
		label var `biassqvar' "`tevarname'"
	}
	if "`telower'" != ""{
		quietly generate `telower' = .
		local tevarname "Lower bound estimate for Farrell `ba'-based TE measure under `rt'"
		label var `telower' "`tevarname'"
	}
	if "`teupper'" != ""{
		quietly generate `teupper' = .
		local tevarname "Upper bound estimate for Farrell `ba'-based TE measure under `rt'"
		label var `teupper' "`tevarname'"
	}

	/*
		check boot options
	*/

	if "`smoothed'" == "smoothed" & "`subsampling'" == "subsampling" {
		di "{p 1 1 7}{error}options 'smoothed' and 'subsampling' cannot be specified together{p_end}"
		exit 198
	}
	if "`subsampling'" == "subsampling" {
		local smoothing = 0
		local homogeneous = 0
	}
	else {
		local smoothing = 1
		if "`heterogeneous'" == "heterogeneous" {
			local homogeneous = 0
		}
		else {
			local homogeneous = 1
		}
	}

	_deaparse       `anything'	// provides s(outps), s(inps), s(outpsRef), s(inpsRef)
	local outps     `s(outps)'
	local inps      `s(inps)'
	local outpsRef  `s(outpsRef)'
	local inpsRef   `s(inpsRef)'
	local refThere  `s(Ref)'

	/*
		check if all numeric
	*/

	confirm numeric variable `outps'
	confirm numeric variable `inps'
	confirm numeric variable `outpsRef'
	confirm numeric variable `inpsRef'

	* check is not constant and positive
	foreach var of varlist `outps' `inps' {
		quietly su `var'
		if r(sd) < 1e-6 {
			capture noi error 409
			di "{p 1 1 7}{error}Variable `var' has no variance{p_end}"
			exit 409
		}
		if r(min) <= 0 {
			capture noi error 411
			di "{p 1 1 7}{error}Variable `var' has nonpositive values{p_end}"
			exit 411
		}
	}
	if "`reference'" != ""{
		foreach var of varlist `outpsRef' `inpsRef' {
			quietly su `var'
			if r(sd) < 1e-6 {
				capture noi error 409
				di "{p 1 1 7}{error}Variable `var' has no variance{p_end}"
				exit 409
			}
			if r(min) <= 0 {
				capture noi error 411
				di "{p 1 1 7}{error}Variable `var' has nonpositive values{p_end}"
				exit 411
			}
		}
	}

	marksample touse
	markout `touse' `outps' `inps'

	scalar nO    = wordcount("`outps'")
	scalar nI    = wordcount("`inps'")
	scalar nORef = wordcount("`outpsRef'")
	scalar nIRef = wordcount("`inpsRef'")

	if nO != nORef{
		di "{p 1 1 7}{error}Number of outputs in the data and reference do not match{p_end}"
		exit 198
	}
	if nI != nIRef{
		di "{p 1 1 7}{error}Number of inputs in the data and reference do not match{p_end}"
		exit 198
	}


	capture matrix drop yobs
	capture matrix drop xobs
	capture matrix drop yxref
	capture matrix drop yref
	capture matrix drop xref

	local yxRef `outpsRef' `inpsRef'

	mkmat `outps' if `touse', matrix(yobs)  nomissing
	mkmat `inps'  if `touse', matrix(xobs)  nomissing
	if "`reference'" == "" {
		mkmat `yxRef' if `touse', matrix(yxref) nomissing
	}
	else {
		mkmat `yxRef' if `reference', matrix(yxref) nomissing
	}

	mat yref = yxref[1..., 1..nO]
	mat xref = yxref[1..., nO+1..nO+nI]

	scalar nobs = rowsof(yobs)
	scalar nref = rowsof(yref)

	if nref <= nO + nI {
		capture noi error 409
		display
		if "`reference'" == "" & `refThere'  == 0 {
			display "{p 1 1 2 `ls'} {error}Number of {error}reference data points (" as result nref "{error}) must be larger than sum of number of outputs (" as result nO "{error}) and inputs(" as result nI "{error}){p_end}"
		}
		else {
			display "{p 1 1 2 `ls'} {error}Number of {input} provided {error}reference data points (" as result nref "{error}) must be larger than sum of number of outputs (" as result nO "{error}) and inputs(" as result nI "{error}){p_end}"
		}
		matrix drop yobs xobs yxref yref xref
		scalar drop nO nI nORef nIRef nref nobs rt ba
		exit 409
	}

	sca ifqh = 1

	matrix te = J(nobs, 1, .)
	if "`tename'" != ""{
		scalar intovar = 1
		plugin call te_radial `tename' if `touse', yobs xobs nO nI nobs yref xref nref rt ba ifqh intovar te
	}
	else {
		scalar intovar = 0
		plugin call te_radial, yobs xobs nO nI nobs yref xref nref rt ba ifqh intovar te
	}

	/*
		handle infeasible
	*/

	mata : chkMiss("te", "touseNoMiss")

	/*
		modify the existing matrices by tossing the missing values
	*/

	mata : selReplace("yobs", "touseNoMiss")
	mata : selReplace("xobs", "touseNoMiss")
	mata : selReplace("te", "touseNoMiss")

	scalar nobs = rowsof(yobs)

	if nobs == 0 {
		di "{p 1 1 7}{error}Could not compute measure of technical efficiency for a single data point{p_end}"
		exit 498
	}

	/*
		redefine `touse'
	*/

	local my2use = "my2use"
	quietly generate `my2use' = .
	mata : st_store(., st_local("my2use"), st_local("touse"), st_matrix("touseNoMiss"))
	quietly replace `touse' = 0 if `my2use' == 0
	quietly drop `my2use'


	if "`print'" != "noprint" {
		display
		display "{p 1 1 7 `ls'} {text}Radial (Debreu-Farrell) {input:`ba'}{text}-based measures of technical efficiency under assumption of {input:`rt'} {text} technology are computed for the following data:{p_end}"
		display
		display as txt "    Number of data points ({input:K}{text}) = " as result nobs
		display as txt "    Number of outputs     ({input:M}{text}) = " as result nO
		display as txt "    Number of inputs      ({input:N}{text}) = " as result nI
		display
		if "`reference'" == ""  & `refThere'  == 0{
			display "{p 1 1 7 `ls'} {text}Data for reference set are not provided. Reference set is formed by " as result nobs "{text} data points, for which measures of technical efficiency are computed.{p_end}"
		}
		else {
			display "{p 1 1 7 `ls'} {text}Reference set is formed by " as result nref "{input} provided {text} reference data points.{p_end}"
		}
		display "{p 1 1 7 `ls'}{text}{hline `lshl'}{p_end}"
	}

	/*
		Begin Bootstrap
	*/

	/*
		Put matrices into mata
		(not to feed them to mata every time in resampling)
	*/

	/*
		Data
	*/

	mata : nobs = st_numscalar("nobs")
	mata : nO   = st_numscalar("nO")
	mata : nI   = st_numscalar("nI")
	mata : nref = st_numscalar("nref")
	mata : yref = st_matrix("yref")
	mata : xref = st_matrix("xref")
	mata : te   = st_matrix("te")

	sca smoothed = `smoothing'
	mata : smoothed = st_numscalar("smoothed")
	sca level = (100-`level')/100
	mata : level = st_numscalar("level")
	*mata : te  = st_matrix("te")

	if `smoothing' {
		/*
			Needed for homogeneous and heterogeneous
		*/

		mata : ba = st_numscalar("ba")

		/*
			Calculate TE for reference data points
		*/

		if "`reference'" == "" {
			mata : teRef = st_matrix("te")
		}
		else {
			matrix teRef = J(nref, 1, .)
			plugin call te_radial, yref xref nO nI nref yref xref nref rt ba ifqh intovar teRef
			* handle infeasible
			mata : chkMiss("teRef", "touseNoMissRef")
			* modify the existing matrices by tossing the missing values
			mata : selReplace("yref", "touseNoMissRef")
			mata : selReplace("xref", "touseNoMissRef")
			mata : selReplace("teRef", "touseNoMissRef")
			scalar nref = rowsof(yref)
			* feed feasible to mata
			mata : teRef = st_matrix("teRef")
		}
		if `homogeneous' {
			/*
				Homogeneous:
				- te and te reflected
				- bandwidth
				- scalar
			*/

			mata : terfl    = teRef \ 2:- teRef
			mata : mybw     = kdens_bw( terfl, 1, "sjpi")
			mata : scVarHom = ( 1 + mybw^2 / variance(terfl) )^(-1/2)
			/*
				end homogeneous
			*/
		}
		else {
			/*
				begin heterogeneous
			*/

			/*
				translation to polar coordinates
			*/

			if (ba == 1) {
				/*
					- angle
					- Z
				*/

				if (nI == 1) {
					mata : Z1 = yref, xref
				}
				else {
					mata : nu = atan(xref[,2..nI] :/ xref[,1])
					mata : nu = editmissing(nu, 0)
					mata : Z1 = yref, nu
					/* modulus */
					*mata: wx = sqrt( nI * mean((xref'):^2) )
				}
			}
			* output based
			else {
				/*
					- angle
					- Z
				*/

				if (nO == 1) {
					mata : Z1 = yref, xref
				}
				else {
					mata : nu = atan(yref[,2.. nO] :/ yref[,1])
					mata : nu = editmissing(nu, 0)
					mata : Z1 = nu, xref
					/* modulus */
					*mata : wy = sqrt( nO * mean((yref'):^2) )
				}
			}
			mata : nZ    = cols(Z1) + 1
			mata : bwZ("bws")
			mata : bws   = st_matrix("bws")
			mata : onesN = J(nref, 1, 1)
			mata : M1    = I(nref) - J(nref, nref, 1/nref)

			/*
				- Augmented Z
				- Variances and Cholesky
				- Bandwidth (univariate)
			*/

			mata : Z        = Z1, teRef
			mata : ZR       = Z1, 2 :- teRef
			mata : Zt       = Z \ ZR
			mata : L1       = cholesky(variance(Z))'
			mata : L2       = cholesky(variance(ZR))'
			mata : BW       = bws, kdens_bw(teRef \ 2:- teRef, 1, "sjpi")
			mata : scVarHet = (1 :+ BW:^2):^(-1/2)
			mata : cons1    = onesN # scVarHet
			mata : cons2    = onesN # BW
			mata : mata drop Z ZR BW scVarHet
		}
		/*
			end heterogeneous
		*/
	}
	else {
		* begin subsampling
		scalar msub  = round(nref^`kappa')
		mata:  msub  = st_numscalar("msub")
		* end subsampling
	}

	capture matrix drop teboot
	capture matrix drop testar

	if "`print'" != "noprint" {
		display
		if "`reference'" == "" {
			display "{p 1 1 7 `ls'} {text}Bootstrapping reference set formed by " as result nobs "{text} data points and computing radial (Debreu-Farrell) {input: `ba'}{text}-based measures of technical efficiency under assumption of {input:`rt'}{text} technology for each of " as result nobs "{text} data points relative to the bootstrapped reference set {text} {p_end}"
			display
		}
		else {
			display "{p 1 1 7 `ls'} {text}Bootstrapping reference set formed by " as result nref "{input} provided {text} reference data points and computing radial (Debreu-Farrell) {input: `ba'}{text}-based measures of technical efficiency under assumption of {input:`rt'}{text} technology for each of " as result nobs "{text} data points relative to the bootstrapped reference set {text} {p_end}"
			display
		}
		if `smoothing'{
			if `homogeneous' {
				display as txt " Smoothed " as input "homogeneous" as txt " bootstrap (" ///
								as res nrep as txt " replications)"
			}
			else {
				display as txt " Smoothed " as input "heterogeneous" as txt " bootstrap (" ///
								as res nrep as txt " replications)"
			}
		}
		else {
			display as input " Subsampling" as txt " bootstrap (" ///
								as res nrep as txt " replications)"
		}
		display
		if "`dots'" != "nodots" {
			_dots 0
		}

	}

	scalar intovar = 0

	matrix teB = J(nobs, 1, .)
	mata: BOOT = J(st_numscalar("nrep"),st_numscalar("nobs"),.)

	forvalues num = 1/`reps'{
		sca ifqh = 0
		sca numsc = `num'
		if `smoothing' {
			* begin smoothed bootstrap
			if `homogeneous' {
				* begin homogeneous
				if ba == 1 {
					* step 1: sampling
					mata: smplHom("xstar")
					*mat list xstar1
					* step 2: applying DEA
					plugin call te_radial, yobs xobs nO nI nobs yref xstar nref rt ba ifqh intovar teB
				}
				else {
					* step 1: sampling
					mata: smplHom("ystar")
					* step 2: applying DEA
					plugin call te_radial, yobs xobs nO nI nobs ystar xref nref rt ba ifqh intovar teB
				}
				* end homogeneous
			}
			else {
				* begin heterogeneous
				* "smplHet" and "refHet" take care of base
				* step 1: sampling
				mata: smplHet("ystar", "xstar", "teStar")
				matrix teRefB = J(nref, 1, .)
				* step 2: get efficiency under H0
				plugin call te_radial, ystar xstar nO nI nref yref xref nref rt ba ifqh intovar teRefB
				* step 3: get efficient xRef or yRef
				if ba == 1 {
					mata: refHet("xstar", "teRefB", "teStar")
				}
				else {
					mata: refHet("ystar", "teRefB", "teStar")
				}
				* handle infeasible
				mata : chkMiss("teRefB", "touseNoMissRefB")
				* modify the existing matrices by tossing the missing values
				mata : selReplace("ystar", "touseNoMissRefB")
				mata : selReplace("xstar", "touseNoMissRefB")
				mata : selReplace("teRefB", "touseNoMissRefB")
				* new number of obs in reference
				scalar nRef = rowsof(ystar)
				* step 4: get the bootstrapped TE
				plugin call te_radial, yobs xobs nO nI nobs ystar xstar nRef rt ba ifqh intovar teB
				* end heterogeneous
			}
			* end smoothed bootstrap
		}
		else {
			* begin subsampling bootstrap
			* step 1: sub-sampling
			mata : smplSub("ystar", "xstar")
			* step 2: applying DEA
			plugin call te_radial, yobs xobs nO nI nobs ystar xstar msub rt ba ifqh intovar teB
			* end subsampling bootstrap
		}
		/*
			posting the results from any type of bootstrap
		*/
		mata: BOOT[st_numscalar("numsc"),] = st_matrix("teB")'
		*mat teboot = (nullmat(teboot) \ teB')
		if "`print'" != "noprint" {
			if "`dots'" != "nodots" {
				_dots `num' 0
			}
		}
	}

	if "`print'" != "noprint" {
		if mod(`reps',50) != 0 {
			display
		}
	}

	/*
		Working with large BOOT matrix
	*/

	/*
		Step 3: bias correction and CIs
	*/

	mata : forceLargerOne = 0
	mata : biasAndCI("realnrep","bias","bootvar","biasOverVar","tebc","telow","teupp")

	/*
		Check if any bias corrected measure is negative;
		If yes, do inference in terms of Shephard
		distance functions (must be input oriented)
	*/
	mata: teBCneg("tebc", "tebcneg")
	if tebcneg == 1 {
		if "`print'" != "noprint" {
			display
			display "{p 1 1 7 `ls'} {result}One or more bias-corrected radial (Debreu-Farrell) {input: `ba'}{result}-based measures of technical efficiency is negative.  Analysis will be done in terms of Shephard distance function, a reciprocal of the Debreu-Farrell measure{p_end}"
		}
		mata : forceLargerOne = 1
		mata : biasAndCI("realnrep","bias","bootvar","biasOverVar","tebc","telow","teupp")
		mata : oneOver("te")
		if "`tename'" != ""{
			quietly replace `tename' = 1/`tename'
		}
	}

	/*
		Check if bias squared over var is larger for all
	*/
	mata: biasSqVarOne("biasOverVar", "bSqVar")
	if bSqVar == 1 {
		if "`print'" != "noprint" {
			display "{p 1 1 7 `ls'}{text}{hline `lshl'}{p_end}"
			display
			display "{p 1 1 7 `ls'} {result}Warning: for one or more data points statistic [(1/3)*(bias)^2/var] is smaller than 1; bias-correction should not be used{p_end}"
		}
	}


	* write to variables if names are provided
	if "`tebc'" != ""{
		mata: st_store(., st_local("tebc"), st_local("touse"), st_matrix("tebc"))
*		plugin call vectovar `tebc' if `touse', tebc
	}
	if "`biasboot'" != ""{
		mata: st_store(., st_local("biasboot"), st_local("touse"), st_matrix("bias"))
	}
	if "`varboot'" != ""{
		mata: st_store(., st_local("varboot"), st_local("touse"), st_matrix("bootvar"))
	}
	if "`biassqvar'" != ""{
		mata: st_store(., st_local("biassqvar"), st_local("touse"), st_matrix("biasOverVar"))
	}
	if "`telower'" != ""{
		mata: st_store(., st_local("telower"), st_local("touse"), st_matrix("telow"))
	}
	if "`teupper'" != ""{
		mata: st_store(., st_local("teupper"), st_local("touse"), st_matrix("teupp"))
	}

	* return

	* macros
	eret post, esample(`touse')

	ereturn local ref_inputs `"`inpsRef'"'
	ereturn local ref_outputs `"`outpsRef'"'
	ereturn local inputs `"`inps'"'
	ereturn local outputs `"`outps'"'
	ereturn local base `"`ba'"'
	ereturn local rts `"`rt'"'
	ereturn local cmdline `"`0'"'
	ereturn local cmd "teradialbc"
	ereturn local title "Bias-correction of `ba'-based measures of technical efficiency"

	* scalars
	ereturn scalar K = nobs
	ereturn scalar M = nO
	ereturn scalar N = nI
	ereturn scalar Kref = nref
	ereturn scalar reps = `reps'

	* matrices
	*ereturn matrix teboot = teboot
	ereturn matrix teupp = teupp
	ereturn matrix telow = telow
	ereturn matrix realreps = realnrep
	ereturn matrix biassqvar = biasOverVar
	ereturn matrix varboot = bootvar
	ereturn matrix biasboot = bias
	ereturn matrix tebc = tebc
	ereturn matrix te = te

	/*
		clean up
	*/

	if `smoothing' {
		if `homogeneous' {
			if ba == 1 {
				matrix drop xstar
			}
			else {
				matrix drop ystar
			}
			mata mata drop mybw scVarHom terfl
		}
		else {
			capture mata mata drop nu
			mata mata drop L1 L2 M1 onesN Z1 Zt nZ bws cons1 cons2
			matrix drop teRefB xstar ystar touseNoMissRefB teStar bws
			scalar drop nRef
		}
		capture mata mata drop msub
		mata mata drop ba teRef
		capture matrix drop touseNoMissRef
		capture matrix drop teRef
	}
	else {
		scalar drop msub
		matrix drop xstar ystar
		mata mata drop msub
	}

	* common
	mata mata drop nobs nO nI nref xref yref te level smoothed ///
							forceLargerOne BOOT
	scalar drop tebcneg level smoothed intovar ifqh ///
							bSqVar nref nobs nIRef nORef nI nO ba rt numsc nrep
	matrix drop teB xref yref xobs yobs touseNoMiss yxref

end