st: dmariano test for significance regarding forecast accuracy of two competing methods

[Oliver Jones <ojones@wiwi.uni-bielefeld.de>]

Hi everybody,

did anyone configure the -dmariano- (from SSC) command
to be able to apply other criteria than mean squared error (MSE) and
mean absolute error (MAE)?
If so, would you like to share your code?

I gave it a try and included some lines in the .ado file which generate the
mean absolute percentage errors (MAPE).


My first question goes to the statisticians:
Dose the dmariano test procedure change when MAPE is used instead of MSE or MAE?

My second question goes to the programmers:
If the statisticians confirm that the test procedure dose not change, is my
configuration of -dmariano- correct?

************* Begin mydmariano.ado *******************
*! version 1.0.2  04Feb2003  CFBaum
* 1.0.1: 04feb2003
* 1.0.2: 27jun2006, corrected for calculation of autocovariances using Mata

program define mydmariano, rclass
	version 9.2
	syntax varlist(ts min=3 max=3) [if] [in] [, Crit(string) Maxlag(integer -1) Kernel(string)]
	marksample touse
	_ts timevar, sort
	markout `touse' `timevar'
* MSE default criterion
	local criterion "MSE"
	if upper("`crit'") == "MSE" | upper("`crit'") == "MAE" | upper("`crit'") == "MAPE" {
		local criterion = upper("`crit'")
		}
	else if "`crit'" != "" {
		dis as err ///
			"The accuracy measure you specified is not available."
		exit
		}
* Uniform default kernel
	local window "uniform"
	if lower("`kernel'") == "uniform" | lower("`kernel'") == "bartlett" {
		local window = lower("`kernel'")
		}
	tempvar  e1 e2 d count
	tempname lrvar T J s1 gamma
* get the three variables
	tokenize `varlist'
	local depvar `1'
	local fcvar1 `2'
	local fcvar2 `3'
* set maxlag
	qui	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'"
		}
* calculate e1 and e2
	if "`criterion'" == "MSE" {
		qui gen double `e1' = (`depvar' - `fcvar1')^2 if `touse'
		qui gen double `e2' = (`depvar' - `fcvar2')^2 if `touse'
		}
	else if "`criterion'" == "MAE" {
		qui gen double `e1' = abs(`depvar' - `fcvar1') if `touse'
		qui gen double `e2' = abs(`depvar' - `fcvar2') if `touse'
		}
	else if "`criterion'" == "MAPE" {
		qui gen double `e1' = abs(`depvar' - `fcvar1')/`depvar' if `touse'
		qui gen double `e2' = abs(`depvar' - `fcvar2')/`depvar' if `touse'
		}
	qui {
		gen double `d' = `e1' - `e2' if `touse'
		summ `e1' if `touse', meanonly
		local e1bar = r(mean)
		summ `e2' if `touse', meanonly
		local e2bar = r(mean)	
		summ `d' if `touse', meanonly
		local T = r(N)
		local dbar = r(mean)
		}

* generate autocovariances of d series
	local varlist2 L(0/`maxlag').`d'
	tsrevar `varlist2'
	local varlist3 `r(varlist)'
	local ml1 = `maxlag' + 1
	mata:  autocov("`varlist3'",`dbar',"`touse'")
    scalar `J' = __gamma[1,1]
    forv l = 1/`maxlag' {
* uniform kernel is unweighted; otherwise use Bartlett kernel
        local w 1
        if "`window'" == "bartlett" {
        	local w = 1 - (`l'/(`maxlag'+1))
        	}
        scalar `J' = `J' + 2*`w'*__gamma[`l'+1,1]
        }
    scalar `lrvar' = sqrt(`J'/`T')
    if `lrvar' == . {
    	di _n "Long-run variance is non-positive for this kernel and truncation lag."
    	error 506
    	}
*
	qui {
		tsset
		summ `r(timevar)' if `touse'
		}	
	return scalar e1bar = `e1bar'
	return scalar e2bar = `e2bar'
 	return scalar dbar = `dbar'
 	return scalar lrsd = `lrvar'
	return scalar s1 = `dbar'/`lrvar'
	return scalar p = 2*normprob(-abs(return(s1)))
	return scalar nobs = `T'
	return scalar maxlag = `maxlag'
	return local criterion  `criterion'
	return local kernel  `window'
	di _n "Diebold-Mariano forecast comparison test for actual : {res:`depvar'}"
*	di in gr "Forecast period : " %tm r(min) " -" %tm r(max) "  N : `nobs'"
	di "Competing forecasts:  {res:`fcvar1'} versus {res:`fcvar2'}"
	di "Criterion: " return(criterion) " over " return(nobs) " observations"
	di in gr "`kmax'   Kernel : " return(kernel)
	di in gr _n "Series" _col(24)  return(criterion)
 	di in gr _dup(30) "_"
	di in ye "{res:`fcvar1'}" _col(20) in gr %9.4g `e1bar'
	di in ye "{res:`fcvar2'}" _col(20) in gr %9.4g `e2bar'
	di "Difference" _col(20) %9.4g `dbar'
	local better `fcvar1'
	if `e1bar'>`e2bar' {
		local better `fcvar2'
		}
	di _n "By this criterion, {res:`better'} is the better forecast"
	di "H0: Forecast accuracy is equal."
	di "S(1) = " %9.4g return(s1)   "  p-value = " %5.4f return(p)  _n
*	di "Forecast comparison statistic : " %9.4g
*	di "H0 : no difference in accuracy, distributed N(0,1) with pval = " %5.4f return(p)
	di _n
*
end

mata:
mata set matastrict on
void autocov(string scalar vname, real scalar dbar, string scalar touse)
{
        real matrix X, gamma
        real scalar T
        string rowvector vars
        string scalar v
   // access the Stata variables in varlist, respecting touse
        vars = tokens(vname)
        v = vars[|1,.|]
        st_view(X,.,v,touse)
// demean by dbar
		X = X :- dbar
// change all missing values to 0, which will be ignored in cross
		_editmissing(X,0)
// apply cross to get the covariances and scale by T
        gamma = cross(X,X) :/ strtoreal(st_local("T"))
// return the __gamma matrix
		st_matrix("__gamma",gamma)
}
end
************* End mydmariano.ado *******************


Thanks,
Oliver
--
Universität Bielefeld
Fakultät für Wirtschaftswissenschaften
Lehrstuhl für Ökonometrie und Statistik
- -
Bielefeld University
Faculty of Business Administration and Economics
Chair of Econometrics and Statistics
- -
Raum / room:  V9-110
Tel / phone:  +49 (0)521 106 4871
---
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/