program define npreg
*! v 0.1 gmh 10 feb 01
version 6
di in red "BETA VERSION"
syntax varname [if] [in], Nonpar(varlist min=1 max=2) Param(varlist) [ BW1(real 4) BW2(real 4) ]
tempfile temp
qui save `temp'
tempvar touse
mark `touse' `if' `in'
markout `touse' `nonpar' `param'
qui keep if `touse'

* unabbreviate the variables
unab nonpar:`nonpar'
unab param: `param'
unab y: `varlist'

tokenize `nonpar'
local z1 "`1'"
local z2 "`2'"

local numx : word count `param'
tokenize `"`param'"'
local i=1
while `i'<=`numx' {
	local x`i' "``i''"
	local i = `i'+1
	}
	

* generate count matrices
qui count
scalar N = r(N)
qui sum `z1'
scalar stdz1=r(sd)
scalar h1 = `bw1' * scalar(stdz1)  * scalar(N) ^ (-.2)
qui sum `z2'
scalar stdz2=r(sd)
if scalar(stdz2)==0 {
	scalar stdz2=1
	di in red "Scalar(stdz2) recoded"
	}
scalar h2 = `bw2' * scalar(stdz2)  * scalar(N) ^ (-.2)
qui ta `z1', matcell(cz1)
qui ta `z2', matcell(cz2)


** z1
tempvar g1 g2
egen `g1' = group(`z1')
qui sum `g1' 
scalar nz1 = r(max)

preserve
collapse (mean) `z1', by(`g1')
mkmat `z1' , mat(z1val)
restore
qui tabulate `g1', matcell(cz1)
** end z1
** z2
egen `g2' = group(`z2')
qui sum `g2' 
scalar nz2 = r(max)

preserve
collapse (mean) `z2', by(`g2')
mkmat `z2' , mat(z2val)
restore
qui tabulate `g2', matcell(cz2)

** end z2



