*********************************************************************************************************************************
*********************************************************************************************************************************
*     						ACREG: Program for computing corrected standard errors for Arbitrary Clustering		   				*
*		   		 					 Copyright: F. Colella, R. Lalive, S.O. Sakalli, M. Thoenig								    *
*																																*
*									Beta Version, please do not circulate - This Version: December 2020							*
*																																*
* 									Before using this program, please read to our companion papers								*
*																																*
*  			"Acreg: arbirtrary correlation regression", F. Colella, R. Lalive, S.O. Sakalli, M. Thoenig (2020)					*
*																																*
*  	"Inference with Arbitrary Clustering", F. Colella, R. Lalive, S.O. Sakalli, M. Thoenig, (2019) IZA Discussion Papaer		*
*																																*
*********************************************************************************************************************************
*********************************************************************************************************************************

capt program drop acreg

*! Version December 2020  (1.1.0)
*! ACREG: Arbitrary Correlation Regression
*! Authors: F. Colella, R. Lalive, S.O. Sakalli, M. Thoenig


program acreg, eclass
	version 12
	if replay() {
	di "No Variables"
	exit
	}
	else {
syntax [anything(name=0)]  [if] [in] [fweight pweight] [,  dist_mat(varlist) links_mat(varlist) ///
	weights(varlist) cluster(varlist)  ///
	latitude(varname) longitude(varname) id(varname) time(varname) ///
	LAGcutoff(integer -1) DISTcutoff(real -1) LAGDISTcutoff(integer -1) ///
	bartlett partial correction network spatial storeweights storedistances small ak0 ///
	pfe1(varname) pfe2(varname) hac correctr2 nbclust(integer -1) dropsingletons ]

	tempname b V weightsmat iddd_var

	local current_dir `c(pwd)'

			local n 0

		gettoken depvar1 0 : 0, parse(" ,[") match(paren)
		PEnd `depvar1'
		if `s(stop)' {
			error 198
		}
		while `s(stop)'==0 {
			if "`paren'"=="(" {
				local n = `n' + 1
				if `n'>1 {
capture noi error 198
di in red `"syntax is "(all instrumented variables = instrument variables)""'
exit 198
				}
				gettoken p depvar1 : depvar1, parse(" =")
				while "`p'"!="=" {
					if "`p'"=="" {
capture noi error 198
di in red `"syntax is "(all instrumented variables = instrument variables)""'
di in red `"the equal sign "=" is required"'
exit 198
					}
					local end1 `end1' `p'
					gettoken p depvar1 : depvar1, parse(" =")
				}
				local temp_ct  : word count `end1'
				if `temp_ct' > 0 {
	//				tsunab end1 : `end1'
					prog_fv_unab `end1'
					local end1 `r(ts_varlist)'
				}
* Allow for empty instrument list
				local temp_ct  : word count `depvar1'
				if `temp_ct' > 0 {
	//				tsunab ivv1 : `depvar1'
					prog_fv_unab `depvar1'
					local ivv1 `r(ts_varlist)'
				}
			}
			else {
				local exog1 `exog1' `depvar1'
			}
			gettoken depvar1 0 : 0, parse(" ,[") match(paren)
			PEnd `depvar1'
		}
		local 0 `"`depvar1' `0'"'

	//	tsunab exog1 : `exog1'
		prog_fv_unab `exog1'
		local exog1 `r(ts_varlist)'
		tokenize `exog1'
		local depvar1 "`1'"
		local 1 " "
		local exog1 `*'

	capture drop touse
	marksample touse

}


/*
*** CONSIDER THE POSSIBILITY TO USE XTSET: do we need FE?
xtset
local iddd = r(panelvar)
di "`iddd'"
*/

* drop repeated variables for each fo the the four groups
if "`exog1'"=="" {
}
else {
	local dupvars : list dups exog1
}
if "`dupvars'"=="" {
}
else {
	local exog1 = subinstr("`exog1'","`dupvars' ","",1)
}

if "`depvar1'"=="" {
}
else {
	local dupvars : list dups depvar1
}
if "`dupvars'"=="" {
}
else {
	local depvar1 = subinstr("`depvar1'","`dupvars' ","",1)
}

if "`end1'"=="" {
}
else {
	local dupvars : list dups end1
}
if "`dupvars'"=="" {
}
else {
	local end1 = subinstr("`end1'","`dupvars' ","",1)
}

if "`ivv1'"=="" {
}
else {
	local dupvars : list dups ivv1
}
if "`dupvars'"=="" {
}
else {
	local ivv1 = subinstr("`ivv1'","`dupvars' ","",1)
}

* check for double inputs
	local allvars `exog1' `depvar1' `end1' `ivv1'
	local dupvars : list dups allvars
	verify "`dupvars'"=="", msg("Error: there are repeated variables: <`dupvars'>")


* check  other variables are non string
	local othervars `dist_mat' `links_mat' `weights' `cluster' `latitude' `longitude'
	check_number `othervars'

* check for missing values
	check_missing `othervars'

* check lat and lon are between -180 and 180
	check_lat_lon `latitude' `longitude'

* check links mat is made by 1 and 0
	check_links `links_mat'

* count number of factor variables
	prog_count_factor_variables  `depvar1' `exog1' `ivv1'
	mata: factor_variables = st_numscalar("factor_variables")

* check for ranktest installed
	capture which ranktest
	if (_rc==111) {
		disp as error  "please install raktest package by typing " in gr "ssc install ranktest"
		 disp as error  "OR install all required packages by typing " ///
		 in gr "acregpackcheck "
		 disp as error  "and run acreg again"
		exit 498
	 }



if (`nbclust' < 0){
	local nbclust = 100
}
else {
	local nbclust = `nbclust'
}


******
***create signals for cutoffs
if (`distcutoff' < 0){
	local distcutoff_yn = 0
	local distcutoff = 0
}
else {
	local distcutoff_yn = 1
}


if (`lagcutoff' < 0){
	local lagcutoff_yn = 0
	local lagcutoff = 0
}
else {
	local lagcutoff_yn = 1
}

if (`lagdistcutoff' < 0){
	local lagdistcutoff_yn = 0
	local lagdistcutoff = 0
}
else {
	local lagdistcutoff_yn = 1
}

if ("`time'" == ""){
	local time_yn = 0
}
else {
	local time_yn = 1
}






