*! version 1.0.2    0724       C F Baum
program define hegy12, rclass
	version 6
* lags numlist rather than integer 0
	syntax varname(ts) [if] [in] [ , Lags(numlist >0) Det(string) ]

   	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
	}
	qui tsset
	if r(unit1) != "m" {
		di in red "hegy12 must be used with monthly data"
		exit
		}
	local tmin=r(tmins)
	local tmax=r(tmaxs)
	tempvar count trd patn seas1 x y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 y11 y12 y13 et
	tempname b v cv r32 r3
* seas (seasonal dummies and constant) default
	local aux "err"
	if "`det'"=="" | "`det'" =="seas" {
		local aux "L(1/11).`seas1' "
		local daux "Seasonal dummies + constant"
		local case 3
		}
	if "`det'"=="none" {
		local aux ", noc"
		local daux "None"
		local case 1
		}
	if "`det'"=="const" {
		local aux " "
		local daux "Constant"
		local case 2
		}
	if "`det'"=="strend" {
		local aux "L(1/11).`seas1' `trd'"
		local daux "Seasonal dummies + constant + trend"
		local case 4
		}
	if "`aux'"=="err" {
		di in red "Error: det must be none, const, seas, or strend"
		exit
		}
	qui gen `count'=sum(`touse')
	local nobs=`count'[_N]
	local aug ""
	if "`lags'" != "" {
	tokenize `lags'
	local i 1
        while "``i''" != "" {
		local aug = "`aug'  L``i''.`y5' "
		local i = `i' +1
		}
	}
	drop `count'
	qui gen `count'=sum(`touse')
	local nobs=`count'[_N]
* nobs reflects the max lag provided or calculated
	gen double `trd' = _n
	egen `patn' = fill(1/12 1/12)
	qui gen `seas1'= (`patn'==1)
