function [dm,pdm] = dbmar(err1,err2,h,lags,loss,wind)
% DIEBMAR.m 
% Purpose: computes the Diebold Mariano statistics with small sample
% correction (Harvey, L and Newbold (1997): compare with a t-stat with N-1
% degrees of freedom
%
% Inputs
% -------  
% err1: the prediction error of model 1            (T-by-1 vector)   
% err2: the prediction error of model 2            (T-by-1 vector)   
% h   : the forecasting horizon                    (scalar)
% lags: the maximum number of lags used to compute 
%       the long run variance (it must be > h-1)   (scalar)
% loss: type of loss function                      (string)
%       'SE' --->  squared error loss (default, if omitted)           
%       'AD'  ---> absolute error loss 
% wind: type of windown for correcting the long-run variance estimator (string)
%       'triang' --->  Newey-West triangular (default, if omitted)           
%       'rect'   --->  equal weights truncated at h-1
% Outputs
% -------  
% dm : the Diebold Mariano statistics              (a scalar) 
% pdm: the p-value associated to the statistics    (a scalar)   
%
% Example:
% --------  
% [dm,pdm] = dbmar(err1,err2,4,5,'AD')  
  
if lags<h-1
  error('lags must be larger than h-1!')
end;
  
if (lags~=h-1) && strcmp(wind,'rect')
%     disp('with rectangular window I set the lags to h-1')
    lags = h-1;
end;

losses = [{'SE'},{'AD'}];

if nargin<4
  error('Not enough arguments');
elseif nargin==4
  typeloss = 1;
  wind = 'triang';
else
  closs = {loss};
  typeloss = find(strcmp(losses,closs));
end;

switch typeloss
 case 1
  lo1 = err1.^2;
  lo2 = err2.^2;
 case 2
  lo1 = abs(err1);
  lo2 = abs(err2);
end;

dt = lo1-lo2;
mdt = mean(dt);

% make it a time series to be used by bits
lrcov = myxcov(dt,dt,0);
%lrcov = cov(dt);
if strcmp(wind,'rect');
    if h>1
        autocov = myxcov(dt,dt,lags);
        autocov = autocov(1:lags);
        correz0 = 1;
        lrcov = lrcov + 2*sum((correz0*autocov));
    end;
elseif strcmp(wind,'triang')
    if h>1
        autocov = myxcov(dt,dt,lags);
        autocov = autocov(1:lags);
        correz0 = 1:-1/lags:1/lags;
        correz0 = 1-correz0;
        lrcov = lrcov + 2*sum((correz0*autocov));
    end;
end
T = max(size(err1));

dm = mdt/sqrt(lrcov/T);
pdm = 2*(1-normcdf(abs(dm),0,1));  