if ("`id'" == ""){
	if (`time_yn' == 0){
		capt drop `iddd_var'
		tempvar iddd_var
		gen `iddd_var' = _n
		local id `iddd_var'
		if (`lagcutoff_yn' == 1){
			disp as error "Lagcutoff requires Id and Time to be specified"
			exit 498
		}
		if (`lagdistcutoff_yn' == 1){
			disp as error "Lagdistcutoff requires Id and Time to be specified"
			exit 498
		}
		if ("`hac'" == ""){
		}
		else {
			disp as error "Hac option requires Id and Time to be specified"
			exit 498
		}
	}
	else {
		disp as error "ID is required if Time dimension is specified"
		exit 498
	}
}
else {
	if (`time_yn' == 0){
		disp as error "Time dimesion is required if id is specified"
		exit 498
	}
	else {
	//	xtset `id' `time'
		if ("`hac'" == ""){
			if (`lagcutoff_yn' == 0){
				*local lagcutoff=10000000000
				local lagcutoff=0
			}
			if (`lagdistcutoff_yn' == 0){
				*local lagdistcutoff=10000000000
				local lagdistcutoff=0
			}
		}
		else {
			if (`lagcutoff_yn' == 0){
				disp as error "Lagcutoff required with the HAC option"
				exit 498
			}
		}
	}
}



if ("`time'" == ""){
	qui sort `id'
}
else {
	qui sort `id' `time'
}


if ("`cluster'" == ""){
	if ("`weights'" == ""){
		if ("`network'" == ""){
			if ("`spatial'" == ""){
//1) NO CLUSTERING
				if ("`links_mat'" == ""){
				}
				else {
					disp as error "Matrix of links requires Network Option to be specified"
					exit 498
				}
				if ("`dist_mat'" == ""){
				}
				else {
					disp as error "Matrix of Distances requires Network or Spatial Option"
					exit 498
				}
				if ("`latitude'" == ""){
				}
				else {
					disp as error "Latitude and Longitude require Spatial Option to be specified"
					exit 498
				}
				if ("`longitude'" == ""){
				}
				else {
					disp as error "Latitude and Longitude require Spatial Option to be specified"
					exit 498
				}
				if (`distcutoff_yn' == 1){
					disp as error "Distance cutoff requires Network or Spatial Option"
					exit 498
				}
				if (`lagdistcutoff_yn' == 1){
					disp as error "Lag Distance cutoff requires Network or Spatial Option"
					exit 498
				}
				if ("`storedistances'" == "") {
				}
				else {
				disp as error "Distances can be stored only with Network or Spatial Option"
				exit 498
				}
				if (`lagcutoff_yn' == 1){
					di in gr "TEMPORAL CORRECTION"
				}
				else {
				di in gr "HETEROSKEDASTICITY ROBUST STANDARD ERRORS"
				}
				if ("`hac'" == ""){
					di in gr "No HAC Correction"
				}
				else {
					di in gr "HAC Correction "
				}
			}
			else {
//2) SPATIAL
				di in gr "SPATIAL CORRECTION"
				if ("`links_mat'" == ""){
				}
				else {
					disp as error "Matrix of Links may not be specified with the Spatial Option"
					exit 498
				}
				if ("`dist_mat'" == ""){
					if ("`latitude'" == ""){
						disp as error "Latitude and Longitude are required with the Spatial Option if Distance Matrix is not specified"
						exit 498
					}
					if ("`longitude'" == ""){
						disp as error "Latitude and Longitude are required with the Spatial Option if Distance Matrix is not specified"
						exit 498
					}
				}
				else {
					if ("`storedistances'" == "") {
					}
					else {
						disp as error "Distances can not be stored if they are given as input"
						exit 498
					}
					if ("`latitude'" == ""){
					}
					else {
						disp as error "Latitude may not be specified with the Spatial Option if Distance Matrix is specified"
						exit 498
					}
					if ("`longitude'" == ""){
					}
					else {
						disp as error "Longitude may not be specified with the Spatial Option  if Distance Matrix is specified"
						exit 498
					}
				}
				if (`distcutoff_yn' == 0){
					disp as error " Distance cutoff is required with the Spatial Option. Distcutoff must be a nonnegative real number"
					exit 498
				}
				di in gr "DistCutoff: " in ye "`distcutoff'"
				di in gr "LagCutoff:  " in ye "`lagcutoff'"
			//	di in gr "LagDistCutoff:  " in ye "`lagdistcutoff'"
				if ("`hac'" == ""){
					di in gr "No HAC Correction"
				}
				else {
				di in gr "HAC Correction"
		*			if ("`lagdistcutoff_yn'" == "0"){
		*			disp as error "Lag Distance cutoff required with HAC Option if Spatial is specified"
		*			exit 498
		*			}
				}
			}
		}
		else {
			if ("`spatial'" == ""){
//3) NETWORK
				di in gr  "NETWORK CORRECTION"
				if ("`latitude'" == ""){
				}
				else {
					disp as error "Latitude may not be specified with the Network Option"
					exit 498
				}
				if ("`longitude'" == ""){
				}
				else {
					disp as error "Longitude may not be specified with the Network Option"
					exit 498
				}
				if ("`dist_mat'" == ""){
					if ("`links_mat'" == ""){
						disp as error "Please specify the Links Matrix or the Distance Matrix when the Network Option is specified"
						exit 498
					}
					else {
						if (`distcutoff_yn' == 0){
							disp  "WARNING: Distance cutoff not Specified, only first degree links will be considered"
							local distcutoff=1
						}
					}
				}
				else {
					if ("`links_mat'" == ""){
						if (`distcutoff_yn' == 0){
							disp as error "Distance cutoff is required with the Network Option when Distance Matrix is specified. Distcutoff must be a nonnegative real number"
							exit 498
						}
					}
					else {
						disp as error "Distance Matrix may not be specified with the Network Option if Links Matrix is specified"
						exit 498
					}
					if ("`storedistances'" == "") {
					}
					else {
						disp as error "Distances can not be stored if they are given as input"
						exit 498
					}
				}
				di in gr "DistCutoff: " in ye "`distcutoff'"
				di in gr "LagCutoff:  " in ye "`lagcutoff'"
		//		di in gr "LagDistCutoff:  " in ye "`lagdistcutoff'"
				if ("`hac'" == ""){
					di in gr "No HAC Correction"
				}
				else {
					di in gr "HAC Correction "
		*			if ("`lagdistcutoff_yn'" == "0"){
		*			disp as error "Lag Distance cutoff required with HAC Option if Network is specified"
		*			exit 498
		*			}
				}
			}
			else {
				disp as error "Only one option between network and spatial is allowed"
				exit 498
			}
		}
	}
	else {
		if ("`spatial'" == "spatial"){
			disp as error "Spatial option with Weighting Matrix is not allowed"
			exit 498
		}
		if ("`network'" == "network"){
			disp as error "Network option with Weighting Matrix is not allowed"
			exit 498
		}
//4) WEIGHTS
		di in gr  "Correction using the Weighting Matrix  provided"
		if ("`time'" == ""){
		}
		else {
			disp as error "Panel Dimension not allowed with the Weighting Matrix""
			exit 498
		}
		if ("`latitude'" == ""){
		}
		else {
			disp as error "Latitude may not be specified with the Weighting Matrix"
			exit 498
		}
		if ("`longitude'" == ""){
		}
		else {
			disp as error "Longitude may not be specified with the Weighting Matrix"
			exit 498
		}
		if ("`links_mat'" == ""){
		}
		else {
			disp as error "Matrix of Links may not be specified with the Weighting Matrix"
			exit 498
		}
		if ("`dist_mat'" == ""){
		}
		else {
			disp as error "Distance Matrix may not be specified with the Weighting Matrix"
			exit 498
		}
		if (`distcutoff_yn' == 1){
			disp as error "Distcutoff may not be specified with the Weighting Matrix"
			exit 498
		}
		if (`lagcutoff_yn' == 1){
			disp as error "Lagcutoff may not be specified with the Weighting Matrix"
			exit 498
		}
		if ("`lagdistcutoff'" == ""){
			disp as error "Lagdistcutoff may not be specified with the Weighting Matrix"
			exit 498
		}
		if ("`hac'" == ""){
			di in gr "No HAC Correction"
		}
		else {
			di in gr "HAC Correction "
		}
	}
}
else {
//5) MULTIWAY CLUSTERING
	di in gr  "MULTIWAY CLUSTERING CORRECTION"
	di in gr "Cluster variable(s): " in ye "`cluster'"
	if ("`time'" == ""){
	}
	else {
		disp as error "Panel Settings not allowed with multiway clustering"
		exit 498
	}
	if ("`weights'" != ""){
		disp as error "Weights option with multiway clustering is not allowed"
		exit 498
	}
	if ("`spatial'" == "spatial"){
		disp as error "Spatial option with multiway clustering is not allowed"
		exit 498
	}
	if ("`network'" == "network"){
		disp as error "Network option with multiway clustering is not allowed"
		exit 498
	}
	if ("`latitude'" == ""){
	}
	else {
		disp as error "Latitude may not be specified with multiway clustering"
		exit 498
	}
	if ("`longitude'" == ""){
	}
	else {
		disp as error "Longitude may not be specified with multiway clustering"
		exit 498
	}
	if ("`links_mat'" == ""){
	}
	else {
		disp as error "Matrix of Links may not be specified with multiway clustering"
		exit 498
	}
	if ("`dist_mat'" == ""){
	}
	else {
		disp as error "Distance Matrix may not be specified with multiway clustering"
		exit 498
	}
		if (`distcutoff_yn' == 1){
			disp as error "Distcutoff may not be specified with multiway clustering"
			exit 498
		}
		if ("`lagdistcutoff'" == ""){
			disp as error "Lagdistcutoff may not be specified with multiway clustering"
			exit 498
		}
		if ("`hac'" == ""){
			di in gr "No HAC Correction"
		}
		else {
			di in gr "HAC Correction"
		}
}

	if ("`small'" == ""){
		}
		else {
		di in gr "Small Sample Correction"
		}
	if ("`ak0'" == ""){
		}
		else {
		di in gr "ak0 SS Correction - Conservative SEs"
		}
	if ("`time'" == ""){
		qui sort `id'
	}
	else {
	qui sort `id' `time'
	}
*




if ("`correctr2'" == ""){
	scalar correctr2_yn = 0
}
else {
	scalar correctr2_yn = 1

}

mata: correctr2_yn_s = st_numscalar("correctr2_yn")






********************************************************************************************


if ("`pfe1'"!=""){

preserve

 capture which hdfe
 if (_rc==111) {
 disp as error  "please install hdfe command by typing " in gr "ssc install hdfe"
 disp as error  "OR install all required packages by typing " ///
 in gr "acregpackcheck "
 disp as error  "and run acreg again"
 capt drop iddd_var
 exit 498
 }

 prog_fv_partial `depvar1' `end1' `exog1' `ivv1'

	if ("`pfe2'"!=""){
	di in gr "Absorbed FE:" in ye " `pfe1' " in gr "and" in ye " `pfe2'"
	local fe_1 "`pfe1'"
	local fe_2 "`pfe2'"
	}
	else {
	di in gr "Absorbed FE:" in ye " `pfe1' "
	local fe_1 "`pfe1'"
	local fe_2 "fake_fe2"
	tempvar  fake_fe2
	capt drop fake_fe2
	gen fake_fe2=0
	}


if ("`correctr2'" == ""){
}
else {
	capture which ivreg2
	if (_rc==111) {
		disp as error  "please install ivreg2 command by typing " in gr "ssc install ivreg2"
		disp as error  "OR install all required packages by typing " ///
		in gr "acregpackcheck "
		disp as error  "and run acreg again"
		capt drop iddd_var
		capt drop fake_fe2
		exit 498
	 }
*	preserve
		capt cd "`current_dir'"

	if "`weight'" == "" {
		if ("`pfe2'"==""){
			acreg_r2partial  `exog1' i.`pfe1' `if' `in' ,  depvar(`depvar1')  end(`end1') iv(`ivv1')
		}
		else {
			acreg_r2partial  `exog1' i.`pfe1' i.`pfe2' `if' `in' ,  depvar(`depvar1')  end(`end1') iv(`ivv1')
		}
	}
	else {
		if ("`pfe2'"==""){
			acreg_r2partial  `exog1' i.`pfe1'  `if' `in'  [`weight'`exp'],  depvar(`depvar1')  end(`end1') iv(`ivv1')
		}
		else {
			acreg_r2partial  `exog1' i.`pfe1' i.`pfe2' `if' `in'  [`weight'`exp'],  depvar(`depvar1')  end(`end1') iv(`ivv1')
		}
	}
*	restore
}

mata: correctr2_yn_s = st_numscalar("correctr2_yn")

if "`dropsingletons'" == "" {
if "`weight'" == "" {
qui hdfe `depvar1' `end1' `exog1' `ivv1' `if' `in' ,  keepsingletons absorb(`fe_1' `fe_2') generate(po_)
}
else {
qui hdfe `depvar1' `end1' `exog1' `ivv1' `if' `in' [`weight'`exp'] ,  keepsingletons  absorb(`fe_1' `fe_2') generate(po_)
}
}
else {
qui reg `depvar1' `end1' `exog1' `ivv1' `if' `in' [`weight'`exp']
qui count if e(sample)==1
local count1 =  r(N)
bys `fe_1':gen NN=_N
qui drop if NN==1
drop NN
if "`weight'" == "" {
qui hdfe `depvar1' `end1' `exog1' `ivv1' `if' `in' ,   absorb(`fe_1' `fe_2') generate(po_)
}
else {
qui hdfe `depvar1' `end1' `exog1' `ivv1' `if' `in' [`weight'`exp'] ,   absorb(`fe_1' `fe_2') generate(po_)
}
qui sum po_`depvar1'
local count2 =  r(N)
if `count2'==`count1' {
di in ye "No singleton observations"
}
else {
local differe = `count1' - `count2'
di in ye "             `differe' singleton observations dropped"
}
}


local count_var = 1
local temp_countvars  : word count `depvar1'
	if `temp_countvars'>0 {
	tokenize `depvar1'
	while  `count_var' <= `temp_countvars' {
	local depvar2  `depvar2' po_``count_var''
	local count_var=`count_var'+1
	}
	}

local count_var = 1
local temp_countvars  : word count `end1'
	if `temp_countvars'>0 {
	tokenize `end1'
	while  `count_var' <= `temp_countvars' {
	local end2 `end2' po_``count_var''
	local count_var=`count_var'+1
	}
	}

local count_var = 1
local temp_countvars  : word count `exog1'
	if `temp_countvars'>0 {
	tokenize `exog1'
	while  `count_var' <= `temp_countvars' {
	local exog2  `exog2' po_``count_var''
	local count_var=`count_var'+1
	}
	}

local count_var = 1
local temp_countvars  : word count `ivv1'
	if `temp_countvars'>0 {
	tokenize `ivv1'
	while  `count_var' <= `temp_countvars' {
	local ivv2 `ivv2' po_``count_var''
	local count_var=`count_var'+1
	}
	}


************************************************************
*** THIS WILL BE NEEDED TO IMPLEMENT THE TESTS ***

* Create group variable
//tempvar group
//qui makegps2, id1(`fe_1') id2(`fe_2') groupid(`group')

* Calculate Degrees of Freedom
//qui count
//local N = r(N)
//local k : word count `end1' `exog1' `ivv1' //Check whether here I need also the instruments or not
//sort `fe_1'
//qui count if `id'!=`id'[_n-1]
//local G1 = r(N)
//sort `fe_2'
//qui count if `time'!=`time'[_n-1]
//local G2 = r(N)
//sort `group'
//qui count if `group'!=`group'[_n-1]
//local M = r(N)
//local kk = `k' + `G1' + `G2' - `M'
//local dof = `N' - `kk'
//local G = `G2'-1
************************************************************


	drop `depvar1' `end1' `exog1' `ivv1'

	foreach var in `depvar2' `end2' `exog2' `ivv2' {
	local var2 = subinstr("`var'", "po_", "", .)
	rename `var' `var2'
	}



	capt cd "`current_dir'"

if "`weight'" == "" {
	acreg_core  `exog1' ,  depvar(`depvar1')  end(`end1') iv(`ivv1') ///
	dist_mat(`dist_mat') links_mat(`links_mat') weights(`weights') ///
	latitude(`latitude') longitude(`longitude') id(`id') time(`time')  ///
	lag(`lagcutoff') dist(`distcutoff') lagdist(`lagdistcutoff') ///
	cluster(`cluster') `hac' nbclust(`nbclust') ///
	`bartlett' `partial' `correction' `network' `spatial' `storeweights' `storedistances' `small' `ak0'
}
else {
	acreg_core  `exog1'  [`weight'`exp'],  depvar(`depvar1')  end(`end1') iv(`ivv1') ///
	dist_mat(`dist_mat') links_mat(`links_mat') weights(`weights') ///
	latitude(`latitude') longitude(`longitude') id(`id') time(`time')  ///
	lag(`lagcutoff') dist(`distcutoff') lagdist(`lagdistcutoff') ///
	cluster(`cluster') `hac' nbclust(`nbclust') ///
	`bartlett' `partial' `correction' `network' `spatial' `storeweights' `storedistances' `small' `ak0'
}

restore
}
else {
	if ("`pfe2'"!=""){
	disp as error "pfe2 requires pfe1." in ye " Please, specify pfe1 instead of pfe2 if you want to partial out only one fixed effect."
	exit 498
	}
}


mata {
if (correctr2_yn_s==0) {
RSSr2 =0
TSScr2 =0
TSSur2 =0
r2r2 =0
r2ur2 =0
}
st_numscalar("r(RSSr2)", RSSr2)
st_numscalar("r(TSScr2)", TSScr2)
st_numscalar("r(TSSur2)", TSSur2)
st_numscalar("r(r2r2)", r2r2)
st_numscalar("r(r2ur2)", r2ur2)
}


if ("`pfe1'"!=""){


if ("`correctr2'" == ""){

***output
di in gr _col(55) "Number of obs = " in ye %8.0f e(N)
*di in gr _c _col(55) "F(" %3.0f e(Fdf1) "," %6.0f e(Fdf2) ") = "
*		if e(F) < 99999 {
*di in ye %8.2f e(F)
*		}
*		else {
*di in ye %8.2e e(F)
*		}
*di in gr _col(55) "Prob > F      = " in ye %8.4f e(Fp)
di in gr "Total (centered) SS     = " in ye %12.0g e(tss) _continue
di in gr _col(55) "Centered R2   = " in ye %8.4f e(r2)
di in gr "Total (uncentered) SS   = " in ye %12.0g e(tssu) _continue
di in gr _col(55) "Uncentered R2 = " in ye %8.4f e(r2u)
di in gr "Residual SS             = " in ye %12.0g e(rss) // _continue
*di in gr _col(55) "Root MSE      = " in ye %8.4g e(rmse)
di

ereturn display

di in gr "nb: total SS, model and R2s are after partialling out."
di in gr "To get the corrected ones use the option correctr2"

}


else {

scalar 	MSSc = TSScr2 - RSSr2
scalar 	MSSu = TSSur2 - RSSr2

ereturn scalar rss=r(RSSr2)
ereturn scalar tss=r(TSScr2)
ereturn scalar tssu=r(TSSur2)
ereturn scalar r2=r(r2r2)
ereturn scalar r2u=r(r2ur2)
ereturn scalar mss=MSSc
ereturn scalar mssu=MSSu

***output
di in gr _col(55) "Number of obs = " in ye %8.0f e(N)
di in gr "Total (centered) SS     = " in ye %12.0g e(tss) _continue
di in gr _col(55) "Centered R2   = " in ye %8.4f e(r2)
di in gr "Total (uncentered) SS   = " in ye %12.0g e(tssu) _continue
di in gr _col(55) "Uncentered R2 = " in ye %8.4f e(r2u)
di in gr "Residual SS             = " in ye %12.0g e(rss) // _continue
*di in gr _col(55) "Root MSE      = " in ye %8.4g e(rmse)
di
ereturn display
}


}
else {
preserve

di in gr "No Absorbed FEs"

if "`weight'" == "" {
acreg_core  `exog1' `if' `in',  depvar(`depvar1') end(`end1') iv(`ivv1') ///
 dist_mat(`dist_mat') links_mat(`links_mat') weights(`weights') ///
	latitude(`latitude') longitude(`longitude') id(`id') time(`time')  ///
	lag(`lagcutoff') dist(`distcutoff') lagdist(`lagdistcutoff') ///
	cluster(`cluster')  `hac' nbclust(`nbclust') ///
	`bartlett' `partial' `correction' `network' `spatial' `storeweights' `storedistances' `small' `ak0'
}
else {
acreg_core  `exog1' `if' `in' [`weight'`exp'],  depvar(`depvar1') end(`end1') iv(`ivv1') ///
 dist_mat(`dist_mat') links_mat(`links_mat') weights(`weights') ///
	latitude(`latitude') longitude(`longitude') id(`id') time(`time')  ///
	lag(`lagcutoff') dist(`distcutoff') lagdist(`lagdistcutoff') ///
	cluster(`cluster')  `hac' nbclust(`nbclust') ///
	`bartlett' `partial' `correction' `network' `spatial' `storeweights' `storedistances' `small' `ak0'
}




***output
di in gr _col(55) "Number of obs = " in ye %8.0f e(N)
di in gr "Total (centered) SS     = " in ye %12.0g e(tss) _continue
di in gr _col(55) "Centered R2   = " in ye %8.4f e(r2)
di in gr "Total (uncentered) SS   = " in ye %12.0g e(tssu) _continue
di in gr _col(55) "Uncentered R2 = " in ye %8.4f e(r2u)
di in gr "Residual SS             = " in ye %12.0g e(rss) // _continue
*di in gr _col(55) "Root MSE      = " in ye %8.4g e(rmse)
di
ereturn display

restore
}

ereturn repost , esample(`touse')

********************************************************************************************

mata: st_numscalar("nofullrankXX", nofullrankXX)
mata: st_numscalar("nofullrankXZZZZX", nofullrankXZZZZX)

if nofullrankXX == 1 {
di in red "Warning: X'X matrix not of full rank. Some variables might be omitted."
di in red "Beta Coefficients and Standard Errors should be interpreted with caution."
}
if nofullrankXZZZZX == 1 {
di in red "Warning: X'PzX matrix not of full rank. Some variables might be omitted."
di in red "Beta Coefficients and Standard Errors should be interpreted with caution."
}

capt scalar drop nofullrankXX
capt scalar drop nofullrankXZZZZX
capt scalar drop factor_variables
capt scalar drop widstat
capt scalar drop hac_signal
capt scalar drop lagdistcutoff
capt scalar drop distcutoff
capt scalar drop lagcutoff
capt scalar drop ntsq
capt scalar drop ntsize
capt scalar drop correctr2_yn

capt scalar drop MSSu
capt scalar drop MSSc
capt scalar drop r2ur2
capt scalar drop r2r2
capt scalar drop TSScr2
capt scalar drop TSSur2
capt scalar drop RSSr2
capt scalar drop sumweigh

capt drop `depvar2'
capt drop `end2'
capt drop `exog2'
capt drop `ivv2'

capt drop fake_fe2
capt drop iddd_var

end
*
*******************************************************************************
*******************************************************************************
*				  **************** SUBROUTINES***************
*******************************************************************************
*******************************************************************************
*
capt program drop acreg_core
program acreg_core, eclass
	version 12
	syntax   [anything(name=regs)]  [if] [in]  [pweight fweight] [, depvar(string)  end(string) iv(string) ///
	dist_mat(varlist) links_mat(varlist) weights(varlist) ///
	latitude(varname) longitude(varname) id(varname) time(varname) ///
	LAGcutoff(integer 0) DISTcutoff(real 1) LAGDISTcutoff(integer 0) ///
	cluster(varlist) hac nbclust(integer -1) ///
	bartlett partial correction network spatial storeweights storedistances small ak0 ]

		tempname b V weightsmat AVCVb_fs

if ("`nbclust'" < "0"){
	local N_clust = 100
}
else {
	local N_clust = `nbclust'
}

*subrutines
*qui do acreg/acreg_no_correction.do
*qui do acreg/acreg_netw_dist_no_bartlett.do
*qui do acreg/acreg_netw_dist_bartlett.do
*qui do acreg/acreg_netw_links_no_bartlett.do
*qui do acreg/acreg_netw_links_bartlett.do
*qui do acreg/acreg_spat_dist_bartlett.do
*qui do acreg/acreg_spat_dist_no_bartlett.do
*qui do acreg/acreg_spat_ll_bartlett.do
*qui do acreg/acreg_spat_ll_no_bartlett.do
*qui do acreg/acreg_weig.do
*qui do acreg/acreg_cluster.do
*qui do acreg/nwacRanktest.ado
*qui do acreg/acreg_ranktest.do

*-------------------------------------------------------
*A) Run reg to check for missing values and etc

if "`weight'" == "" {
qui reg `depvar' `regs' `end' `iv'  `if' `in'
local wtexp = ""
 }
 else {
qui reg `depvar' `regs' `end' `iv'  `if' `in' [`weight'`exp']
local wtexp = "[`weight'`exp']"
}

qui keep if e(sample)==1



* check if links_mat or dist_mat is dimension N
check_n_dim `links_mat' , id(`id')
check_n_dim `dist_mat' , id(`id')

* check if weights_mat is dimension NT
check_nt_dim `weights' , id(`id')



tempvar wvar wf wff wvarfake

		if "`weight'" == "fweight" | "`weight'"=="aweight" {
			local wtexp `"[`weight' `exp']"'
			qui gen double `wvar' `exp'
			qui gen byte `wvarfake' = 1
		}
		if "`weight'" == "fweight" & "`kernel'" !="" {
			di in red "fweights not allowed (data are -tsset-)"
			exit 101
		}
		if "`weight'" == "fweight" & "`sw'" != "" {
			di in red "fweights currently not supported with -sw- option"
			exit 101
		}
		if "`weight'" == "iweight" {
			if "`robust'`cluster'`gmm2s'`kernel'" !="" {
				di in red "iweights not allowed with robust or gmm"
				exit 101
			}
			else {
				local wtexp `"[`weight' `exp']"'
				qui gen double `wvar' `exp'
			}
		}
		if "`weight'" == "pweight" {
			local wtexp `"[aweight `exp']"'
			qui gen double `wvar' `exp'
			qui gen double `wvarfake' `exp'
			local robust "robust"
		}
		if "`weight'" == "" {
* If no weights, define neutral weight variable
			qui gen byte `wvar'=1
			qui gen byte `wvarfake' = 1
		}

********************************************************************************
// weight factor and sample size
// Every time a weight is used, must multiply by scalar wf ("weight factor")
// wf=1 for no weights, fw and iw, wf = scalar that normalizes sum to be N if aw or pw

	qui sum `wvar' `wtexp', meanonly

// Weight statement
		if "`weight'" ~= "" {
di in gr "(sum of wgt is " %14.4e `r(sum_w)' ")"
scalar sumweigh = `r(sum_w)'
		}
		if "`weight'"=="" | "`weight'"=="fweight" | "`weight'"=="iweight" {
// Effective number of observations is sum of weight variable.
// If weight is "", weight var must be column of ones and N is number of rows
			scalar `wf'=1
			local NNN=r(sum_w)
		}
		else if "`weight'"=="aweight" | "`weight'"=="pweight" {
			scalar `wf'=r(N)/r(sum_w)
			local NNN=r(N)
		}
		else {
// Should never reach here
di as err " error - misspecified weights"
			exit 198
		}
		if `NNN'==0 {
di as err "no observations"
			exit 2000
		}

*-------------------------------------------------------
*B) Construct the matrices in Mata

//if `touse' == 1 {
	if ("`time'" == ""){
		qui sort `id'
	}
	else {
	qui sort `id' `time'
	}

	scalar define ntsize = _N
	scalar define ntsq = ntsize*ntsize
	scalar lagcutoff = `lagcutoff'
	scalar distcutoff = `distcutoff'
	scalar lagdistcutoff = `lagdistcutoff'
	if ("`hac'" == ""){
		scalar hac_signal = 0
	}
	else {
		scalar hac_signal = 1
	}


* Destring ID (if needed)
local vartype: type `id'
	if substr("`vartype'",1,3)=="str" {
		qui tempvar id_enc
		encode `id', generate(`id_enc')
	*	drop `id'
		qui tempvar id
		gen  `id'=`id_enc'
	}
*
	mata {
	lagcutoff_s = st_numscalar("lagcutoff")
	distcutoff_s = st_numscalar("distcutoff")
	lagdistcutoff_s = st_numscalar("lagdistcutoff")
	hac_signal_s = st_numscalar("hac_signal")

	*st_view(id_vec=.,.,"`id_enc'")
	id_vec = st_data(., "`id'")
	}
	if ("`time'" == ""){
	mata: time=J(rows(id_vec),1,1)
	}
	else {
	mata: st_view(time=.,.,"`time'")
	}

//	qui sort `id' `time'

	scalar `wff' = `wf'^(1/2)
	tempvar wvar2
	qui gen `wvar2' = `wf'*`wvar'
	qui	replace `wvar'  = sqrt(`wvar')
	qui	replace `wvar' = `wvar' * `wff'
	qui	replace `wvarfake'  = sqrt(`wvarfake')
	qui	replace `wvarfake' = `wvarfake' * `wff'

	mata {
	ntsq = st_numscalar("ntsq")
	st_view(Y=., ., "`depvar'")
	st_view(X=., ., "`end' `regs'")
	st_view(Z=., ., "`iv' `regs'")
	st_view(KX=., ., "`regs'") 			// -------------------------- KP
	st_view(KZ=., ., "`iv'")			// -------------------------- KP
	st_view(KY=., ., "`end'")        	// -------------------------- KP
	st_view(wvar=., ., "`wvar'")
	st_view(wvar2=., ., "`wvar2'")
	st_view(wvarfake=., ., "`wvarfake'")
	}

*	qui drop id_enc
mata {
X=(X, J(rows(Y),1,1))
Z=(Z, J(rows(Y),1,1))
K = cols(X)

Yn = Y
Y=	Y:*wvar
Zfake = Z:*wvarfake
X=	X:*wvar
Z=	Z:*wvar

//crossproducts
ZZ = quadcross(Z,Z)
ZZin = invsym(ZZ)
XZ = quadcross(X,Z)
ZX = quadcross(Z,X)
ZY = quadcross(Z,Y)

nofullrankXX=0
nofullrankXZZZZX=0


if (Z==X) {
// OLS: parameters
b=ZZin*ZY
P=ZZin*Z'
Pfake=ZZin*Zfake'

ranzz = rank(ZZ)
colzz = cols(ZZ) - factor_variables
if (ranzz == colzz) {
}
else {
nofullrankXX=1
}

}
else {
// 2SLS: parameters
XZZZZX = makesymmetric(XZ*ZZin*ZX)
XZZZZXin = invsym(XZZZZX)
XZZZZY = XZ*ZZin*ZY
b = XZZZZXin*XZZZZY
P = XZZZZXin*XZ*ZZin*Z'
Pfake = XZZZZXin*XZ*ZZin*Zfake'

ranxzx = rank(XZZZZX)
colxzx = cols(XZZZZX) - factor_variables
if (ranxzx == colxzx) {
}
else {
nofullrankXZZZZX=1
}
}


res = Y-X*b
N=rows(res)

Nend = cols(KY) 	// -------------------------- KP
Niv = cols(KZ) 		// -------------------------- KP
Nie = Nend*Niv		// -------------------------- KP

*------------------------------------------------------------------------
*        -------------------------- KP --------------------------	 	*
* Partial out X, p refers to partialling out
KX = (KX, J(rows(Y),1,1))
//KX = (KX, J(rows(Y),1,1))

KX=	KX:*wvar 	// -------------------------- KP
KZ=	KZ:*wvar 	// -------------------------- KP
KY= KY:*wvar 	// -------------------------- KP

Zp1 = KX'*KZ
Zp2 = KX*invsym(KX'KX)*Zp1
Zp = KZ - Zp2

Yp1 = KX'*KY
Yp2 = KX*invsym(KX'KX)*Yp1
Yp = KY - Yp2

yp1 = KX'*Y
yp2 = KX*invsym(KX'KX)*yp1
yp = Y - yp2

* VCV first stage parameters
mata: u1 = (invsym(Zp'Zp)*Zp') * Yp

u1 = (invsym(Zp'Zp)*Zp') * Yp
u2 = Zp * u1
u = Yp - u2

Sz = 1/sqrt(N)*Zp'
SzN = Sz

*------------------------------------------------------------------------
}




//}
*


//ak0 small sample correction
if ("`ak0'" == ""){
}
else {
mata: avgY = colsum(Y) :/ colnonmissing(Y)
mata: res = Y :- avgY
}


if ("`cluster'" == ""){
	if ("`weights'" == ""){
		if ("`network'" == ""){
			if ("`spatial'" == ""){
//1) NO CLUSTERING
				if ("`time'" == ""){
					mata: P = Pfake
				}
				else {
					mata: P = P
				}
				acreg_no_corr
			}
			else {
//2) SPATIAL
			if ("`time'" == ""){
			if (`distcutoff'>0) {
				mata: P = P
			}
			else {
				mata: P = Pfake
			}
			}
			else {
				mata: P = P
			}
				if ("`dist_mat'"==""){
					mata: st_view(lat=.,.,"`latitude'")
					mata: st_view(lon=.,.,"`longitude'")
					if ("`bartlett'"=="bartlett"){
						acreg_spat_ll_bart ,  `storeweights' `storedistances'
					}
					else {
						acreg_spat_ll_no_bart ,  `storeweights' `storedistances'
					}

				}
				else{
					mata: st_view(distance=., ., "`dist_mat'")
					if ("`bartlett'"=="bartlett"){
						acreg_spat_dist_bart , `storeweights' `storedistances'
					}
					else {
						acreg_spat_dist_no_bart , `storeweights' `storedistances'
					}
				}

			}
		}
		else {
//3) NETWORK
			if ("`time'" == ""){
			if (`distcutoff'>0) {
				mata: P = P
			}
			else {
				mata: P = Pfake
			}
			}
			else {
				mata: P = P
			}
			if ("`dist_mat'"==""){
				mata: st_view(links=., ., "`links_mat'")
				if ("`bartlett'"=="bartlett"){
					acreg_netw_links_bart , `storeweights' `storedistances'
				}
				else {
					acreg_netw_links_no_bart , `storeweights' `storedistances'
				}
			}
			else{
				mata: st_view(distance=., ., "`dist_mat'")
				if ("`bartlett'"=="bartlett"){
					acreg_netw_dist_bart ,	 `storeweights' `storedistances'
				}
				else {
					acreg_netw_dist_no_bart , `storeweights' `storedistances'
				}
			}
		}
	}
	else {
//4) WEIGHTS
		mata: weig = st_data(., "`weights'")
		acreg_weig
	}
}
else {
//5) MULTIWAY CLUSTERING
mata: st_view(cluster_mat=., ., "`cluster'")
acreg_cluster 	,	 `storeweights'
}


*if ("`correction'"=="correction"){
*		mata: mean_nonzero = nonzero / ntsq
*		mata: Nclus = 1 / mean_nonzero
*		mata: AVCVb_fs = (Nclus/(Nclus -1)) * AVCVb_fs
*	}



// export
	local temp_ct  : word count `regs'
	if `temp_ct' > 0 {
		prog_fvexpand `regs'
		local regss `r(newfvvars)'
	di in gr "Included instruments:" in ye " `regss'"
	local fvops `r(fvops)'
	}
	local temp_ct  : word count `end'
	if `temp_ct' > 0 {
		prog_fvexpand `end'
		local ends `r(newfvvars)'
	di in gr "Instrumented:"in ye " `ends'"
	if "`fvops'" == "true" {
		}
		else {
		local fvops `r(fvops)'
		}
	}
	local temp_ct  : word count `iv'
	if `temp_ct' > 0 {
		prog_fvexpand `iv'
		local ivs `r(newfvvars)'
	di in gr "Excluded instruments:" in ye " `ivs'"
	if "`fvops'" == "true" {
		}
		else {
		local fvops `r(fvops)'
		}
	}

	//small sample correction
if ("`small'" == ""){
local n_cluster=1000000
}
else {
	local iv_ct : word count `ivs'
	local sdofminus = 0
	local partial_ct : word count `regss'
	local partial_ct = `partial_ct' + 1
	local sdofminus =`sdofminus'+`partial_ct'
	mata: mean_nonzero = nonzero / ntsq
	mata: n_cluster = 1 / mean_nonzero
	mata: st_local("N", strofreal(N))
	mata: st_local("n_cluster", strofreal(n_cluster))
	local n_cluster `n_cluster'
	mata: rankxx = rows(W) - diag0cnt(W)
	mata: st_local("rankxx", strofreal(rankxx))
	local ss_corr = (`N'-1) / (`N'-`rankxx'-`sdofminus') *`n_cluster'/(`n_cluster'-1)
	mata ss_corr = st_numscalar("ss_corr")
	mata: V = V * `ss_corr'
	mata: AVCVb_fs = (n_cluster/(n_cluster -1)) * AVCVb_fs
}


mata {
	b=b'
	st_matrix("r(V)", V)
	st_matrix("r(b)", b)
	st_numscalar("r(N)", N)

	RSS=(res'res)
	//TSSc = sum((Y:-mean(Y)):^2)
	TSSc = sum(wvar2'*((Yn:-mean(Yn, wvar2)):^2))
	TSSu = sum((Y):^2)

	Rsqc = 1 - RSS/TSSc
	Rsqu = 1 - RSS/TSSu

	Fdf1 = K-1
	Fdf2 = N - K

	*Fsta = r(chi2)/`Fdf1' * `df_r'/(`N'-`dofminus')

	MSSc = TSSc - RSS
	MSSu = TSSu - RSS

	st_numscalar("r(RSS)", RSS)

	st_numscalar("r(TSSc)", TSSc)
	st_numscalar("r(TSSu)", TSSu)

	st_numscalar("r(MSSc)", MSSc)
	st_numscalar("r(MSSu)", MSSu)

	st_numscalar("r(Rsqc)", Rsqc)
	st_numscalar("r(Rsqu)", Rsqu)

	st_numscalar("r(Fdf1)", Fdf1)
	st_numscalar("r(Fdf2)", Fdf2)

*	st_matrix("r(AVCVb_fs)", AVCVb_fs)

}


mat `b'=r(b)
mat `V'=r(V)

matname `V' `ends' `regss' "_cons" , e
mat colnames `b' = `ends' `regss' _cons

ereturn post `b' `V'
ereturn local depvar "`depvar'"

ereturn scalar N=r(N)
if ("`weight'" == "fweight") {
ereturn scalar N = sumweigh
}
ereturn local cmd "acreg"

ereturn scalar rss=r(RSS)

ereturn scalar tss=r(TSSc)
ereturn scalar tssu=r(TSSu)

ereturn scalar mss=r(MSSc)
ereturn scalar mssu=r(MSSu)

ereturn scalar r2=r(Rsqc)
ereturn scalar r2u=r(Rsqu)

ereturn scalar Fdf1=r(Fdf1)
ereturn scalar Fdf2=r(Fdf2)

*mat `AVCVb_fs'=r(AVCVb_fs)

*
*-------------------------------------------------------
*Z) KP TEST
********************************************************************


/*
if ("`iv'" == ""){
}
else {

acreg_ranktest , end(`end') iv(`iv')


*	mata: st_matrix("r(AVCVb_fs)", AVCVb_fs)
*	mat AVCVb_fs=r(AVCVb_fs)

	*** We retrieve code for IVreg2

	* Stata convention is to exclude constant from instrument list
	* Need word option so that varnames with "_cons" in them aren't zapped
	qui count
	local N = r(N)
	local iv_ct : word count `ivs'
	local iv_ct  = `iv_ct' + 1
	local sdofminus = 0
	local partial_ct : word count `regs'
	local partial_ct = `partial_ct' + 1
	local sdofminus =`sdofminus'+`partial_ct'
	local N_clust=1000	// November 2017, beta version: this will become endogenous in the future
*	local N_clust "`n_cluster'"
	local exex1_ct     : word count `iv'
	*local noconstant "noconstant"
*	local noconstant ""
*	local robust "robust"
	*local robust ""
*	local CLUS="AVCVb_fs"

	* (line 1735 in ivreg2)
	*Need only test of full rank
*	qui nwacRanktest (`end') (`iv') , vb("`CLUS'") partial(`regs') full wald `noconstant' `robust'
					/* `robust' `clopt' `bwopt' `kernopt'*/

	* sdofminus used here so that F-stat matches test stat from regression with no partial
	scalar rkf=chi2/(`N'-1) *(`N'-`iv_ct'-`sdofminus')  *(`N_clust'-1)/`N_clust' /`exex1_ct'

	scalar widstat=rkf
	scalar KPstat=widstat

	di "Kleibergen Paap rk statistic =  " KPstat
	********************************************************************
}
*/

*******************************************************************************************
* Rank identification and redundancy block
*******************************************************************************************
*Note:  included instruments (or included exogenous) = regs
*		endogenous = end
*		excluded exogenous = iv

	local iv_ct 		: word count `iv' `regs'
	local iv_ct  = `iv_ct' + 1
	local sdofminus = 0
	local partial_ct : word count `regs'
	local partial_ct = `partial_ct' + 1
	local sdofminus =`sdofminus'+`partial_ct'
	local dofminus = 0
	local exex_ct     : word count `iv'
	local endo1_ct     : word count `end'
	local N = r(N)

	*local N_clust = 100

	mata: st_matrix("r(AVCVb_fs)", AVCVb_fs)
	mat AVCVb_fs=r(AVCVb_fs)
	local AVCVb_fs="AVCVb_fs"


if ("`cluster'" == ""){
	if ("`weights'" == ""){
		if ("`network'" == ""){
			if ("`spatial'" == ""){
//1) NO CLUSTERING
			local ncc = 1
			}
		}
	}
//5) CLUSTERING
local ncc = 0
}
else {
//5) MULTIWAY CLUSTERING
local ncc = 0
}



if `endo1_ct' > 0  {
// id=underidentification statistic, wid=weak identification statistic
			tempname idrkstat widrkstat iddf idp
			tempname ccf cdf rkf cceval cdeval cd
			tempname idstat

// WEAK IDENTIFICATION
// Weak id statistic is Cragg-Donald F stat, rk Wald F stat if not
// ranktest exits with error if not full rank so can use iv1_ct and rhs1_ct etc.
	//	if "`robust'`cluster'`kernel'"=="" {
	//			scalar `widstat'=`cdf'
	//		}
	//		else {
// Need only test of full rank
				cap  acreg_ranktest						///
								(`end')			///
								(`iv')			///
								,				///
								partial(`regs')	///
								full					///
								wald					///
								mat_clst("`AVCVb_fs'")	///
								robust

							*	noconstant				///
							*	`wtexp'					///
							*	`clopt'					///
							*	`bwopt'					///
							*	`kernopt'
// Canonical correlations returned in r(ccorr), sorted in descending order.
// If largest = 1, collinearities so enter error block.
				local rkerror		= _rc>0 | r(chi2)==.
				if ~`rkerror' {
					local rkerror	= el(r(ccorr),1,1)==1
				}
				if `rkerror' {
di in ye "warning: -ranktest- error in calculating KP weak identification test statistics;"
di in ye  "may be caused by collinearities or the presence of factor variables or weights"
					scalar `rkf'		= .
					scalar widstat	= .
				}
				else {
// sdofminus used here so that F-stat matches test stat from regression with no partial
					if `ncc' == 1 {
						scalar `rkf'=r(chi2)/r(N)*(`N'-`iv_ct'-`sdofminus'-`dofminus')/`exex_ct'
					}
					else {
						scalar `rkf' =	r(chi2)/(`N'-1) *				///
										(`N'-`iv_ct'-`sdofminus') *	///
										(`N_clust'-1)/`N_clust' /		///
										`exex_ct'
					}
					scalar widstat=`rkf'
				}

	//		}	// end weak-identification stat

		ereturn scalar widstat=widstat
	*	di "chi2 = " r(chi2)
	*	di "N = " r(N)
	*	di "N = " `N'
	*	di "iv ct = " `iv_ct'
	*	di "sdof = "`sdofminus'
	*	di "dof = "`dofminus'
	*	di "ex = "`exex_ct'
	*	di "rkf = " `rkf'
*di in gr "Kleibergen-Paap rk Wald F statistic:" in ye _col(71) %8.3f e(widstat)
di in gr "Kleibergen-Paap rk Wald F statistic:" in ye  %8.3f e(widstat)

		ereturn scalar widstat=e(widstat)
}	// end under- and weak-identification stats
end


**************************************************************************************
capt program drop PEnd
program define PEnd, sclass
	version 8.2
	if `"`0'"' == "[" {
		sret local stop 1
		exit
	}
	if `"`0'"' == "," {
		sret local stop 1
		exit
	}
	if `"`0'"' == "if" {
		sret local stop 1
		exit
	}
	if substr(`"`0'"',1,3) == "if(" {
		sret local stop 1
		exit
	}
	if `"`0'"' == "in" {
		sret local stop 1
		exit
	}
	if `"`0'"' == "" {
		sret local stop 1
		exit
	}
	else	sret local stop 0
end


*******************************************************************************************
capt program drop prog_fvexpand
program prog_fvexpand, rclass
		syntax varlist(numeric fv)
		fvexpand `varlist'
		local fvops = r(fvops)
		if "`fvops'" == "true" {
		local newfvvarlist = "`r(varlist)'"
		return local newfvvars `"`newfvvarlist'"'
		return local fvops `"`fvops'"'
		}
		else {
		return local newfvvars `"`varlist'"'
		}
	end


*******************************************************************************************
capt program drop prog_fv_unab
program prog_fv_unab, rclass
	syntax varlist(numeric fv)
	local temp_vt  : word count `varlist'
	forval j=1/`temp_vt' {
		local this_var `: word `j' of `varlist''
		if substr("`this_var'", 1, 2) == "i." {
			fvrevar `this_var', list
			local this_var = "`r(varlist)'"
			return local this_var `"`this_var'"'
			tsunab this_var : `this_var'
			local this_var = "i.`this_var'"
		}
		else {
		tsunab this_var : `this_var'
		}
		local ts_varlist `ts_varlist' `this_var'
	}
	return local ts_varlist `"`ts_varlist'"'
end


*********************************************************************************************
capt program drop acreg_r2partial
program acreg_r2partial, eclass
	version 12
	syntax   [anything(name=regs)]  [if] [in]  [pweight] [, depvar(string)  end(string) iv(string) ]

if "`weight'" == "" {
qui ivreg2 `depvar' `regs' (`end'=`iv')  `if' `in'
 }
 else {
qui ivreg2  `depvar' `regs' (`end'=`iv')  `if' `in' [`weight'`exp']
}

scalar RSSr2=e(rss)
scalar TSSur2= e(yy)
scalar TSScr2= e(yyc)

scalar r2r2=e(r2)
scalar r2ur2=e(r2u)


mata {
	RSSr2 = st_numscalar("RSSr2")
	TSScr2 = st_numscalar("TSScr2")
	TSSur2 = st_numscalar("TSSur2")
	r2r2 = st_numscalar("r2r2")
	r2ur2 = st_numscalar("r2ur2")
	}

end


*********************************************************************************************
capt program drop verify
program define verify
    syntax anything(everything equalok) [, MSG(string asis) RC(integer 198)]
    if !(`anything') {
        di as error `msg'
        exit 198
    }
end

*********************************************************************************************
capt program drop check_missing
program define check_missing
syntax [anything(name=anything)]
	if "`anything'" == "" {
	}
	else {
		foreach cv of varlist `anything' {
			qui count if missing(`cv')
			local missing_val = r(N)
			if (`missing_val' != 0) {
				di as error "`cv' contains missing values"
				exit 198
			}
		}
	}
end


*********************************************************************************************
capt program drop check_number
program define check_number
syntax [anything(name=anything)]
	if "`anything'" == "" {
	}
	else {
		foreach cv of varlist `anything' {
			capture confirm numeric variable `cv'
			if !_rc {
			}
			else {
				di as error "`cv' is not a numeric variable"
				exit 198
			}
		}
	}
end


*********************************************************************************************
capt program drop check_lat_lon
program define check_lat_lon
syntax [anything(name=anything)]
	if "`anything'" == "" {
	}
	else {
		foreach cv of varlist `anything' {
			qui sum `cv'
			local max_cv = r(max)
			local min_cv = r(min)
			if (`max_cv' > 180 | `min_cv' < -180) {
				di as error "`cv' must be between -180 and 180"
				exit 198
			}
		}
	}
end

*********************************************************************************************
capt program drop check_links
program define check_links
syntax [anything(name=anything)]
	if "`anything'" == "" {
	}
	else {
		foreach cv of varlist `anything' {
			qui sum `cv'
			local max_cv = r(max)
			local min_cv = r(min)
			if (`max_cv' > 1 | `min_cv' < 0) {
				di as error "`cv' must be a binary variable"
				exit 198
			}
		}
	}
end

*********************************************************************************************
capt program drop check_link_dim
program define check_link_dim
syntax [anything(name=anything)] , id(varname)
	if "`anything'" == "" {
	}
	else {
		local nu_links = 0
		foreach cv of varlist `anything' {
			local nu_links = `nu_links' + 1
			}
		qui tab `id'
		local nu_ind = r(r)
		if (`nu_ind' !=  `nu_links') {
			di as error "number of link variables should coincide with number of individuals"
			di in ye "this can be due to if/in function"
			exit 198
		}
	}
end

*********************************************************************************************
capt program drop check_n_dim
program define check_n_dim
syntax [anything(name=anything)] , id(varname)
	if "`anything'" == "" {
	}
	else {
		local nu_vars = 0
		foreach cv of varlist `anything' {
			local nu_vars = `nu_vars' + 1
			}
		qui tab `id'
		local nu_ind = r(r)
		if (`nu_ind' !=  `nu_vars') {
			di as error "number of link/distance variables should coincide with number of individuals"
			di in ye "this can be due to if/in function"
			capt drop iddd_var
			exit 198
		}
	}
end

*********************************************************************************************
capt program drop check_nt_dim
program define check_nt_dim
syntax [anything(name=anything)] , id(varname)
	if "`anything'" == "" {
	}
	else {
		local nu_vars = 0
		foreach cv of varlist `anything' {
			local nu_vars = `nu_vars' + 1
			}
		qui sum `id'
		local nu_ind = r(N)
		if (`nu_ind' !=  `nu_vars') {
			di as error "number of weight variables should coincide with number of observations"
			di in ye "this can be due to if/in function"
			capt drop iddd_var
			exit 198
		}
	}
end


*******************************************************************************************
capt program drop prog_fv_partial
program prog_fv_partial, rclass
	syntax varlist(numeric fv)
	local temp_vt  : word count `varlist'
	forval j=1/`temp_vt' {
		local this_var `: word `j' of `varlist''
		if substr("`this_var'", 1, 2) == "i." {
		di in red "factor variables not allowed with pfe1 option, please specify" in ye " xi: " in red "before acreg"
		exit 198
		}
	}
end


*******************************************************************************************
capt program drop prog_count_factor_variables
program prog_count_factor_variables, rclass
	syntax varlist(numeric fv)
	scalar factor_variables=0
	local temp_vt  : word count `varlist'
	forval j=1/`temp_vt' {
		local this_var `: word `j' of `varlist''
		if substr("`this_var'", 1, 2) == "i." {
		scalar factor_variables=factor_variables+1
		}
	}
end


*********************************************************************************************
*********************************************************************************************
*********************************************************************************************

*********************************************************************************************

capt program drop acreg_weig
program acreg_weig, eclass
	version 12
	syntax [, ]

mata {
weight_cols = cols(weig)
weight_rows = rows(weig)

st_local("N_o", strofreal(N))
st_local("weight_cols", strofreal(weight_cols))
st_local("weight_rows", strofreal(weight_rows))

}

if `N_o' != `weight_cols' {
	disp as error "The number of columns in the wights matrix is different than the number of observations"
exit 498
}


mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0

AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP // -------------------------- KP

		V =  P * (res * res' :* weig) * P'
		weigz = weig :> 0
		nonzero = sum(weigz)

		AVCVb_fs = (I(Nend)#SzN) * ((vec(u) * (vec(u))') :* (J(Nend,Nend, 1)#weig)) * (I(Nend)#SzN)'		// -------------------------- KP

//*mata: V[1..5,1..5]

_makesymmetric(V)

}

end


*********************************************************************************************

capt program drop acreg_spat_ll_no_bart
program acreg_spat_ll_no_bart, eclass
	version 12
	syntax [,  storeweights storedistances ] // storedistances to be implemented


*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

// values
a = .01745329
c = 6371
pi = 3.141593

Idunique = uniqrows(id_vec)
NId = rows(Idunique)

*some vars for checks
**dist= J(N,NId, 0)
**adjj= J(N,N, 0)



if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			lon_scale = cos(lat[s,1]*pi()/180)*111
			lat_scale = 111
			distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
					  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

**dist[.,i] = distance_i // check

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		else {
			weight_dis = (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP

**adjj[.,ind+1] = w[.,1] // check
	}

	ind = ind + t1
}

*
st_matrix("weightsmat", weigg)
}
else {

if ("`storedistances'" == "storedistances"){
distances = J(N,NId, 0)

for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			lon_scale = cos(lat[s,1]*pi()/180)*111
			lat_scale = 111
			distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
					  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

	distances[.,i] = distance_i[.,1]

	for (t = 1; t <= t1; t++){
	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		else {
			weight_dis = (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP

**adjj[.,ind+1] = w[.,1] // check
	}

	ind = ind + t1
}

st_matrix("distancesmat", distances)
}
else {
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			lon_scale = cos(lat[s,1]*pi()/180)*111
			lat_scale = 111
			distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
					  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}


	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		else {
			weight_dis = (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP

**adjj[.,ind+1] = w[.,1] // check
	}

	ind = ind + t1
}

*
}
}
*

_makesymmetric(V)

}
end


*********************************************************************************************

capt program drop acreg_spat_ll_bart
program acreg_spat_ll_bart, eclass
	version 12
	syntax  [, storeweights storedistances ] // storedistances to be implemented

*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

// values
a = .01745329
c = 6371
pi = 3.141593

Idunique = uniqrows(id_vec)
NId = rows(Idunique)

*some vars for checks
**dist= J(N,NId, 0)
**adjj= J(N,N, 0)

if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			lon_scale = cos(lat[s,1]*pi()/180)*111
			lat_scale = 111
			distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
					  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

*
st_matrix("weightsmat", weigg)
}
else {

if ("`storedistances'" == "storedistances"){
distances = J(N,NId, 0)

for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			lon_scale = cos(lat[s,1]*pi()/180)*111
			lat_scale = 111
			distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
					  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

	distances[.,i] = distance_i[.,1]

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}



st_matrix("distancesmat", distances)
}
else{
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			lon_scale = cos(lat[s,1]*pi()/180)*111
			lat_scale = 111
			distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
					  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

*
}
}
*

_makesymmetric(V)

}
end

*********************************************************************************************

capt program drop acreg_spat_dist_no_bart
program acreg_spat_dist_no_bart, eclass
	version 12
	syntax [, storeweights  ]


*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

// values
a = .01745329
c = 6371
pi = 3.141593

Idunique = uniqrows(id_vec)
NId = rows(Idunique)
d_cols = cols(distance)

st_local("NId", strofreal(NId))
st_local("d_cols", strofreal(d_cols))

}

if `NId' != `d_cols' {
	disp as error "The number of columns in the distance matrix is different than the number of individuals"
exit 498
}

mata {

if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			distance_i = distance[.,i]
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		else {
			weight_dis = (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}
*
st_matrix("weightsmat", weigg)
}
else {
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			distance_i = distance[.,i]
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		else {
			weight_dis = (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}


*
}
*

_makesymmetric(V)

}

end

*********************************************************************************************

capt program drop acreg_spat_dist_bart
program acreg_spat_dist_bart, eclass
	version 12
	syntax  [,storeweights ]

*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

// values
a = .01745329
c = 6371
pi = 3.141593

Idunique = uniqrows(id_vec)
NId = rows(Idunique)
d_cols = cols(distance)

st_local("NId", strofreal(NId))
st_local("d_cols", strofreal(d_cols))

}

if `NId' != `d_cols' {
	disp as error "The number of columns in the distance matrix is different than the number of individuals"
exit 498
}

mata {

if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			distance_i = distance[.,i]
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

*
st_matrix("weightsmat", weigg)
}
else {
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector

			distance_i = distance[.,i]

			*lon_scale = cos(lat[s,1]*pi()/180)*111
			*lat_scale = 111
			*distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
			*		  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s))
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

