*! version 1.0.1     0411     C F Baum  (STB-57: sts15)
*  from dfgls 1.2.0, 0405 and J.Barkoulas RATS procedure kpss.src
*  mod to use gen, predict double
program define kpss, rclass
	version 6

	syntax varname(ts) [if] [in] [ , Maxlag(integer -1) noTrend ]  

   	marksample touse
			/* get time variables */
	_ts timevar, sort
	markout `touse' `timevar'
	tsreport if `touse', report
	if r(N_gaps) {
		di in red "sample may not contain gaps"
		exit
	}
	tempvar trd y dty psum psum2 
	tempname nobs stat count kmax s20 sumst2 kpss0 A p1 p5 p10 yvar
		
	if "`trend'"=="notrend" {
		local trend 0
		local stat "level stationary"
		local p1 0.739
		local p2h 0.574
		local p5 0.463
		local p10 0.347
		}
	else {
		local trend 1
		local stat "trend stationary"
		local p1 0.216
		local p2h 0.176
		local p5 0.146
		local p10 0.119
		}
	gen `count'=sum(`touse')
	local nobs=`count'[_N]
	if `maxlag'==-1 {
* set maxlag via Schwert criterion (Ng/Perron JASA 1995)
		local maxlag = int(12*(`nobs'/100)^0.25)
		local kmax = "Maxlag = `maxlag' chosen by Schwert criterion" 
		}
	else {
		local kmax "Maxlag = `maxlag'"
		}
	markout `touse' L(1/`maxlag').`touse'
	drop `count'
	gen `count'=sum(`touse')
	local nobs=`count'[_N]
* nobs reflects the max lag provided or calculated 
	if `trend' {gen `trd' = _n}
	qui {
		gen double `y' = `varlist' 
		local yvar  `varlist'
* run the OLS regression to detrend (or demean) the data, calc resids
		if `trend' {reg  `y'  `trd' if `touse'}
		else {reg `y' if `touse'}
		predict double `dty' if `touse',r
		local s20=e(df_r)/e(N)*e(rmse)^2
		local nobs=e(N)
* calculate partial sum series and its squares
		gen double `psum' = sum(`dty') if `touse'
		gen double `psum2' = `psum'^2  if `touse'
		summ `psum2',meanonly
* numerator: average value of squared partial sum series / T
		local sumst2 = r(mean)/`nobs'
		local kpss0=`sumst2'/`s20'		
		if `maxlag'>0 {
			local l 1
* accumulate cross products of autocovariances
			mat accum A = `dty' L(1/`maxlag').`dty',noc
			while `l'<=`maxlag' {
				local s2`l'=0
				local s 1
				while `s'<=`l' {
					local w=1-(`s'/(`l'+1))
* denominator: accumulate long run variance using Newey-West weighting
					local s2`l'=`s2`l''+`w'*A[`s'+1,1]
					local s=`s'+1
					}
* s2l values differ marginally from RATS since RATS uses all available
* observations to compute mcov
				local s2`l'=2.0*`s2`l''/`nobs'+`s20'
				local kpss`l'=`sumst2'/`s2`l''
				local l=`l'+1
				}
			}
		}
	di " "
	di in gr "KPSS test for `yvar'"
	di " "
	di in gr "`kmax'"
	di " "
	di in gr "Critical values for H0: `yvar' is `stat'"
	di " "
	di in gr "10%: `p10'  5% : `p5'  2.5%: `p2h'  1% : `p1'"
	di " "
	di in gr "Lag order    Test statistic"
	di       "    0       " %8.3g `kpss0'
	return scalar N = `nobs'
	return scalar kpss0 = `kpss0'
	if `maxlag'>0 {
		local l 1
		while `l'<=`maxlag' {
		di  _col(4) %2.0f `l' _col(13) %8.3g `kpss`l''
		return scalar kpss`l'=`kpss`l''
		local l=`l'+1
			}
		}
	end
	exit