** This version implements the Loops In MATA!!
capture program drop vc_test
program vc_test, rclass sortpreserve
syntax varlist(numeric fv) [if] [in], [ vcoeff(varname) bw(real -1) seed(str) ///
									knots(real -1) km(real 1) vcoeff_par(varname) kernel(str) wbsrep(real 50) degree(real 0) ]

	tokenize `varlist'
	local y `1'
	//get the rest of the vars
	macro shift
	local varx `*'
	if "`seed'"!="" set seed `seed'
	** Defines Sample, Works with and in addition to IF
	marksample touse
	markout `touse' `sample' `varlist' `vcoeff'
	sort `touse'
	**
	*********************************************
	** checking if vcoeff was defined or using "system"
	if "`vcoeff'"=="" & "$vcoeff_"=="" {
	   display "Need to define vcoeff. variable for varying coefficients"
	   exit
	}
	*** this uses system
	if "`vcoeff'"=="" & "$vcoeff_"!="" {
		local vcoeff $vcoeff_
    }
	local full_x `vcoeff'

    ** Defining the running variable. Or the Bins
	if `knots'>=0 {
		tempvar _kn
		egen `_kn'=vbin(`vcoeff') if `touse', knot(`knots') km(`km')
		local par_x `_kn'
	}
	else {
		if "`vcoeff_par'"=="" {
		local par_x `vcoeff'
		}
		else {
		local par_x `vcoeff_par'
		}
	}

    ** For future applications that will incorporate Sampling Weights
	local sweight  1
	if "`weight'"!="" {
	local sweight  `weight'
	}

    ** Uses Kernel Saved in memory
	if "$kernel_"!="" & "`kernel'"=="" {
		local kernel "$kernel_"
	}
	if  "$kernel_"=="" & "`kernel'"=="" {
		local kernel gaussian
	}
	** as well as BW
	if "$opbw_"!="" & `bw'==-1 {
		local bw $opbw_
	}

	numlist "`degree'", int range(>=0 <4)
	numlist "`wbsrep'", int range(>=1)
	if `degree'==0      local xdepvar `varx' c.`full_x'
	else if `degree'==1 local xdepvar `varx' c.(`varx')#c.`full_x' c.`full_x'
	else if `degree'==2 local xdepvar `varx' c.(`varx')#(c.`full_x'##c.`full_x') c.`full_x'##c.`full_x'
	else if `degree'==3 local xdepvar `varx' c.(`varx')#(c.`full_x'##c.`full_x'##c.`full_x') c.`full_x'##c.`full_x'##c.`full_x'

	******** Basic test: Against 3 types of polynomials.
	tempvar res_p res_vc

	******** No interaction

	qui:reg `y' `xdepvar'    if `touse'==1

	qui:predict double `res_p',res

	******** VC_coeff
	qui:vc_predict `y' `varx' if `touse'==1, vcoeff(`full_x') vcoeff_par(`par_x') bw(`bw')  kernel(`kernel') res(`res_vc') nstat
	******** Estimation Test
	tempname M t
	qui:matrix accum `M'    = `res_p' `res_vc'
	matrix `t'=vecdiag(`M')

	*** Here we will do the "Wild bootstrap"

	*local bsmax=(40*(`wbsrep')-1)
	tempvar y_p bwres_p  bwres_vc
	tempname mata_t aux_t
	mata:`mata_t'= st_matrix("`t'")
	display "Estimating J statistic CI using `wbsrep' Reps"
	display "----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5"
	forvalues i=1/`wbsrep' {
		*small mata program to run regressions? or just here
		local cnt=`cnt'+1
		if `cnt'<50 display _c "."
		else {
		local cnt=0
		display _c "." _n
		}

		capture drop `y_p' `bwres_p' `bwres_vc'
		qui:egen double `y_p'=wberr(`y'), err(`res_p')
		qui:reg `y_p' `varx' `xdepvar' if `touse'==1
		qui:predict double `bwres_p',res
		******** VC_coeff
		qui:vc_predict `y_p' `varx' if `touse'==1, vcoeff(`full_x') vcoeff_par(`par_x') bw(`bw') kernel(`kernel') res(`bwres_vc') nstat
		qui:matrix accum `M'    = `bwres_p' `bwres_vc'
		matrix `M'=vecdiag(`M')
		mata:`mata_t'= `mata_t'\st_matrix("`M'")
	}
	local p975=ceil(`wbsrep'*.975)
	local p950=ceil(`wbsrep'*.950)
	local p900=ceil(`wbsrep'*.900)

	mata:`mata_t'=(`mata_t'[,1]-`mata_t'[,2]):/`mata_t'[,2]
 	mata:`aux_t'=sort(`mata_t'[2..rows(`mata_t'),],1)
	mata:`aux_t'=`mata_t'[1], `aux_t'[`p900'],`aux_t'[`p950'],`aux_t'[`p975']
	mata:st_matrix("`mata_t'",`aux_t')
 	*** To display CV. For now lets just GET CV and save it in the other site
	ereturn clear
	mata:mata drop `mata_t' `aux_t'
	display  _n as text "Specification test."
	if      `degree'==0 display as text "H0: y=x*b0+g*z+e"
	else if `degree'==1 display as text  "H0: y=x*b0+g*z+(z*x)*b1+e"
	else if `degree'==2 display as text  "H0: y=x*b0+g*z+(z*x)*b1+(z^2*x)*b2+e"
	else if `degree'==3 display as text  "H0: y=x*b0+g*z+(z*x)*b1+(z^2*x)*b2+(z^3*x)*b3+e"
	display as text  "H1: y=x*b(z)+e"
	display as text  "J-Statistic      :" as result %6.5f `mata_t'[1,1]
	display as text  "Critical Values"
	display as text  "90th   Percentile:" as result %6.5f `mata_t'[1,2]
	display as text  "95th   Percentile:" as result %6.5f `mata_t'[1,3]
	display as text  "97.5th Percentile:" as result %6.5f `mata_t'[1,4]
end

/*
mata:
void cvreg( string scalar nloerr,string scalar lev,string scalar xbhat)
{
	external x,y,w
	n    = rows(y)
	ixwx=invsym(quadcross(x,w,x))
	xwy=quadcross(x,w,y)
	wx=sqrt(w):*x
	tr=rowsum(wx*ixwx:*wx)
	tr2=(tr:*-1):+1
	yhat=x*ixwx*xwy
	loerr=(y-yhat):/tr2
	st_addvar("double",nloerr)
	st_store((1,n),nloerr,loerr)
	st_addvar("double",lev)
	st_store((1,n),lev,tr)
	st_addvar("double",xbhat)
	st_store((1,n),xbhat,yhat)
}

void cvreg2(string scalar nloerr,string scalar nlev,string scalar nyhat, string scalar nobs)
{
	external y,x,fvc,bw, pvc, vc
	external krn
	// which kernel?

	if (krn=="gaussian") k=1
	else if (krn=="epan") k=2
	else if (krn=="epan2") k=3
	else if (krn=="biweight") k=4
	else if (krn=="cosine") k=5
	else if (krn=="cosine2") k=6
	else if (krn=="parzen") k=7
	else if (krn=="rectan") k=8
	else if (krn=="trian") k=9
	else if (krn=="logistic") k=10
	else if (krn=="tricube") k=11
	else if (krn=="triweight") k=12

	n    = rows(y)
	yhat2=J(rows(y),1,0)
	loerr2=J(rows(y),1,0)
	lev2=J(rows(y),1,0)
	nnobs=J(rows(y),1,0)
 	for(i=1;i<=rows(vc);i++) {
 		z=(fvc:-vc[i]):/bw
	    w=kweight(z,k)
		ixwx=invsym(quadcross(x,w,x))
		xwy=quadcross(x,w,y)
		wx=sqrt(w):*x
		lrv=quadrowsum(wx*ixwx:*wx)
		yhat=x*ixwx*xwy
		loerr=(y-yhat):/((lrv:*-1):+1)
		loerr=loerr:*(loerr:!=.)
		lev2 =lev2+lrv:*(pvc:==vc[i])
		yhat2 =yhat2+yhat:*(pvc:==vc[i])
		loerr2=loerr2+loerr:*(pvc:==vc[i])
		nnobs=nnobs+(J(rows(y),1,0):+quadsum(w)):*(pvc:==vc[i])
	}
	stata("sum z1")
 	st_addvar("double",nloerr)
	st_store((1,n),nloerr,loerr2)
	st_addvar("double",nlev)
	st_store((1,n),nlev,lev2)
	st_addvar("double",nyhat)
	st_store((1,n),nyhat,yhat2)
	st_addvar("double",nobs)
	st_store((1,n),nobs,nnobs)
	stata("sum z1")
}


real matrix kweight(real matrix z,real scalar k)
{
	if (k==1) {
		kz=normalden(z):/normalden(0)
	}
	else if (k==2) {
		kz=(-0.2*(z:^2):+1):*(abs(z):<(5^.5))
	}
	else if (k==3) {
		kz= (-z:^2:+1):*(abs(z):<1)
	}
	else if (k==4) {
		kz=(-z:^2:+1):^2:*(abs(z):<1)
	}
	else if (k==5) {
		kz=(cos(2*pi()*z):+1):/2:*(abs(z):<0.5)
  	}
	else if (k==6) {
		kz=cos(pi()/2*z):*(abs(z):<1)
	}
	else if (k==7) {
		kz=((-6*z:^2+6*abs(z):^3):+1):*(abs(z):<=0.5)+(-abs(z):+1):^3:*2:*((abs(z):>0.5):*(abs(z):<=1))
	}
	else if (k==8) {
		kz=abs(z):<=1
	}
	else if (k==9) {
		kz=(-abs(z):+1):*(abs(z):<=1)
	}
	else if (k==10) {
		kz=(exp(z)+exp(z:*-1):+2):^-1:*4
	}
	else if (k==11) {
		kz=(-abs(z):^3:+1):^3:*(abs(z):<1)
	}
	else if (k==12) {
		kz=(-abs(z):^2:+1):^3:*(abs(z):<1)
	}
	return(kz)
}

end

 */