*
}
*
//*mata: V[1..5,1..5]

_makesymmetric(V)

}

end


*********************************************************************************************

capt program drop acreg_no_corr
program acreg_no_corr, eclass
	version 12
	syntax [,  storeweights  ]


*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/


mata {
// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

Idunique = uniqrows(id_vec)
NId = rows(Idunique)

	if ("`storeweights'" == "storeweights"){
		weigg = J(N,N, 0)
		for (i = 1; i <= NId; i++){
			rows_ni = id_vec:==Idunique[i,1]
			// __ Get subsets of variables for ID i (without changing original matrix)
			time1 = select(time, rows_ni)
			t1 = length(time1)
			// __ Weights ID i vector
			w = J(N,1, 0)
			// __ Indicator variables (start, end)
			s = ind + 1
			e = ind + t1

			for (t = 1; t <= t1; t++){
				// __ Indicator variables (position in the vector)
				pos  = ind + t
				if (lagcutoff_s>0) {
				// __ Vector of weights (same ID)
					time_gap_i = abs(time1[t,1] :- time1)
					if (hac_signal_s==1) {
						weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
					}
					else {
						weight = (time_gap_i :<= lagcutoff_s)
					}
					w[s..e,1] = weight
				}
				else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
				}
				// __ VCV
				wz = w :> 0
				nonzero = nonzero + sum(wz)
				V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
				V = V + V1
				weigg[.,pos] = w[.,1]
				AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP // -------------------------- KP
				AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
			}
		ind = ind + t1
		}
	*
		st_matrix("weightsmat", weigg)
	}
	else {
		for (i = 1; i <= NId; i++){
			rows_ni = id_vec:==Idunique[i,1]
			// __ Get subsets of variables for ID i (without changing original matrix)
			time1 = select(time, rows_ni)
			t1 = length(time1)
			// __ Weights ID i vector
			w = J(N,1, 0)
			// __ Indicator variables (start, end)
			s = ind + 1
			e = ind + t1

			for (t = 1; t <= t1; t++){
				// __ Indicator variables (position in the vector)
				pos  = ind + t

				if (lagcutoff_s>0) {
					// __ Vector of weights (same ID)
					time_gap_i = abs(time1[t,1] :- time1)
					if (hac_signal_s==1) {
						weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
					}
					else {
						weight = (time_gap_i :<= lagcutoff_s)
					}
					w[s..e,1] = weight
				}
				else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
				}
				// __ VCV
				wz = w :> 0
				nonzero = nonzero + sum(wz)
				V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
				V = V + V1
				AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
				AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
			}
			ind = ind + t1
		}

	}
