*
* @DMariano( options ) actual f1 f2 start end
* Diebold-Mariano forecast comparison test
*
* Parameters:
*   actual     series of actual data
*   f1         first set of forecasts
*   f2         second set of forecasts
*   start end  range to analyze (common range of f1, f2 and actual by default)
*
* Options:
*   CRITERION=[MSE]/MAE
*     Chooses whether the error criterion will be squared or absolute errors
*   LWINDOW=NEWEYWEST/BARTLETT/[TRUNCATED]/PARZEN/QUADRATIC
*     Chooses the lag window form (NEWEYWEST and BARTLETT are synonyms).
*   LAGS=number of lags or bandwidth (for QUADRATIC)
*
* Reference:
*  Diebold, Francis X., and Roberto S. Mariano(1995), "Comparing Predictive Accuracy,"
*  JBES, v.13, no.3, pp. 253-63.
*
* Revision Schedule:
*  08/2003 Written by Tom Doan, Estima
*
procedure DMariano actual f1 f2 start end
type series actual f1 f2
type integer start end
*
option choice criterion 1 mse mae
option choice lwindow   3 newey bartlett truncated parzen quadratic
option integer lags     0
*
local integer startl endl
local series crit1 crit2 d
local string clabel
local real c1 c2 dmstat
*
*  Figure out the range if not supplied
*
inquire(reglist) startl<<start endl<<end
# actual f1 f2
*
*  Compute the criterion for the forecast errors
*
if criterion==1 {
   set crit1 startl endl = (actual-f1)**2
   set crit2 startl endl = (actual-f2)**2
   compute clabel="MSE"
}
else {
   set crit1 startl endl = abs(actual-f1)
   set crit2 startl endl = abs(actual-f2)
   compute clabel="MAE"
}
*
* Compute the gap
*
set d startl endl = crit1-crit2
*
* The Diebold-Mariano statistic is simply the "t-stat" on a regression of d on a
* constant with HAC standard errors
*
linreg(robusterrors,lwindow=lwindow,lags=lags,noprint) d startl endl
# constant
compute dmstat=%tstats(1)
*
* Compute the means of the forecast criteria for the two sets of forecasts
*
stats(noprint) crit1 startl endl
compute c1 = %mean
stats(noprint) crit2 startl endl
compute c2 = %mean
*
display "Diebold-Mariano Forecast Comparison Test"
display "Forecasts of" %l(actual) "over" %datelabel(startl) "to" %datelabel(endl)
if lags>0
   display "Test Statistics Corrected for Serial Correlation of" lags "lags"
*
display "Forecast 1:" %l(f1) clabel+"=" c1
display "Test Statistic" dmstat "P-Value for 1=2 vs 1 better than 2" %cdf(dmstat)
display "Forecast 2:" %l(f2) clabel+"=" c2
display "Test Statistic" -dmstat "P-Value for 1=2 vs 2 better than 1" %cdf(-dmstat)
compute %cdstat=dmstat
end

