*! version 0.0.2  28nov2007
*!   bug fix in the standard error calculation for xthp_nox
*!   'if' is handled correctly
*! version 0.0.1  16feb2007

program define xthp, eclass byable(recall) sortpreserve
	version 6
	if replay() {
		if `"`e(cmd)'"' != `"xthp"' { error 301 }
		if _by() { error 190 }
		error 301
		/*NOTREACHED*/
	}
	syntax varlist(ts) [if] [in] [, TRend FE RE BE ITER(integer 100) EPSilon(real 1.0e-6) * ]

	if ("`fe'"!="")+("`re'"!="")+("`be'"!="")>1 {
		di as err "choose only one of fe, re or be"
		exit 198
	}

	capture tsset

	local tmin = r(tmin)
	local tmax = r(tmax)

	capture local id "`r(panelvar)'"
	if "`id'"=="" {
		di as err "must tsset data and specify panelvar"
		exit 459
	}
	capture local t "`r(timevar)'"
	if "`t'"=="" {
		di as err "must tsset data and specify timevar"
		exit 459
	}

	local tvar `t'

	tempvar touse touseb
	mark `touse' `if' `in'
	markout `touse' `id'

	if _by() {
		qui replace `touse' = 0 if `_byindex' != _byindex()
	}
	qui sum `touse'
	if r(max) == 0 {
		di as err "no observations"
		exit 2000
	}
	qui gen byte `touseb' = `touse'

	preserve
	if ("`if'"!="") {
          quietly keep `if'
	}

	tempvar orig
	qui gen byte `orig' = 1
	/* fill in missing times with missing observations */
	tsfill, full

	sort `id' `t'
	qui replace `orig' = 0 if `orig' >= .
	qui count if `t'==.
	if r(N)>0 {
		di as err "no missing values in time variable (`t') allowed"
		exit 459
	}

	local varlist "`varlist' "
	tokenize `varlist'
	local depvar `1'
	local depvar: subinstr local depvar "." ".", count(local legal)
	if `legal' > 0 {
		di as err "`depvar' invalid in this context"
		di as err /*
		*/"dependent variable may not contain time-series operators"
		exit 198
	}
	macro shift
	local indeps "`*'"
	local exog `indeps'

	tempvar uniqi uniqt
	tempname inum tnum
	/* Count the number of individual units */
	qui by `id': gen byte `uniqi' = (`touseb' & _n==_N)
	qui summ `uniqi'
	scalar `inum' = _result(18)

	/* Count the number of temporal units */
	sort `t'
	qui by `t': gen byte `uniqt' = (`touseb' & _n==_N)
	qui summ `uniqt'
	scalar `tnum' = _result(18)

	sort `id' `t'

	tempvar dy dy1
	tempname s2m se wm t NT

	qui gen `dy' = D.`depvar'
	qui gen `dy1' = L.`dy' `if'

	local opt = "`fe'`re'`be'"

	if "`exog'" == "" { /* no X variables */
		if "`fe'"!="" | "`re'"!="" | "`be'"!="" {
			di as err "Option `opt' is irrelevant. Ignored."
		}
		if "`trend'" == "" { /* no trend */
			xthp_nox `depvar' if `touseb', i(`id')
		}
		else { /* with trend */
			/* This is to be implemented later. */
			di as err "Trend is not supported yet. Ignored."
			xthp_nox `depvar' if `touseb', i(`id')
		}
	}
	else { /* with X variables */
		if "`trend'" == "" { /* no trend */
			xthp_x `varlist' if `touseb', i(`id') `fe' `re' `be' iter(`iter') epsilon(`epsilon')
		}
	}

	est local ivar "`id'"
	est local tvar "`tvar'"
	est local depvar "`depvar'"

	qui summ `dy1', meanonly
	scalar `NT' = r(N)
	est scalar N = `NT'
	est scalar N_g = `inum'
	est scalar N_t = `tnum'

	restore

	Disp

end

program define xthp_nox, eclass
	syntax varname [if] [, I(varname) * ]

	tempvar dy x z resid xe xe1
	tempname B V xz bhat wm sxz NT

	local depvar `1'
	local id `i'

	qui gen `dy' = D.`depvar'
	qui gen `x' = L.`dy'
	qui gen `z' = 2*`dy' + `x'

	qui est clear

	qui matrix accum `sxz' = `x' `z' `if', noc
	matrix `B' = `sxz'[2,1]/`sxz'[1,1]
	scalar `bhat' = `B'[1,1]

	qui gen `resid' = `z' - `bhat' * `x' `if'
	qui gen `xe' = `x' * `resid'
	qui by `id': egen `xe1' = sum(`xe')
	qui by `id': replace `xe1' = `xe1'/sqrt(_N)
	qui matrix accum `wm' = `xe1' `if', noc
	matrix `V' = `wm'[1,1] / (`sxz'[1,1]^2)
	qui summ `x' `if', meanonly
	scalar `NT' = r(N)

	mat colnames `B' = "L.`depvar'"
	mat rownames `B' = "L.`depvar'"
	mat colnames `V' = "L.`depvar'"
	mat rownames `V' = "L.`depvar'"

	est post `B' `V', depname("`depvar'")