* generate y1..y13 per Beaulieu/Miron p.308, backshifted one
	scalar `r32'=sqrt(3.0)/2.0
	scalar `r3' =sqrt(3.0)
	qui { gen double `x' = `varlist' if `touse'
	gen double `y1' =  L.`x' + L2.`x' + L3.`x' + L4.`x' + L5.`x' + L6.`x' + /*
	*/ L7.`x' + L8.`x' + L9.`x' + L10.`x' + L11.`x' + L12.`x'
	gen double `y2' = -L.`x' + L2.`x' - L3.`x' + L4.`x' - L5.`x' + L6.`x' - /*
	*/ L7.`x' + L8.`x' - L9.`x' + L10.`x' - L11.`x' + L12.`x'
	gen double `y3' = -L2.`x' + L4.`x' - L6.`x' + L8.`x' - L10.`x' + L12.`x'
	gen double `y4' =  -L.`x' + L3.`x' - L5.`x' + L7.`x' -  L9.`x' + L11.`x'
	gen double `y5' = -0.5 *(L.`x' + L2.`x' -2.0*L3.`x' + L4.`x' + L5.`x' -2.0*L6.`x' + /*
	*/ L7.`x' + L8.`x' - 2.0*L9.`x' + L10.`x' + L11.`x' -2.0*L12.`x' )
	gen double `y6' = `r32'*(L.`x' - L2.`x' + L4.`x' - L5.`x' + L7.`x' - L8.`x' + L10.`x' - L11.`x')
	gen double `y7' = 0.5 *(L.`x' - L2.`x' -2.0*L3.`x' - L4.`x' + L5.`x' + 2.0*L6.`x' + /*
	*/ L7.`x' - L8.`x' - 2.0*L9.`x' - L10.`x' + L11.`x' + 2.0*L12.`x' )
	gen double `y8' = -`r32'*(L.`x' + L2.`x' - L4.`x' - L5.`x' + L7.`x' + L8.`x' - L10.`x' - L11.`x')
	gen double `y9' = -0.5 *(`r3'*L.`x' - L2.`x' + L4.`x' -`r3'*L5.`x' + 2.0*L6.`x'  /*
	*/ -`r3'*L7.`x' + L8.`x' - L10.`x' + `r3'*L11.`x' -2.0*L12.`x' )
	gen double `y10' = 0.5 *(L.`x' -`r3'*L2.`x' +2.0*L3.`x' -`r3'*L4.`x' + L5.`x'  /*
	*/ - L7.`x' + `r3'*L8.`x' -2.0*L9.`x' + `r3'*L10.`x' - L11.`x' )
	gen double `y11' = 0.5 *(`r3'*L.`x' + L2.`x' - L4.`x' -`r3'*L5.`x' - 2.0*L6.`x'  /*
	*/ -`r3'*L7.`x' - L8.`x' + L10.`x' + `r3'*L11.`x' + 2.0*L12.`x' )
	gen double `y12' = -0.5 *(L.`x' + `r3'*L2.`x' +2.0*L3.`x' + `r3'*L4.`x' + L5.`x'  /*
	*/ - L7.`x' - `r3'*L8.`x' -2.0*L9.`x' - `r3'*L10.`x' - L11.`x' )
	gen double `y13' = `x' - L12.`x'
	}
* run the basic regression to generate pi(1)..pi(12)
 	qui {
	reg `y13' `y1' `y2' `y3' `y4' `y5' `y6' `y7' `y8' `y9' `y10' `y11' `y12' `aug' `aux'
	predict double `et',r
	mat `b'=e(b)
	mat `v'=e(V)
	local n=e(N)
	local j 1
	while `j'<13 {
		 return scalar pi`j' =`b'[1,`j']/sqrt(`v'[`j',`j'])
		 local j = `j'+1
		 }

* joint (F) test on pi(3), pi(4)
	test `y3' `y4'
	return scalar F34 = `r(F)'
* joint (F) test on pi(2) - pi(4)
	test `y2' `y3' `y4'
	return scalar F234 = `r(F)'
* joint (F) test on pi(1) - pi(4)
	test `y1' `y2' `y3' `y4'
	return scalar F1234 = `r(F)'
	GetCritH `case'
	mat `cv'= r(cv)
	return scalar N=`n'
*	local n4=`n'+4+`lags'
        local n4=`n'
  	}
* for now, skip run the auxiliary regression (only for diagnostics)
*	reg `et' L(1/4).`et' `y1' `y2' `y3' `y4' `aug' `aux'
*

	di in gr " "
	di in gr "HEGY Monthly seasonal unit root test for `varlist'"
*	summ `varlist' if e(sample)
*	summ `y5' `y1' `y2' `y3' `y4' `aug' `aux' if e(sample)
	di in gr " "
	di in gr "Number of observations  : `n4' " /* ( based on `tmin' - `tmax')*/
	di in gr "Deterministic variables : `daux'"
	if "`aug'" !="" { di in gr "Augmented by lags : `lags'"}
	di in gr " "
* could plug in 1% criticals here
	di in gr _col(11) "Stat" _col(20) "5% critical" _col(34) "10% critical"
	di in gr _dup(78) "-"
	di "t[Pi1]" _col(8) %9.3f `return(pi1)' _col(19) %9.3f `cv'[1,3] /*
	*/ _col(33) %9.3f `cv'[1,4]
	di "t[Pi2]" _col(8) %9.3f `return(pi2)' _col(19) %9.3f `cv'[1,7] /*
	*/ _col(33) %9.3f `cv'[1,8]
	di "t[Pi3]" _col(8) %9.3f `return(pi3)' _col(19) %9.3f `cv'[1,11] /*
	*/ _col(33) %9.3f `cv'[1,12]
	di "t[Pi4]" _col(8) %9.3f `return(pi4)' _col(19) %9.3f `cv'[1,15] /*
	*/ _col(33) %9.3f `cv'[1,16]
	di "t[Pi5]" _col(8) %9.3f `return(pi5)' _col(19) %9.3f `cv'[1,11] /*
	*/ _col(33) %9.3f `cv'[1,12]
	di "t[Pi6]" _col(8) %9.3f `return(pi6)' _col(19) %9.3f `cv'[1,15] /*
	*/ _col(33) %9.3f `cv'[1,16]
	di "t[Pi7]" _col(8) %9.3f `return(pi7)' _col(19) %9.3f `cv'[1,11] /*
	*/ _col(33) %9.3f `cv'[1,12]
	di "t[Pi8]" _col(8) %9.3f `return(pi8)' _col(19) %9.3f `cv'[1,15] /*
	*/ _col(33) %9.3f `cv'[1,16]
	di "t[Pi9]" _col(8) %9.3f `return(pi9)' _col(19) %9.3f `cv'[1,11] /*
	*/ _col(33) %9.3f `cv'[1,12]
	di "t[Pi10]" _col(8) %9.3f `return(pi10)' _col(19) %9.3f `cv'[1,15] /*
	*/ _col(33) %9.3f `cv'[1,16]
	di "t[Pi11]" _col(8) %9.3f `return(pi11)' _col(19) %9.3f `cv'[1,11] /*
	*/ _col(33) %9.3f `cv'[1,12]
	di "t[Pi12]" _col(8) %9.3f `return(pi12)' _col(19) %9.3f `cv'[1,15] /*
	*/ _col(33) %9.3f `cv'[1,16]
	di " "
* figure out how to deal with Fs!
/* 	di "F[3-4]" _col(8) %9.3f `return(F34)' _col(19) %9.3f `cv'[1,9] /*
	*/ _col(33) %9.3f `cv'[1,10]
	di "F[2-4]" _col(8) %9.3f `return(F234)' _col(19) %9.3f
	di "F[1-4]" _col(8) %9.3f `return(F1234)' _col(19) %9.3f
*/
end

program define GetCritH, rclass
* values from Beaulieu/Miron Table A.1 (blocks correspond to cases 1,2,3,-,4) for 240 obs.
* for fractiles 0.01, 0.025, 0.05, 0.10 (for t), .9, 0.95, 0.975 (for F)
        args c
	tempname zt cv
	mat `zt' = ( /*
        */         -2.51, -2.18,  -1.89,  -1.58,   -2.53, -2.16,  -1.87,  -1.57,   -2.50, -2.16,  -1.88,  -1.55,   -2.31, -1.95,  -1.63,  -1.27,   2.34,  3.03,  3.71,    4.6 \ /*
        */         -3.35, -3.06,  -2.80,  -2.51,   -2.48, -2.15,  -1.89,  -1.57,   -2.51, -2.16,  -1.87,  -1.54,   -2.30, -1.93,  -1.62,  -1.27,   2.32,  3.01,  3.68,    4.6 \ /*
        */         -3.32, -3.02,  -2.76,  -2.47,   -3.28, -3.01,  -2.76,  -2.48,   -3.83, -3.51,  -3.25,  -2.95,   -2.61, -2.21,  -1.85,  -1.45,   5.27,  6.26,  7.19,    8.35 \ /*
        */         -3.83, -3.54,  -3.28,  -2.99,   -3.31, -3.02,  -2.75,  -2.47,   -3.79, -3.50,  -3.24,  -2.95,   -2.57, -2.18,  -1.85,  -1.45,   5.25,  6.23,  7.14,    8.33 )

        matrix `cv'=`zt'[`c',1...]
	return matrix cv  `cv'
end
exit