*! version 1.0.4 20sept2019
		
		program define rho_xtregar, rclass
		version 10
		tempvar K Ki prem NbrindvN ni lemaxni num yhat utilde tempdnum0 /*
			*/ tempdnum1 dnum tempddenom1 ddnom WIndiv1 WIndiv2
        tempname tdelta final Windiv1 Windiv2
		
        syntax varlist(fv ts) [if] [using/] [, approx approx_balanced approx_unbalanced option_nodisplay]
		
				
		if "`weight'" != "" {
                local weights "weights"
                        di as error "no option for weights so far"
                        exit 198
        }
        else {
                local weights ""
        }
	
		
		//The preserve is not done for the if condition in itself
		quietly: preserve		
		if "`using'"!="" {
			quietly: use `using', clear
		}	
		if "`if'"!="" {
			quietly: keep `if'
		}
				
		 _xt, trequired
        local id "`r(ivar)'"
        local t  "`r(tvar)'"
        scalar `tdelta' = r(tdelta)
		
		gettoken depvar xvars : varlist

		
		quietly {
		
		
		*****************************
		//definition of the Ki
			bysort `id': gen K = (`t'[_n-1]==`t'[_n]-`tdelta') if  _n!=1
			bysort `id': egen Ki=total(K)
			bysort `id': gen prem = (_n==1) //premier means first in French
		}
		
		if "`option_nodisplay'"=="" {
			quietly {
			egen WIndiv1 = total(prem)
			local Windiv1 = WIndiv1[1]
			egen WIndiv2 = total(prem*(Ki>0 & Ki!=.))
			local Windiv2 = WIndiv2[1]
			}
		dis "There are `Windiv1' units, over which `Windiv2' have at least two successive observations"
		}
		//
		
		quietly {
		*****************************
		//erasing those that never have two successive observations		
		drop if Ki==0 
		
		
		//Total Number of individuals
			egen double NbrindvN= total(prem)
			local the_NbrindvN = NbrindvN[1]
		//sometimes useful: ni and maxni 
		bysort `id': gen ni = _N
		egen double lemaxni = max(ni)
			local the_maxni = lemaxni[1]
		//and to get vecN2:=(n_1^2, (...), n_N^2):
			gsort -prem +`id'
			order ni
			mata: vecN1 = st_data((1,`the_NbrindvN'), 1) 
			mata: vecN2=power(2,vecN1)  
			sort `id' `t'
		
		*****************************
		*term at the numerator (save 1-rho)
		egen double num = total(prem * Ki /(1+Ki))
		replace num =  num / `the_NbrindvN'
		local the_num = num[1]
		
		*****************************
		*estimate of r_d		
		//THESE TWO LINES SHOULD BE SPEED UP
		xtreg `depvar' `xvars', fe	
		predict yhat, xbu 
		gen utilde = `depvar' - yhat	
		//the numerateur	
			gen double tempdnum0 = 0
			bysort `id': replace tempdnum0 = ((utilde[_n]-utilde[_n-1])^2) ///
						if _n!=1 & (`t'[_n-1]==`t'[_n]-`tdelta')
			bysort `id': egen tempdnum1 = total(tempdnum0)
			bysort `id': replace tempdnum1 = tempdnum1 * prem / (1+Ki)
			egen double dnum = total(tempdnum1)
			local the_dnum = dnum[1]
		//the denominateor
			gen double tempddenom0 = utilde*utilde	
				bysort `id': egen tempddenom1 = total(tempddenom0)
				bysort `id': replace tempddenom1 = tempddenom1 * prem / ni
			egen double ddnom = total(tempddenom1)
			local the_ddnom = ddnom[1]
		// final: r_d = 1-d/2
		local the_rhod= 	1-0.5*`the_dnum' / `the_ddnom'
				
		//the general case
			if ("`approx'"=="") & ("`approx_unbalanced'"=="") & ("`approx_balanced'"=="") {
				
				***********
				//estimate of A, mata matrix of size(`the_NbrindvN',`the_maxni'*`the_maxni')
				//whith the exponents for the rho	
				keep `id' `t' ni lemaxni
				forvalues i = 1/`the_maxni' {
					forvalues j = 1/`the_maxni' {
					quietly{
						gen double coeff_`i'_`j' = .
							bysort `id': replace coeff_`i'_`j' = abs(`t'[`i']-`t'[`j'])/`tdelta' if max(`i',`j')<=ni
					}
					//
					}
					//
				}
				//
				bysort `id': keep if _n==1
				drop `id' `t' ni lemaxni
				mata: A = st_data(., .)	
				**********
				//almost done:
				//definition of MATA function pourE below the main program
			
				**********
				//preparing the optimisation
					mata: myrhoD = `the_rhod'
					mata: mynum = `the_num'
					mata: myNbrindvN  = `the_NbrindvN'	
				//definition of MATA function funcfinal below the main program
				//the optimisation
				mata: mm_root(r=.,&funcfinal(), -1.1, 1.1, 0, 1000,myrhoD,mynum,myNbrindvN,A,vecN2) 
				mata: r
				mata: st_numscalar("final",r)
				restore
			}
			//
		
		
		//with option approx 
			if ("`approx'"!="") | ("`approx_unbalanced'"!="") | ("`approx_balanced'"!="") {			
				restore 
				preserve		
				
					if "`using'"!="" {
						quietly: use `using', clear
					}	
					if "`if'"!="" {
						quietly: keep `if'
					}
				
				xtset 
				local balancing = r(balanced)
				local tmax=r(tmaxs)
				local tmin=r(tmins)
				local thegrandT = `tmax'-`tmin'
			
				//balanced case
				if inlist("`balancing'","weakly balanced","strongly balanced") | "`approx_balanced'"!="" {
					scalar final = `the_rhod'/(1-2/`thegrandT')
				}
				//	
			
				//unbalanced case
				if inlist("`balancing'","unbalanced") | "`approx_unbalanced'"!="" {	
						scalar final = (`the_rhod'-1+`the_num')/`the_num'			
				}
				//
			restore
			}
			//
		
		}
		//end of the quietly
			
		dis "*******   rho_BFN  =  " final
		
		return scalar rho_BFN = final 
		end
	

		
		//MATA functions used in the program
		mata: function power(a,r) return (exp(ln(r)*a))		
		mata: function pourE(A,r,vecN2,N) return (sum(power(A,r):/vecN2)/N)
		mata: function funcfinal(r,rhoD,lenum,leNbrindvN,A,vecN2) ///
			return (rhoD-1+((1-r)*lenum)/(1-pourE(A,r,vecN2,leNbrindvN)))