tempvar kern1 kern2 wil y1 y2
tempname kern1 kern2
qui gen `kern1'=.
qui gen `kern2'=.
qui gen `wil'=.
qui gen `y1'=.
qui gen `y2'=.

/*
label var x1_1 "x1 for alpha = 1"
label var x2_1 "x2 for alpha = 1"
label var x1_2 "x1 for alpha = 2"
label var x2_2 "x2 for alpha = 2"
*/

local i = 1
while `i'<=`numx' {
	tempvar x`i'_1 x`i'_2
	qui gen `x`i'_1'=.
	qui gen `x`i'_2'=.
	local i = `i'+1
	}

* z1
di _n " "
di in gr "Matrix loop for variable " in ye "`z1'"
			di _col(1) "0%" _con _col(18) "25%" _con _col(39) "50%" /*
 			*/ _con _col(59) "75%" _con _col(75) "100%" _new
			local repcol=scalar(N)/77
			local lastcol=0

local i = 1
while `i' <= scalar(N) {
	local l = 1
	while `l' <=scalar(N) {
		scalar `kern1' = exp(-.5 * ((`z1'[`l'] - `z1'[`i'])/scalar(h1) ) ^2)
		 qui replace `kern2' = scalar(`kern1')*exp(-.5 * ((`z2'[`l'] - `z2')/scalar(h2) ) ^2)
		qui replace `kern2'=0 in `l'
		sum `kern2', mean
		qui replace `wil'=r(sum) in `l'
		local l = `l'+1
		}
	* rescale to sum to 1
	sum `wil', mean
	qui replace `wil'=`wil'/r(sum)
	sum `y' [w=`wil'], meanonly
	qui replace `y1' = r(sum) in `i'

	local n=1
	while `n'<=`numx' {
	sum `x`n'' [w=`wil'], meanonly
	qui replace `x`n'_1' = r(sum) in `i'
	local n= `n'+1
	}



				while `i'/`repcol'>`lastcol' {
					di _con in green "*" 
					local lastcol=`lastcol'+1
				}

	local i = `i'+1
}

* z2

di _n " "
di in gr "Matrix loop for variable " in ye "`z2'"
			di _col(1) "0%" _con _col(18) "25%" _con _col(39) "50%" /*
			*/ _con _col(59) "75%" _con _col(75) "100%" _new
			local lastcol=0


qui replace `wil'=.
local i = 1
while `i' <= scalar(N) {
	local l = 1
	while `l' <=scalar(N) {
		scalar `kern2' = exp(-.5 * ((`z2'[`l'] - `z2'[`i'])/scalar(h2) ) ^2)
		qui replace `kern1' = scalar(`kern2')*exp(-.5 * ((`z1'[`l'] - `z1')/scalar(h1) ) ^2)
		qui replace `kern1'=0 in `l'
		sum `kern1', mean
		qui replace `wil'=r(sum) in `l'
		local l = `l'+1
		}
	* rescale to sum to 1
	qui sum `wil', mean
	qui replace `wil'=`wil'/r(sum)
	sum `y' [w=`wil'], meanonly
	qui replace `y2' = r(sum) in `i'


	local n=1
	while `n'<=`numx' {
	sum `x`n'' [w=`wil'], meanonly
	qui replace `x`n'_2' = r(sum) in `i'
	local n= `n'+1
	}


				while `i' / `repcol' >`lastcol' {
					di _con in green "*" 
					local lastcol=`lastcol'+1
				}

	local i = `i'+1
}

** now create new variables: step 1
tempvar ytil 

qui gen `ytil'=`y' - `y1' - `y2'
local n=1
while `n'<=`numx' {
	tempvar x`n'til
	qui gen `x`n'til'=`x`n'' - `x`n'_1' - `x`n'_2'
	local rhslist "`rhslist' `x`n'til' "
	local n=`n'+1
	}

qui reg `ytil' `rhslist'
mat b1=e(b)

*** make the PHI matrix
mat phi = J(`numx'+1,`numx'+1,0)
/*
local i =1
tempname tempphi
while `i'<=scalar(N) {
	qui mat accum `tempphi' = `rhslist' in `i'
	mat phi=phi+`tempphi'
	local i = `i'+1
	}
*/
qui mat accum phi=`rhslist'
mat phi = phi * 1/scalar(N)
mat rowname phi = `param' _cons
mat colname phi = `param' _cons


* do second step
mat bstar = b1[1,1..`numx']
mat colname bstar =`param'
tempvar xb
matrix score `xb' = bstar
sum `y', meanonly
scalar b0 = r(mean)
sum `xb', meanonly
scalar b0=scalar(b0) - r(mean)

qui replace `xb' = `xb'-scalar(b0)
tempvar ystage2 wt zstage2
qui gen `ystage2'=`y'-`xb'

qui gen `wt'=.
qui gen `zstage2'=.
* for z1
mat a1 = J(scalar(nz1), scalar(nz2),0)
local i = 1
while `i' <= scalar(nz1) {
	local j = 1
	while `j'<=scalar(nz2) {
		qui replace `wt' = exp(-.5 * ((`z1' - z1val[`i',1])/scalar(h1))^2) * exp(-.5 * ((`z2' - z2val[`j',1])/scalar(h2))^2) 
		qui replace `zstage2'=`z1'-z1val[`i',1]

		qui reg `ystage2' `zstage2' [w=`wt']
		mat a1[`i',`j'] = _b[_cons]
		local j = `j'+1
		}
local i = `i'+1
}

* for z2
mat a2 = J(scalar(nz2), scalar(nz1),0)
local i = 1
while `i' <= scalar(nz2) {
	local j = 1
	while `j'<=scalar(nz1) {
		qui replace `wt' = exp(-.5 * ((`z2' - z2val[`i',1])/scalar(h2))^2) * exp(-.5 * ((`z1' - z1val[`j',1])/scalar(h1))^2) 
		qui replace `zstage2'=`z2'-z2val[`i',1]

		qui reg `ystage2' `zstage2' [w=`wt']
		mat a2[`i',`j'] = _b[_cons]
		local j = `j'+1
		}
local i = `i'+1
}



* now we have the alpha matrices
* make the m matrices
matrix m1 = a1 * cz2 / scalar(N)
matrix m2 = a2 * cz1 / scalar(N)

* now make the g's
matrix g1 = m1 - J(scalar(nz1),1,1) *m1'*cz1  /scalar(N)
matrix g2 = m2 - J(scalar(nz2),1,1) * m2'*cz2  /scalar(N)

matrix finalz1=g1,z1val
matrix finalz2=g2,z2val

mat colname finalz1 = ghat1 z1val
mat colname finalz2=ghat2 z2val

svmat finalz1, n(col)
svmat finalz2,n(col)

gr ghat1 z1val, c(s) saving(g1,replace) t1("Predictions for `z1': bandwidth = `bw1'")
gr ghat2 z2val, c(s) saving(g2,replace) t1("Predictions for `z2': bandwidth = `bw2'")
gr using g1 g2

* still need to fix intercept

qui reg `ytil' `rhslist'
tempvar pred
qui predict `pred'
qui replace `pred'=`pred'-_b[_cons]+scalar(b0)
tempvar u
local dim1=rowsof(matrix(z1val))
local dim2=rowsof(matrix(z2val))
local mindim=min(`dim1',`dim2')

local i =1
tempvar gz1 gz2
qui gen gz1=0
qui gen gz2=0
while `i'<= `dim1' {
	qui replace gz1=g1[`i',1] if `g1'==`i'
	local i =`i'+1
	}
local i=1
while `i'<= `dim2' {
	qui replace gz2=g2[`i',1] if `g2'==`i'
	local i =`i'+1
	}
		

qui gen uhat=`y' -`pred'- gz1 - gz2


mat b = e(b)
mat V = e(V)
mat coln b = `param'  _cons
mat coln V = `param' _cons
mat rown V = `param' _cons
mat b[1,`numx'+1]=scalar(b0)
estimates post b V, depname(`y')

di _n " "
di "Fan-Li-Li's Non-parametric regression"
di " "
estimates display

local i=1

di " "
di in gr       "Estimated values for z1"  _col(30)    "Estimated values for z2"
di " "
di in gr      " z1         pred val " _col(30)    "  z2         pred val "
di in gr      "----------------------------------------------------------"
	while `i'<=`mindim' {
di in ye _col(2) %4.2f z1val[`i',1] _col(7) in gr "|" _col(9) in ye %12.8f g1[`i',1] _col(32) %4.2f  z2val[`i',1] _col(37) in gr "|" _col(39) in ye %12.8f g2[`i',1]
	local i = `i'+1
	}
while `i'<=`dim1' {
di in ye _col(2)  %4.2f   z1val[`i',1] _col(7) in gr "|" _col(9) in ye %12.8f g1[`i',1] 
	local i=`i'+1
	}	
	while `i'<=`dim2' {
di in ye                 _col(32)  %4.2f  z2val[`i',1] _col(37) in gr "|" _col(39) in ye %12.8f g2[`i',1]
	local i = `i'+1
	}

use `temp', clear
end