_makesymmetric(V)
}
end


*********************************************************************************************

capt program drop acreg_netw_links_no_bart
program acreg_netw_links_no_bart, eclass
	version 12
	syntax [, storeweights storedistances ]

*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

// values
a = .01745329
c = 6371
pi = 3.141593

Idunique = uniqrows(id_vec)
NId = rows(Idunique)
l_cols = cols(links)

st_local("NId", strofreal(NId))
st_local("l_cols", strofreal(l_cols))

}

if `NId' != `l_cols' {
	disp as error "The number of columns in the links matrix is different than the number of individuals"
exit 498
}

mata {


				//************************************************************
				// __ Creating Distances from links
				if (distcutoff_s>0) {
				// __ a) From NTxN to NxN

					distance = J(N,NId, 0)

					if (distcutoff_s==1) {
						distance = links
					}
					else {

						B_ind = 0
						NN_dist = J(NId,NId, 0)
						NN_links = J(NId,NId, 0)

						for (i = 1; i <= NId; i++){
							B_rows_ni = id_vec:==Idunique[i,1]

							B_time1 = select(time, B_rows_ni)
							B_t1 = length(B_time1)

							first_ind_line = B_ind + 1
							NN_links[i,.] = links[first_ind_line,.]

							B_ind = B_ind + B_t1
						}

				// __ b) Taking power of a matrix


						NN_links_pow = NN_links
						NN_dist = NN_links
						for (k = 1; k <= 100; k++){
				//powers
							if (distcutoff_s>k) {
								NN_links_pow = NN_links_pow * NN_links
								NN_links_pow_ind = (NN_links_pow:!= 0) :* (k+1) :* (NN_dist :== 0)
								NN_dist = NN_dist + NN_links_pow_ind
							}
						}

					// __ c) Form NxN to NTxN
						B_ind = 0
						for (i = 1; i <= NId; i++){
							B_rows_ni = id_vec:==Idunique[i,1]
							B_time1 = select(time, B_rows_ni)
							B_t1 = length(B_time1)

							for (t = 1; t <= B_t1; t++){
								ind_line = B_ind + t
								distance[ind_line,.] = NN_dist[i,.]
							}
							B_ind = B_ind + B_t1
						}
					}
				}
				//************************************************************

if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)

