*
* @DMariano( options ) actual f1 f2 start end
* Diebold-Mariano forecast comparison test
*
* Note the Diebold-Mariano test does not apply when the forecasting
* procedures are "nested" (one is a special case of the other).
*
* 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)
*   MODIFIED/[NOMODIFIED]
*      MODIFIED chooses the small sample correction from Harvey, Leybourne and Newbold.
*   TITLE=title for report ["Diebold-Mariano Forecast Comparison Test"]
*   [PRINT]/NOPRINT
*
* Reference:
*   Diebold, Francis X., and Roberto S. Mariano(1995), "Comparing
*    Predictive Accuracy," JBES, v.13, no.3, pp. 253-63.
*
*   Harvey, Leybourne, Newbold(1997), "Testing the equality of prediction mean
*    squared errors", International Journal of Forecasting, vol 13, pp 289-291.
*
* Variables Defined:
*   %CDSTAT  Test statistic for f1=f2 vs f2 better than f1. Flipping the
*      sign gives the test statistic for f2=f1 vs f1 better than f2.
*
* Revision Schedule:
*  08/2003 Written by Tom Doan, Estima
*  06/2010 Modernize code, add TITLE and PRINT options
*  05/2014 Add MODIFIED option
*
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
option string  title
option switch  print     1
option switch  modified  0
*
local integer 	startl endl
local series   crit1  crit2 d
local string   clabel stitle
local real     c1 c2
local report   dreport
local integer  lwindowl
*
*  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. The MDM requires using the TRUNCATED.
*
if modified
   compute lwindowl=3
else
   compute lwindowl=lwindow
*
linreg(robusterrors,lwindow=lwindowl,lags=lags,noprint) d startl endl
# constant
if modified
   *
   * Need /sqrt(%nobs) to correct the formula
   *
   compute %cdstat=%tstats(1)/sqrt(%nobs)*sqrt(%nobs+1-2*lags+lags*(lags-1.0)/%nobs)
else
   compute %cdstat=%tstats(1)
*
if .not.print
   return
*
* Compute the means of the forecast criteria for the two sets of
* forecasts.
*
sstats(mean) startl endl crit1>>c1 crit2>>c2
*
if %defined(title)
   compute stitle=title
else
   compute stitle="Diebold-Mariano Forecast Comparison Test"

report(use=dreport,action=define,title=stitle)
report(use=dreport,atrow=1,atcol=1,span) stitle
report(use=dreport,atrow=2,atcol=1,span) $
   "Forecasts of "+%l(actual)+" over "+%datelabel(startl)+" to "+%datelabel(endl)
if lags>0
   report(use=dreport,atrow=3,atcol=1,span) "Test Statistics Corrected for Serial Correlation of "+lags+" lags"
*
report(use=dreport,atrow=4,atcol=1,align=center) "Forecast" clabel "Test Stat" "P(DM>x)"
report(use=dreport,atrow=5,atcol=1) %l(f1)	  c1      %cdstat   %cdf(-%cdstat)
report(use=dreport,atrow=6,atcol=1) %l(f2)     c2     -%cdstat   %cdf(+%cdstat)
report(use=dreport,action=format,width=10,atcol=2,tocol=2)
report(use=dreport,action=format,picture="*.####",align=decimal,atcol=3,tocol=3,atrow=5)
report(use=dreport,action=format,picture="#.#####",atcol=4,tocol=4)
report(use=dreport,action=show)
end