end

program define xthp_x, eclass
	syntax varlist [if] [, I(varname) FE RE BE ITER(integer 100) EPSilon(real 1.0e-6) * ]

	tempname b_y b_x v_x v_y bhat0 bhat1
	local id `i'
	scalar `bhat0' = .
	scalar `bhat1' = 0

    local varlist "`varlist' "
    tokenize `varlist'
    local depvar `1'
    macro shift
    local indeps "`*'"

	local nindeps = 0
	foreach x of local indeps {
		local nindeps = `nindeps' + 1
	}

	local j = 1
	while `j' <= `iter' & abs(`bhat1'-`bhat0')>`epsilon' {
		tempvar y res
		scalar `bhat0' = `bhat1'
		qui gen `y' = `depvar' - `bhat0' * L.`depvar'
		local xvar
		foreach x of local indeps {
			tempvar xtmp
			qui gen `xtmp' = `x' - `bhat0' * L.`x'
			local xvar "`xvar' `xtmp' "
		}
		if "`re'"=="" & "`be'"=="" {
			qui xtreg `y' `xvar' `if', fe
		}
		else {
			qui xtreg `y' `xvar' `if', `fe'`re'`be'
		}
		matrix `b_x' = e(b)
		matrix `v_x' = e(V)
		qui gen `res' = `depvar'
		local k = 1
		foreach x of local indeps {
			qui replace `res' = `res' - `b_x'[1,`k'] * `x'
			local k = `k' + 1
		}
		drop `xvar'
		qui xthp_nox `res' `if', i(`id')
		matrix `b_y' = e(b)
		matrix `v_y' = e(V)
		scalar `bhat1' = `b_y'[1,1]

		local j = `j' + 1
	}
	tempname b_xy v_xy z
	matrix `b_xy' = ( `b_y'[1,1], `b_x'[1,1..`nindeps'] )
	matrix `z' = J(1,`nindeps',0)
	matrix `v_xy' = `v_y',`z' \ `z'',(`v_x'[1..`nindeps',1..`nindeps'])
	tokenize "L.`depvar' `indeps' "
	mat colnames `b_xy' = `*'
	mat colnames `v_xy' = `*'
	mat rownames `v_xy' = `*'

	est clear
	est post `b_xy' `v_xy', depname("`depvar'")
end

program define Disp, eclass

	#delim ;
		di _n as txt "Han-Phillips dynamic panel-data estimation"
			_col(49) as txt "Number of obs" _col(68) "="
			_col(70) as res %9.0f e(N) ;
		di as txt "Group variable (i): " as res abbrev("`e(ivar)'", 12) as txt
			_col(49) "Number of groups" _col(68) "="
			_col(70) as res %9.0g e(N_g) ;
		di as txt "Time  variable (t): " as res abbrev("`e(tvar)'", 12) as txt
			_col(49) "Number of periods" _col(68) "="
			_col(70) as res %9.0g e(N_t) _n ;

	#delim cr

	est di
end

exit