for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector

			distance_i = distance[.,i]

			*lon_scale = cos(lat[s,1]*pi()/180)*111
			*lat_scale = 111
			*distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
			*		  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		else {
			weight_dis =  (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}
	ind = ind + t1
}

*
st_matrix("weightsmat", weigg)
}
else {
if ("`storedistances'" == "storedistances"){
distances = J(N,NId, 0)

for (i = 1; i <= NId; i++){
	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector

			distance_i = distance[.,i]

			*lon_scale = cos(lat[s,1]*pi()/180)*111
			*lat_scale = 111
			*distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
			*		  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

	distances[.,i] = distance_i[.,1]

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		else {
			weight_dis =  (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}
	ind = ind + t1
}
*

st_matrix("distancesmat", distances)
}
else{


for (i = 1; i <= NId; i++){
	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector

			distance_i = distance[.,i]

			*lon_scale = cos(lat[s,1]*pi()/180)*111
			*lat_scale = 111
			*distance_i = ((lat_scale*(lat[s,1]:-lat[.,1])):^2 + ///
			*		  (lon_scale*(lon[s,1]:-lon[.,1])):^2):^0.5
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		else {
			weight_dis =  (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}
	ind = ind + t1
}
*
}
}
*
//*mata: V[1..5,1..5]

_makesymmetric(V)

}

end


*********************************************************************************************

capt program drop acreg_netw_links_bart
program acreg_netw_links_bart, eclass
	version 12
	syntax  [, storeweights storedistances]
*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

// values
a = .01745329
c = 6371
pi = 3.141593

Idunique = uniqrows(id_vec)
NId = rows(Idunique)
l_cols = cols(links)

st_local("NId", strofreal(NId))
st_local("l_cols", strofreal(l_cols))

}

