*! version 1.0.0  16Feb2003  CFBaum
*********************************************************************
*       THIS PROGRAM COMPUTES THE POSITIVE DEFINITE, ROBUST         *
*       STATISTIC FOR FORECAST COMPARISON                           *
*       EQUATIONS (42) OR (43) OF "FORECAST COMPARISON IN L2"       *
*       Rutgers University Working Paper 1995-24                    *
*       Bruce Mizrach                                               *
*********************************************************************
program define mizrachl2, rclass
	version 7
	syntax varlist(ts min=2 max=2) [if] [in] [, k(integer -1) ]
	marksample touse
	_ts timevar, sort
	markout `touse' `timevar'
	tempvar  count u v uv dot
	tempname T suvuv uvuvtot stat1
* get the two forecast variables
	tokenize `varlist'
	local fcvar1 `1'
	local fcvar2 `2'
* set k = order of dependence
	qui	gen `count'=sum(`touse')
	local nobs=`count'[_N]
	if `k'==-1 {
* set k from sample size
		local k = int((`nobs')^0.3333)+1
		local kchoix = "Order of dependence = `k' chosen from sample size"
		}
	else {
		local kchoix "Order of dependence = `k'"
		}
* calculate u, v, uv
	qui {
		gen double `u' = `fcvar1' - `fcvar2' if `touse'
		gen double `v' = `fcvar1' + `fcvar2' if `touse'
		gen double `uv' = `u'*`v'
		summ `uv', meanonly
		local sn1 = r(mean)
		scalar `T' = r(N)
		}
* generate raw moments of uv series over all available obs,
* apply Bartlett kernel weights
	mat `suvuv' = J(`k',1,0)
	qui gen double `dot' = `uv'^2 if `touse'
	qui summ `dot', meanonly
	scalar `uvuvtot' = r(sum)
    forv l = 1/`k' {
        	local w = 1 - (`l'/(`k'+1))
        	qui replace `dot' = `uv'*L`l'.`uv'
        	qui summ `dot', meanonly
			scalar `uvuvtot' = `uvuvtot' + 2*`w'*r(sum)
        	}
	local sd1 = `uvuvtot'/(`T'^2)
	scalar `stat1' = `sn1'/sqrt(`sd1')
	qui {
		tsset
		summ `r(timevar)' if `touse'
		}
	return scalar stat1 = `stat1'
	return scalar p = 2*normprob(-abs(return(stat1)))
	return scalar nobs = `T'
	return scalar k = `k'

	di _n "Mizrach L2 forecast comparison test"
*	di in gr "Forecast period : " %tm r(min) " -" %tm r(max) "  N : `nobs'"
	di "Competing forecasts:  {res:`fcvar1'} versus {res:`fcvar2'}"
	di "over " return(nobs) " observations"
	di in gr "`kchoix' "
	di "S(1) = " %9.4g return(stat1)   "  p-value = " %5.4f return(p)  _n
	di _n
*
end