function [test, reject] = dm_fsa(d, wceband, wpeband)

% This function performs the tests for equal predictive accuracy using 
% fixed-smoothing asymptotics, as in Laura Coroneoy and Fabrizio Iacone,
% "Comparing predictive accuracy in small samples using fixed-smoothing 
% asymptotics", Journal of Applied Econometrics, 2020.

% INPUTS: 
% d         loss differential (T x 1)
% wceband   bandwidth for Weighted Covariance Estimate (default = T^(1/2))
% wpeband   bandwidth for Weighted Periodogram Estimate (default = T^(1/3))


% OUTPUTS: 

% test     2x1 array with test statistic values: 
%            1) Test statistic with Weighted Covariance Estimate and Bartlett kernel
%            2) Test statistic with Weighted Periodogram Estimate and Daniell kernel

% reject    2x1 array (10 = null of equal predictive accuracy rejected at 10%;
%                   5 = null of equal predictive accuracy rejected at 5%;
%                   0 = null equal predictive accuracy NOT rejected) for the:
%            1) Fixed-b test
%            2) Fixed-m test

test = NaN(2,1);
reject = NaN(2,1);

T = size(d, 1);

num = sqrt(T)* mean(d);

%% Fixed-b test using Weighted Covariance Estimate and Bartlett kernel

if exist ('wceband', 'var') == 0
    wceband = floor(T^(1/2));
end

WCEb = NeweyWest(d,wceband,1);
test(1) = num / sqrt(WCEb);

b = wceband/T;

CV_WCE_b_97 = 1.9600 + 2.9694*b + 0.4160*b^2 - 0.5324*b^3; 
CV_WCE_b_95 = 1.6449 + 2.1859*b + 0.3142*b^2 - 0.3427*b^3; 

if abs(test(1))>CV_WCE_b_97
    reject(1)= 5;
elseif abs(test(1))>CV_WCE_b_95
    reject(1)= 10;
else
    reject(1)= 0;
end
  
%% Fixed-m test using Weighted Periodogram Estimate and Daniell kernel

if exist('wpeband', 'var') == 0
    wpeband = floor(T^(1/3));
end

AA = dftmtx(T);
w = AA*d; % discrete fourier transform of d
pdm = w.* conj(w)/ (2*pi*T);
WPE = sum(pdm(2:wpeband+1,:))/wpeband * 2*pi;

test(2) = num / sqrt(WPE);

if abs(test(2))>icdf('t', 0.975, 2*wpeband)
    reject(2)= 5; 
elseif abs(test(2))>icdf('t', 0.95, 2*wpeband)
    reject(2)= 10;
else
    reject(2)= 0;
end

function omegahat = NeweyWest(Z,nlags,bartlett)

% Returns the Newey-West estimator of the asymptotic variance matrix

% INPUTS
% Z         = nxk matrix 
% nlags     = number of lags
% bartlett  = 1 if bartlett kernel, 0 rectangular

% OUTPUTS: 
% omegahat  = Newey-West estimator of the covariance matrix

[T,~] = size(Z);

% de-mean the variables
Z = Z - ones(size(Z,1),1)*mean(Z);

% sample variance
samplevar = Z'*Z/T; 

omegahat = samplevar;

% sample autocovariances
if nlags > 0
   for ii = 1 : nlags
      gamma = (Z(1+ii:T,:)'*Z(1:T-ii,:) +Z(1:T-ii,:)'*Z(1+ii:T,:))/T;
      if bartlett
          weights = 1 - (ii/nlags);
      else
          weights = 1;
      end
      omegahat = omegahat + weights*gamma;
   end
end