if `NId' != `l_cols' {
	disp as error "The number of columns in the links matrix is different than the number of individuals"
exit 498
}

mata {
				//************************************************************
				// __ Creating Distances from links
				if (distcutoff_s>0) {
				// __ a) Form NTxN to NxN

					distance = J(N,NId, 0)

					if (distcutoff_s==1) {
						distance = links
					}
					else {

						B_ind = 0
						NN_dist = J(NId,NId, 0)
						NN_links = J(NId,NId, 0)

						for (i = 1; i <= NId; i++){
							B_rows_ni = id_vec:==Idunique[i,1]

							B_time1 = select(time, B_rows_ni)
							B_t1 = length(B_time1)

							first_ind_line = B_ind + 1
							NN_links[i,.] = links[first_ind_line,.]

							B_ind = B_ind + B_t1
						}

				// __ b) Taking power of a matrix


						NN_links_pow = NN_links
						NN_dist = NN_links
						for (k = 1; k <= 100; k++){
				//powers
							if (distcutoff_s>k) {
								NN_links_pow = NN_links_pow * NN_links
								NN_links_pow_ind = (NN_links_pow:!= 0) :* (k+1) :* (NN_dist :== 0)
								NN_dist = NN_dist + NN_links_pow_ind
							}
						}

					// __ c) Form NxN to NTxN
						B_ind = 0
						for (i = 1; i <= NId; i++){
							B_rows_ni = id_vec:==Idunique[i,1]
							B_time1 = select(time, B_rows_ni)
							B_t1 = length(B_time1)

							for (t = 1; t <= B_t1; t++){
								ind_line = B_ind + t
								distance[ind_line,.] = NN_dist[i,.]
							}
							B_ind = B_ind + B_t1
						}
					}
				}
				//************************************************************


if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			distance_i = distance[.,i]
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

*
st_matrix("weightsmat", weigg)
}
else {
if ("`storedistances'" == "storedistances"){
distances = J(N,NId, 0)

for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			distance_i = distance[.,i]
		}

	distances[.,i] = distance_i[.,1]

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

st_matrix("distancesmat", distances)
}
else{
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			distance_i = distance[.,i]
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

*
}
}
*
//*mata: V[1..5,1..5]

_makesymmetric(V)

}

end


*********************************************************************************************

capt program drop acreg_netw_dist_no_bart
program acreg_netw_dist_no_bart, eclass
	version 12
	syntax  [, storeweights ]

*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

// values
a = .01745329
c = 6371
pi = 3.141593

Idunique = uniqrows(id_vec)
NId = rows(Idunique)
d_cols = cols(distance)

st_local("NId", strofreal(NId))
st_local("d_cols", strofreal(d_cols))

}

if `NId' != `d_cols' {
	disp as error "The number of columns in the distance matrix is different than the number of individuals"
exit 498
}

mata {

if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector

			distance_i = distance[.,i]

		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		else {
			weight_dis = (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}
	ind = ind + t1
}

*
st_matrix("weightsmat", weigg)
}
else {
for (i = 1; i <= NId; i++){
	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector

			distance_i = distance[.,i]
		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		else {
			weight_dis =  (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight = (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}
	ind = ind + t1
}
*
}
*

_makesymmetric(V)

}

end


*********************************************************************************************

capt program drop acreg_netw_dist_bart
program acreg_netw_dist_bart, eclass
	version 12
	syntax [,  storeweights ]


*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP

// values
a = .01745329
c = 6371
pi = 3.141593

Idunique = uniqrows(id_vec)
NId = rows(Idunique)
d_cols = cols(distance)

st_local("NId", strofreal(NId))
st_local("d_cols", strofreal(d_cols))

}

if `NId' != `d_cols' {
	disp as error "The number of columns in the distance matrix is different than the number of individuals"
exit 498
}

mata {

if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			distance_i = distance[.,i]

		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

*
st_matrix("weightsmat", weigg)
}
else {
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)
	distance_i = J(N,1, 0)

	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)
		if (distcutoff_s>0) {

		// a) Computing Distance vector
			distance_i = distance[.,i]

		}

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
	if (distcutoff_s>0) {
	// a) Computing Weights
		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
		if (hac_signal_s==1) {
			weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)):*(1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		else {
			weight_dis = (1:-abs(distance_i :/ distcutoff_s)) :* ( (time_dis_gap_i :<= lagdistcutoff_s) :* (distance_i :<= distcutoff_s)) :* (distance_i :!= 0)
		}
		w[.,1] = weight_dis
	}
	if (lagcutoff_s>0) {
	// __ Vector of weights (same ID)
		time_gap_i = abs(time1[t,1] :- time1)
		if (hac_signal_s==1) {
			weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
		}
		else {
			weight =  (time_gap_i :<= lagcutoff_s)
		}
		w[s..e,1] = weight
	}
	else {
	w[s..e,1] = time1:==time[pos,1]
	//w[s..e,1] = J(t1,1,0)
	//w[pos,1] = 1
	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP
	}

	ind = ind + t1
}

*
}
*
//*mata: V[1..5,1..5]

_makesymmetric(V)

}

end


*********************************************************************************************

capt program drop acreg_cluster
program acreg_cluster, eclass
	version 12
	syntax [,  storeweights ] // storedistances to be implemented


*
/*
* First stage
mata: PZ=Z*invsym(Z'Z)*Z'
mata: MZ=I(rows(PZ))-PZ


mata: X1=Xe[.,1]
mata: X2=Xe[.,2]

mata: gamma1=PZ * X1
*mata: rows(gamma1)
mata: resX1 = MZ * X1

mata: VCV1 = PZ * (resX1 * resX1' :* cluster) * PZ'
mata: VCV1 = VCV1[|1,1\16,16|]

mata: F1 = (gamma1[.,1..Niv]'*invsym(VCV1)*gamma1[.,1..Niv])/Niv
*/

mata {

// 2SLS: VCV
V = J(K,K, 0)
nonzero = 0
ind = 0
pos = 0

AVCVb_fs_1 = J(Nie,Nie, 0) 		// -------------------------- KP
AVCVb_fs = J(Nie,Nie, 0)   		// -------------------------- KP


Idunique = uniqrows(id_vec)
NId = rows(Idunique)

if ("`storeweights'" == "storeweights"){
weigg = J(N,N, 0)
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)


	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
		// a) Computing link vector
		link_vec_temp = J(N,1, 0)
		link_vec_i = J(N,1, 0)
		cols_clus_mat = cols(cluster_mat)
		q=1
		for (q = 1; q <= cols_clus_mat; q++){
			link_vec_temp = link_vec_temp :+ (cluster_mat[pos,q] :== cluster_mat[.,q])
			link_vec_i = link_vec_temp:>0
		}


//	if (distcutoff_s>0) {
	// a) Computing Weights
//		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
//		weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* link_vec_i
//		w[.,1] = weight_dis
//	}
	w[.,1] = link_vec_i
//	if (lagcutoff_s>0) {
//	// __ Vector of weights (same ID)
//		time_gap_i = abs(time1[t,1] :- time1)
//		time_gap_cut_i = time_gap_i :<= lagcutoff_s
//		weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
//		w[s..e,1] = weight
//	}
//	else {
//	w[s..e,1] = J(t1,1,0)

	w[s..e,1] = time1:==time[pos,1]
	*w[pos,1] = 1

//	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		weigg[.,pos] = w[.,1]
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP

**adjj[.,ind+1] = w[.,1] // check
	}

	ind = ind + t1
}
st_matrix("weightsmat", weigg)
}
else {
for (i = 1; i <= NId; i++){

	rows_ni = id_vec:==Idunique[i,1]

	// __ Get subsets of variables for ID i (without changing original matrix)
	time1 = select(time, rows_ni)
	t1 = length(time1)

	// __ Weights ID i vector
	w = J(N,1, 0)


	// __ Indicator variables (start, end)
	s = ind + 1
	e = ind + t1

	// __ Vector of weights (different ID)

	for (t = 1; t <= t1; t++){

	// __ Indicator variables (position in the vector)
	pos  = ind + t
		// a) Computing link vector
		link_vec_temp = J(N,1, 0)
		link_vec_i = J(N,1, 0)
		cols_clus_mat = cols(cluster_mat)
		q=1
		for (q = 1; q <= cols_clus_mat; q++){
			link_vec_temp = link_vec_temp :+ (cluster_mat[pos,q] :== cluster_mat[.,q])
			link_vec_i = link_vec_temp:>0
		}


//	if (distcutoff_s>0) {
	// a) Computing Weights
//		time_dis_gap_i = abs(time1[t,1] :- time[.,1])
//		weight_dis = (1 :- time_dis_gap_i :/ (lagdistcutoff_s +1)) :* (time_dis_gap_i :<= lagdistcutoff_s) :* link_vec_i
//		w[.,1] = weight_dis
//	}
	w[.,1] = link_vec_i
//	if (lagcutoff_s>0) {
//	// __ Vector of weights (same ID)
//		time_gap_i = abs(time1[t,1] :- time1)
//		time_gap_cut_i = time_gap_i :<= lagcutoff_s
//		weight = (1 :- time_gap_i :/ (lagcutoff_s +1)):* (time_gap_i :<= lagcutoff_s)
//		w[s..e,1] = weight
//	}
//	else {
//	w[s..e,1] = J(t1,1,0)

	w[s..e,1] = time1:==time[pos,1]
	*w[pos,1] = 1

//	}
	// __ VCV
		wz = w :> 0
		nonzero = nonzero + sum(wz)
		V1 =  P * (res :* res[pos,1] * w[.,1]) * (P[.,pos])'
		V = V + V1
		AVCVb_fs_1 = (I(Nend)#SzN) * ((vec(u) * u[pos,.]) :* (J(Nend,1, 1)# w[.,1])) * (I(Nend)#SzN[.,pos])' // -------------------------- KP
		AVCVb_fs =  AVCVb_fs + AVCVb_fs_1																	 // -------------------------- KP

**adjj[.,ind+1] = w[.,1] // check
	}

	ind = ind + t1
}
}

_makesymmetric(V)

}
end





























*********************************************************************************************

capt program drop  acreg_ranktest

* ranktest 1.4.01  18aug2015
* author mes, based on code by fk
* see end of file for version comments

*if c(version) < 12 {
* ranktest uses livreg2 Mata library.
* Ensure Mata library is indexed if new install.
* Not needed for Stata 12+ since ssc.ado does this when installing.
	capture mata: mata drop m_calckw()
	capture mata: mata drop m_omega()
	capture mata: mata drop ms_vcvorthog()
	capture mata: mata drop s_vkernel()
	capture mata: mata drop  acreg_s_rkstat()
	capture mata: mata drop acreg_cholqrsolve()
	mata: mata mlib index
*}

program define acreg_ranktest, rclass sortpreserve

	local lversion 01.4.01

*	if _caller() < 11 {
*		ranktest9 `0'
*		return add						//  otherwise all the ranktest9 results are zapped
*		return local ranktestcmd		ranktest9
*		return local cmd				ranktest
*		return local version			`lversion'
*		exit
*	}
	version 11.2

*	if substr("`1'",1,1)== "," {
*		if "`2'"=="version" {
*			di in ye "`lversion'"
*			return local version `lversion'
*			exit
*		}
*		else {
* 	di as err "invalid syntax"
*			exit 198
*		}
*	}

* If varlist 1 or varlist 2 have a single element, parentheses optional

	if substr("`1'",1,1)=="(" {
		acreg_GetVarlist `0'
		local y `s(varlist)'
		local 0 `"`s(rest)'"'
		sret clear
	}
	else {
		local y `1'
		mac shift 1
		local 0 `"`*'"'
	}

	if substr("`1'",1,1)=="(" {
		acreg_GetVarlist `0'
		local z `s(varlist)'
		local 0 `"`s(rest)'"'
		sret clear
	}
	else {
		local z `1'
		mac shift 1
* Need to reinsert comma before options (if any) for -syntax- command to work
		local 0 `", `*'"'
	}

// Note that y or z could be a varlist, e.g., "y1-y3", so they need to be unab-ed.
	tsunab y : `y'
	local K : word count `y'
	tsunab z : `z'
	local L : word count `z'

* Option version ignored here if varlists were provided
	syntax [if] [in] [aw fw pw iw/]				///
		[,										///
		partial(varlist ts)						///
		fwl(varlist ts)							///
		NOConstant								///
		wald									///
		ALLrank									///
		NULLrank								///
		FULLrank								///
		ROBust									///
		cluster(varlist)						///
		BW(string)								///
		kernel(string)							///
		Tvar(varname)							///
		Ivar(varname)							///
		sw										///
		psd0									///
		psda									///
		version									///
		dofminus(integer 0)						///
		mat_clst(string)						///
		]






	local partial		"`partial' `fwl'"
	local partial		: list retokenize partial

	local cons		= ("`noconstant'"=="")

	mat mat_clst=`mat_clst'


	if "`wald'"~="" {
		local LMWald "Wald"
	}
	else {
		local LMWald "LM"
	}

	local optct : word count `allrank' `nullrank' `fullrank'
	if `optct' > 1 {
di as err "Incompatible options: `allrank' `nullrank' `fullrank'"
		error 198
	}
	else if `optct' == 0 {
* Default
		local allrank "allrank"
	}

	local optct : word count `psd0' `psda'
	if `optct' > 1 {
di as err "Incompatible options: `psd0' `psda'"
		error 198
	}
	local psd	"`psd0' `psda'"
	local psd	: list retokenize psd

* Note that by tsrevar-ing here, subsequent disruption to the sort doesn't matter
* for TS operators.
	tsrevar `y'
	local vl1 `r(varlist)'
	tsrevar `z'
	local vl2 `r(varlist)'
	tsrevar `partial'
	local partial `r(varlist)'

	foreach vn of varlist `vl1' {
		tempvar tv
		qui gen double `tv' = .
		local tempvl1 "`tempvl1' `tv'"
	}
	foreach vn of varlist `vl2' {
		tempvar tv
		qui gen double `tv' = .
		local tempvl2 "`tempvl2' `tv'"
	}

	marksample touse
	markout `touse' `vl1' `vl2' `partial' `cluster', strok

* Stock-Watson and cluster imply robust.
	if "`sw'`cluster'" ~= "" {
		local robust "robust"
	}

	tempvar wvar
	if "`weight'" == "fweight" | "`weight'"=="aweight" {
		local wtexp `"[`weight'=`exp']"'
		gen double `wvar'=`exp'
	}
	if "`fsqrt(wf)*(wvar^0.5):*'" == "fweight" & "`kernel'" !="" {
		di in red "fweights not allowed (data are -tsset-)"
		exit 101
	}
	if "`weight'" == "fweight" & "`sw'" != "" {
		di in red "fweights currently not supported with -sw- option"
		exit 101
	}
	if "`weight'" == "iweight" {
		if "`robust'`cluster'`bw'" !="" {
			di in red "iweights not allowed with robust, cluster, AC or HAC"
			exit 101
		}
		else {
			local wtexp `"[`weight'=`exp']"'
			gen double `wvar'=`exp'
		}
	}
	if "`weight'" == "pweight" {
		local wtexp `"[aweight=`exp']"'
		gen double `wvar'=`exp'
		local robust "robust"
	}
	if "`weight'" == "" {
* If no weights, define neutral weight variable
		qui gen byte `wvar'=1
	}


* Every time a weight is used, must multiply by scalar wf ("weight factor")
* wf=1 for no weights, fw and iw, wf = scalar that normalizes sum to be N if aw or pw
		sum `wvar' if `touse' `wtexp', meanonly
* Weight statement
		if "`weight'" ~= "" {
di in gr "(sum of wgt is " %14.4e `r(sum_w)' ")"
		}
		if "`weight'"=="" | "`weight'"=="fweight" | "`weight'"=="iweight" {
* If weight is "", weight var must be column of ones and N is number of rows.
* With fw and iw, effective number of observations is sum of weight variable.
			local wf=1
			local N=r(sum_w)
		}
		else if "`weight'"=="aweight" | "`weight'"=="pweight" {
* With aw and pw, N is number of obs, unadjusted.
			local wf=r(N)/r(sum_w)
			local N=r(N)
		}
		else {
* Should never reach here
di as err "ranktest error - misspecified weights"
			exit 198
		}

* HAC estimation.
* If bw is omitted, default `bw' is empty string.
* If bw or kernel supplied, check/set `kernel'.
* Macro `kernel' is also used for indicating HAC in use.
	if "`bw'" == "" & "`kernel'" == "" {
		local bw=0
	}
	else {
* Need tvar for markout with time-series stuff
* Data must be tsset for time-series operators in code to work
* User-supplied tvar checked if consistent with tsset
		capture tsset
		if "`r(timevar)'" == "" {
di as err "must tsset data and specify timevar"
			exit 5
		}
		if "`tvar'" == "" {
			local tvar "`r(timevar)'"
		}
		else if "`tvar'"!="`r(timevar)'" {
di as err "invalid tvar() option - data already -tsset-"
			exit 5
		}
* If no panel data, ivar will still be empty
		if "`ivar'" == "" {
			local ivar "`r(panelvar)'"
		}
		else if "`ivar'"!="`r(panelvar)'" {
di as err "invalid ivar() option - data already -tsset-"
			exit 5
		}
		local tdelta `r(tdelta)'
		tsreport if `touse', panel
		if `r(N_gaps)' != 0 {
di in gr "Warning: time variable " in ye "`tvar'" in gr " has " /*
	*/ in ye "`r(N_gaps)'" in gr " gap(s) in relevant range"
		}

* Check it's a valid kernel and replace with unabbreviated kernel name; check bw.
* Automatic kernel selection allowed by ivreg2 but not ranktest so must trap.
* s_vkernel is in livreg2 mlib.
		if "`bw'"=="auto" {
di as err "invalid bandwidth in option bw() - must be real > 0"
			exit 198
		}
		mata: s_vkernel("`kernel'", "`bw'", "`ivar'")
		local kernel `r(kernel)'
		local bw = `r(bw)'
	}

* tdelta missing if version 9 or if not tsset
	if "`tdelta'"=="" {
		local tdelta=1
	}

	if "`sw'"~="" {
		capture xtset
		if "`ivar'" == "" {
			local ivar "`r(panelvar)'"
		}
		else if "`ivar'"!="`r(panelvar)'" {
di as err "invalid ivar() option - data already tsset or xtset"
			exit 5
		}
* Exit with error if ivar is neither supplied nor tsset nor xtset
		if "`ivar'"=="" {
di as err "Must -xtset- or -tsset- data or specify -ivar- with -sw- option"
			exit 198
		}
		qui describe, short varlist
		local sortlist "`r(sortlist)'"
		tokenize `sortlist'
		if "`ivar'"~="`1'" {
di as err "Error - dataset must be sorted on panel var with -sw- option"
			exit 198
		}
	}

* Create variable used for getting lags etc. in Mata
	tempvar tindex
	qui gen `tindex'=1 if `touse'
	qui replace `tindex'=sum(`tindex') if `touse'

********** CLUSTER SETUP **********************************************

* Mata code requires data are sorted on (1) the first var cluster if there
* is only one cluster var; (2) on the 3rd and then 1st if two-way clustering,
* unless (3) two-way clustering is combined with kernel option, in which case
* the data are tsset and sorted on panel id (first cluster variable) and time
* id (second cluster variable).
* Second cluster var is optional and requires an identifier numbered 1..N_clust2,
* unless combined with kernel option, in which case it's the time variable.
* Third cluster var is the intersection of 1 and 2, unless combined with kernel
* opt, in which case it's unnecessary.
* Sorting on "cluster3 cluster1" means that in Mata, panelsetup works for
* both, since cluster1 nests cluster3.
* Note that it is possible to cluster on time but not panel, in which case
* cluster1 is time, cluster2 is empty and data are sorted on panel-time.
* Note also that if data are sorted here but happen to be tsset, will need
* to be re-tsset after estimation code concludes.


// No cluster options or only 1-way clustering
// but for Mata and other purposes, set N_clust vars =0
	local N_clust=0
	local N_clust1=0
	local N_clust2=0
	if "`cluster'"!="" {
		local clopt "cluster(`cluster')"
		tokenize `cluster'
		local cluster1 "`1'"
		local cluster2 "`2'"
		if "`kernel'"~="" {
* kernel requires either that cluster1 is time var and cluster2 is empty
* or that cluster1 is panel var and cluster2 is time var.
* Either way, data must be tsset and sorted for panel data.
			if "`cluster2'"~="" {
* Allow backwards order
				if "`cluster1'"=="`tvar'" & "`cluster2'"=="`ivar'" {
					local cluster1 "`2'"
					local cluster2 "`1'"
				}
				if "`cluster1'"~="`ivar'" | "`cluster2'"~="`tvar'" {
di as err "Error: cluster kernel-robust requires clustering on tsset panel & time vars."
di as err "       tsset panel var=`ivar'; tsset time var=`tvar'; cluster vars=`cluster1',`cluster2'"
					exit 198
				}
			}
			else {
				if "`cluster1'"~="`tvar'" {
di as err "Error: cluster kernel-robust requires clustering on tsset time variable."
di as err "       tsset time var=`tvar'; cluster var=`cluster1'"
					exit 198
				}
			}
		}
* Simple way to get quick count of 1st cluster variable without disrupting sort
* clusterid1 is numbered 1.._Nclust1.
		tempvar clusterid1
		qui egen `clusterid1'=group(`cluster1') if `touse'
		sum `clusterid1' if `touse', meanonly
		if "`cluster2'"=="" {
			local N_clust=r(max)
			local N_clust1=`N_clust'
			if "`kernel'"=="" {
* Single level of clustering and no kernel-robust, so sort on single cluster var.
* kernel-robust already sorted via tsset.
				sort `cluster1'
			}
		}
		else {
			local N_clust1=r(max)
			if "`kernel'"=="" {
				tempvar clusterid2 clusterid3
* New cluster id vars are numbered 1..N_clust2 and 1..N_clust3
				qui egen `clusterid2'=group(`cluster2') if `touse'
				qui egen `clusterid3'=group(`cluster1' `cluster2') if `touse'
* Two levels of clustering and no kernel-robust, so sort on cluster3/nested in/cluster1
* kernel-robust already sorted via tsset.
				sort `clusterid3' `cluster1'
				sum `clusterid2' if `touse', meanonly
				local N_clust2=r(max)
			}
			else {
* Need to create this only to count the number of clusters
				tempvar clusterid2
				qui egen `clusterid2'=group(`cluster2') if `touse'
				sum `clusterid2' if `touse', meanonly
				local N_clust2=r(max)
* Now replace with original variable
				local clusterid2 `cluster2'
			}

			local N_clust=min(`N_clust1',`N_clust2')

		}		// end 2-way cluster block
	}		// end cluster block

************************************************************************************************

* Note that bw is passed as a value, not as a string

	mata: acreg_s_rkstat(						///
					"`vl1'",			///
					"`vl2'",			///
					"`partial'",		///
					"`wvar'",			///
					"`weight'",			///
					`wf',				///
					`N',				///
					`cons',				///
					"`touse'",			///
					"`LMWald'",			///
					"`allrank'",		///
					"`nullrank'",		///
					"`fullrank'",		///
					"`robust'",			///
					"`clusterid1'",		///
					"`clusterid2'",		///
					"`clusterid3'",		///
					`bw',				///
					"`tvar'",			///
					"`ivar'",			///
					"`tindex'",			///
					`tdelta',			///
					`dofminus',			///
					"`kernel'",			///
					"`sw'",				///
					"`psd'",			///
					"`tempvl1'",		///
					"`tempvl2'",		///
					"`mat_clst'"		///
					)

	tempname rkmatrix chi2 df df_r p rank ccorr eval mat_clst
	mat `rkmatrix'=r(rkmatrix)
	mat `ccorr'=r(ccorr)
	mat `eval'=r(eval)
	mat colnames `rkmatrix' = "rk" "df" "p" "rank" "eval" "ccorr"
	mat `mat_clst'=r(mat_clst)

di
di "Kleibergen-Paap rk `LMWald' test of rank of matrix"
	if "`robust'"~="" & "`kernel'"~= "" & "`cluster'"=="" {
di "  Test statistic robust to heteroskedasticity and autocorrelation"
di "  Kernel: `kernel'   Bandwidth: `bw'"
	}
	else if "`kernel'"~="" & "`cluster'"=="" {
di "  Test statistic robust to autocorrelation"
di "  Kernel: `kernel'   Bandwidth: `bw'"
	}
	else if "`cluster'"~="" {
di "  Test statistic robust to heteroskedasticity and clustering on `cluster'"
		if "`kernel'"~="" {
di "  and kernel-robust to common correlated disturbances"
di "  Kernel: `kernel'   Bandwidth: `bw'"
		}
	}
	else if "`robust'"~="" {
di "  Test statistic robust to heteroskedasticity"
	}
	else if "`LMWald'"=="LM" {
di "  Test assumes homoskedasticity (Anderson canonical correlations test)"
	}
	else {
di "  Test assumes homoskedasticity (Cragg-Donald test)"
	}

	local numtests = rowsof(`rkmatrix')
	forvalues i=1(1)`numtests' {
di "Test of rank=" %3.0f `rkmatrix'[`i',4] "  rk=" %8.2f `rkmatrix'[`i',1] /*
	*/	"  Chi-sq(" %3.0f `rkmatrix'[`i',2] ") pvalue=" %8.6f `rkmatrix'[`i',3]
	}
	scalar `chi2' = `rkmatrix'[`numtests',1]
	scalar `p' = `rkmatrix'[`numtests',3]
	scalar `df' = `rkmatrix'[`numtests',2]
	scalar `rank' = `rkmatrix'[`numtests',4]
	local N `r(N)'
	return scalar df = `df'
	return scalar chi2 = `chi2'
	return scalar p = `p'
	return scalar rank = `rank'
	if "`cluster'"~="" {
		return scalar N_clust = `N_clust'
	}
	if "`cluster2'"~="" {
		return scalar N_clust1 = `N_clust1'
		return scalar N_clust2 = `N_clust2'
	}
	return scalar N = `N'
	return matrix rkmatrix `rkmatrix'
	return matrix ccorr `ccorr'
	return matrix eval `eval'

	return matrix mat_clst `mat_clst'

	tempname S V Omega
	if `K' > 1 {
		foreach en of local y {
* Remove "." from equation name
			local en1 : subinstr local en "." "_", all
			foreach vn of local z {
				local cn "`cn' `en1':`vn'"
			}
		}
	}
	else {
		foreach vn of local z {
		local cn "`cn' `vn'"
		}
	}

	mat `V'=r(V)
	matrix colnames `V' = `cn'
	matrix rownames `V' = `cn'
	return matrix V `V'
	mat `S'=r(S)
	matrix colnames `S' = `cn'
	matrix rownames `S' = `cn'
	return matrix S `S'

	return local cmd		"ranktest"
	return local version	`lversion'
end


* Adopted from -canon-
capt program drop acreg_GetVarlist
program define acreg_GetVarlist, sclass
	version 11.2
	sret clear
	gettoken open 0 : 0, parse("(")
	if `"`open'"' != "(" {
		error 198
	}
	gettoken next 0 : 0, parse(")")
	while `"`next'"' != ")" {
		if `"`next'"'=="" {
			error 198
		}
		local list `list'`next'
		gettoken next 0 : 0, parse(")")
	}
	sret local rest `"`0'"'
	tokenize `list'
	local 0 `*'
	sret local varlist "`0'"
end


********************* EXIT IF STATA VERSION < 11 ********************************

* When do file is loaded, exit here if Stata version calling program is < 11.
* Prevents loading of rest of program file (would cause e.g. Stata 10 to crash at Mata).

if c(stata_version) < 11 {
	exit
}

******************** END EXIT IF STATA VERSION < 9 *****************************

*******************************************************************************
*************************** BEGIN MATA CODE ***********************************
*******************************************************************************

version 11.2
mata:

// ********* MATA CODE SHARED BY ivreg2 AND ranktest       *************** //
// ********* 1. struct ms_vcvorthog                        *************** //
// ********* 2. m_omega                                    *************** //
// ********* 3. m_calckw                                   *************** //
// ********* 4. s_vkernel                                  *************** //
// *********************************************************************** //

// For reference:
 struct ms_vcvorthog {
 	string scalar	ename, Znames, touse, weight, wvarname
 	string scalar	robust, clustvarname, clustvarname2, clustvarname3, kernel
 	string scalar	sw, psd, ivarname, tvarname, tindexname
 	real scalar		wf, N, bw, tdelta, dofminus
    real scalar		center
 	real matrix		ZZ
 	pointer matrix	e
 	pointer matrix	Z
 	pointer matrix	wvar
 }

void acreg_s_rkstat(	string scalar vl1,
				string scalar vl2,
				string scalar partial,
				string scalar wvarname,
				string scalar weight,
				scalar wf,
				scalar N,
				scalar cons,
				string scalar touse,
				string scalar LMWald,
				string scalar allrank,
				string scalar nullrank,
				string scalar fullrank,
				string scalar robust,
				string scalar clustvarname,
				string scalar clustvarname2,
				string scalar clustvarname3,
				bw,
				string scalar tvarname,
				string scalar ivarname,
				string scalar tindexname,
				tdelta,
				dofminus,
				string scalar kernel,
				string scalar sw,
				string scalar psd,
				string scalar tempvl1,
				string scalar tempvl2,
				string scalar mat_clst)
{

// iid flag used below
	iid = ((kernel=="") & (robust=="") & (clustvarname==""))

// tempx, tempy and tempz are the Stata names of temporary variables that will be changed by acreg_s_rkstat
	tempy=tokens(tempvl1)
	tempz=tokens(tempvl2)
	tempx=tokens(partial)

	st_view(y=.,.,tokens(vl1),touse)
	st_view(z=.,.,tokens(vl2),touse)
	st_view(yhat=.,.,tempy,touse)
	st_view(zhat=.,.,tempz,touse)
	if (partial~="") {
		st_view(x=.,.,tempx,touse)
	}
	st_view(mtouse=.,.,tokens(touse),touse)
	st_view(wvar=.,.,tokens(wvarname),touse)
	noweight=(st_vartype(wvarname)=="byte")

	K=cols(y)							//  count of vars in first varlist
	L=cols(z)							//  count of vars in second varlist
	P=cols(x)							//  count of vars to be partialled out (excluding constant)

// Note that we now use wf*wvar instead of wvar
// because wvar is raw weighting variable and
// wf*wvar normalizes so that sum(wf*wvar)=N.

// Partial out the X variables.
// Note that this includes demeaning if there is a constant,
//   i.e., variables are centered.
	if (cons & P>0) {					//  Vars to partial out including constant
		ymeans = mean(y,wf*wvar)
		zmeans = mean(z,wf*wvar)
		xmeans = mean(x,wf*wvar)
		xy = quadcrossdev(x, xmeans, wf*wvar, y, ymeans)
		xz = quadcrossdev(x, xmeans, wf*wvar, z, zmeans)
		xx = quadcrossdev(x, xmeans, wf*wvar, x, xmeans)
	}
	else if (!cons & P>0) {				//  Vars to partial out NOT including constant
		xy = quadcross(x, wf*wvar, y)
		xz = quadcross(x, wf*wvar, z)
		xx = quadcross(x, wf*wvar, x)
	}
	else {								//  Only constant to partial out = demean
		ymeans = mean(y,wf*wvar)
		zmeans = mean(z,wf*wvar)
	}
//	Partial-out coeffs. Default Cholesky; use QR if not full rank and collinearities present.
//	Not necessary if no vars other than constant
	if (P>0) {
		by = acreg_cholqrsolve(xx, xy)
		bz = acreg_cholqrsolve(xx, xz)
	}
//	Replace with residuals
	if (cons & P>0) {					//  Vars to partial out including constant
		yhat[.,.] = (y :- ymeans) - (x :- xmeans)*by
		zhat[.,.] = (z :- zmeans) - (x :- xmeans)*bz
	}
	else if (!cons & P>0) {				//  Vars to partial out NOT including constant
		yhat[.,.] = y - x*by
		zhat[.,.] = z - x*bz
	}
	else if (cons) {					//  Only constant to partial out = demean
		yhat[.,.] = (y :- ymeans)
		zhat[.,.] = (z :- zmeans)
	}
	else {								//  no transformations required
		yhat[.,.] = y
		zhat[.,.] = z
	}

	zhzh = quadcross(zhat, wf*wvar, zhat)
	zhyh = quadcross(zhat, wf*wvar, yhat)
	yhyh = quadcross(yhat, wf*wvar, yhat)

//	pihat = invsym(zhzh)*zhyh
	pihat = acreg_cholqrsolve(zhzh, zhyh)

// rzhat is F in paper (p. 103)
// iryhat is G in paper (p. 103)
	ryhat=cholesky(yhyh)
	rzhat=cholesky(zhzh)
	iryhat=luinv(ryhat')
	irzhat=luinv(rzhat')
	that=rzhat'*pihat*iryhat

// cc is canonical correlations.  Squared cc is eigenvalues.
	fullsvd(that, ut, cc, vt)
	vt=vt'
	vecth=vec(that)
	ev = cc:^2
// S matrix in paper (p. 100).  Not used in code below.
//	smat=fullsdiag(cc, rows(that)-cols(that))

	if (abs(1-cc[1,1])<1e-10) {
printf("\n{text:Warning: collinearities detected between (varlist1) and (varlist2)}\n")
	}
	if ((missing(ryhat)>0) | (missing(iryhat)>0) | (missing(rzhat)>0) | (missing(irzhat)>0)) {
printf("\n{error:Error: non-positive-definite matrix. May be caused by collinearities.}\n")
		exit(error(3351))
	}

// If Wald, yhat is residuals
	if (LMWald=="Wald") {
		yhat[.,.]=yhat-zhat*pihat
		yhyh = quadcross(yhat, wvar, yhat)
	}

// Covariance matrices
// vhat is W in paper (eqn below equation 17, p. 103)
// shat is V in paper (eqn below eqn 15, p. 103)



// ************************************************************************************* //
// shat calculated using struct and programs m_omega, m_calckw shared with ivreg2        //

	struct ms_vcvorthog scalar vcvo


	vcvo.ename			= tempy		// ivreg2 has = ename //
	vcvo.Znames			= tempz		// ivreg2 has = Znames //
	vcvo.touse			= touse
	vcvo.weight			= weight
	vcvo.wvarname		= wvarname
	vcvo.robust			= robust
	vcvo.clustvarname	= clustvarname
	vcvo.clustvarname2	= clustvarname2
	vcvo.clustvarname3	= clustvarname3
	vcvo.kernel			= kernel
	vcvo.sw				= sw
	vcvo.psd			= psd
	vcvo.ivarname		= ivarname
	vcvo.tvarname		= tvarname
	vcvo.tindexname		= tindexname
	vcvo.wf				= wf
	vcvo.N				= N
	vcvo.bw				= bw
	vcvo.tdelta			= tdelta
	vcvo.dofminus		= dofminus
	vcvo.ZZ				= zhzh		// ivreg2 has = st_matrix(ZZmatrix) //

	vcvo.e		= &yhat				// ivreg2 has = &e	//
	vcvo.Z		= &zhat				// ivreg2 has = &Z //
	vcvo.wvar	= &wvar

	shat=m_omega(vcvo)

	mat_clst=st_matrix(mat_clst)

	shat = mat_clst

// ***************************************************************************************

// prepare to start collecting test stats
	if (allrank~="") {
		firstrank=1
		lastrank=min((K,L))
	}
	else if (nullrank~="") {
		firstrank=1
		lastrank=1
	}
	else if (fullrank~="") {
		firstrank=min((K,L))
		lastrank=min((K,L))
	}
	else {
// should never reach this point
printf("ranktest error\n")
		exit
	}

// where results will go
	rkmatrix=J(lastrank-firstrank+1,6,.)

// ***************************************************************************************
// Calculate vector of canonical correlations test statistics.
// All we need if iid case.
	rkvec = ev									//  Initialize vector with individual eigenvalues.
	if (LMWald~="LM") {							//  LM is sum of min evals, Wald is sum of eval/(1-eval)
		rkvec = rkvec :/ (1 :- rkvec)
	}
	for (i=(rows(rkvec)-1); i>=1; i--) {		//  Now loop through and sum the eigenvalues.
		rkvec[i,1] = rkvec[i+1,1] + rkvec[i,1]
	}
	rkvec = N*rkvec								//  Multiply by N to get the test statistics.

// ***************************************************************************************

// Finally, calcluate vhat
	if ((LMWald=="LM") & (iid)) {
// Homoskedastic, iid LM case means vcv is identity matrix
// Generates canonical correlation stats.  Default.
		vhat=I(L*K,L*K)/N
	}
	else {
		vhat=(iryhat'#irzhat')*shat*(iryhat'#irzhat')' * N
		_makesymmetric(vhat)
// Homoskedastic iid Wald case means vcv has block-diag identity matrix structure.
// Enforce this by setting ~0 entries to 0.  If iid, vhat not used in calcs, for reporting only.
		if ((LMWald=="Wald") & (iid)) {
			vhat = vhat :* (J(K,K,1)#I(L))
		}
	}

// ***************************************************************************************
// Loop through ranks and collect test stats, dfs, p-values, ranks, evs and ev^2 (=ccs)

	for (i=firstrank; i<=lastrank; i++) {
		if (iid) {							//  iid case = canonical correlations test
			rk = rkvec[i,1]
			}
		else {								//  non-iid case
			if (i>1) {
				u12=ut[(1::i-1),(i..L)]
				v12=vt[(1::i-1),(i..K)]
			}
			u22=ut[(i::L),(i..L)]
			v22=vt[(i::K),(i..K)]

			symeigensystem(u22*u22', evec, eval)
			u22v=evec
			u22d=diag(eval)
			u22h=u22v*(u22d:^0.5)*u22v'

			symeigensystem(v22*v22', evec, eval)
			v22v=evec
			v22d=diag(eval)
			v22h=v22v*(v22d:^0.5)*v22v'

			if (i>1) {
				aq=(u12 \ u22)*luinv(u22)*u22h
				bq=v22h*luinv(v22')*(v12 \ v22)'
			}
			else {
				aq=u22*luinv(u22)*u22h
				bq=v22h*luinv(v22')*v22'
			}

// lab is lambda_q in paper (eqn below equation 21, p. 104)
// vlab is omega_q in paper (eqn 19 in paper, p. 104)
			lab=(bq#aq')*vecth
			vlab=(bq#aq')*vhat*(bq#aq')'

// Symmetrize if numerical inaccuracy means it isn't
			_makesymmetric(vlab)
			vlabinv=invsym(vlab)
// rk stat Assumption 2: vlab (omega_q in paper) is nonsingular.  Detected by a zero on the diagonal,
// since when returning a generalized inverse, Stata/Mata choose the generalized inverse that
// sets entire column(s)/row(s) to zeros.
			if (diag0cnt(vlabinv)>0) {
				rk = .
printf("\n{text:Warning: covariance matrix omega_%f}", i-1)
printf("{text: not full rank; test of rank %f}", i-1)
printf("{text: unavailable}\n")
			}
// Note not multiplying by N - already incorporated in vhat.
			else {
				rk=lab'*vlabinv*lab
			}
		}												//  end non-iid case
// at this point rk has value of test stat
// fill out rest of row of rkmatrix
// save df, rank, etc. even if test stat not available.
		df=(L-i+1)*(K-i+1)
		pvalue=chi2tail(df, rk)
		rkmatrix[i-firstrank+1,1]=rk
		rkmatrix[i-firstrank+1,2]=df
		rkmatrix[i-firstrank+1,3]=pvalue
		rkmatrix[i-firstrank+1,4]=i-1
		rkmatrix[i-firstrank+1,5]=ev[i-firstrank+1,1]
		rkmatrix[i-firstrank+1,6]=cc[i-firstrank+1,1]
// end of test loop
	}

// ***************************************************************************************
// Finish up and return results

	st_matrix("r(rkmatrix)", rkmatrix)
	st_matrix("r(ccorr)", cc')
	st_matrix("r(eval)",ev')
// Save V matrix as in paper, without factor of 1/N
	vhat=N*vhat*wf
	st_matrix("r(V)", vhat)
// Save S matrix as in ivreg2, with factor of 1/N
	st_matrix("r(S)", shat)
	st_numscalar("r(N)", N)
	if (clustvarname~="") {
		st_numscalar("r(N_clust)", N_clust)
	}
	if (clustvarname2~="") {
		st_numscalar("r(N_clust2)", N_clust2)
	}
// end of program
	st_matrix("r(mat_clst)", mat_clst)



}

// Mata utility for sequential use of solvers
// Default is cholesky;
// if that fails, use QR;
// if overridden, use QR.

function acreg_cholqrsolve (	numeric matrix A,
						numeric matrix B,
						| real scalar useqr)
{
	if (args()==2) useqr = 0

	real matrix C

	if (!useqr) {
		C = cholsolve(A, B)
		if (C[1,1]==.) {
			C = qrsolve(A, B)
		}
	}
	else {
		C = qrsolve(A, B)
	}

	return(C)

}

end