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

% 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, forthcoming.

% 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))
% tsize     level of significance of the two-sided test (5% (default) or 10%)

% 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 (0 = null of equal predictive accuracynot rejected;
%                     1 = null equal predictive accuracy 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);

if exist('tsize', 'var') == 0
    tsize = 0.05;  
end

%% 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;

if tsize == 0.05
    crit = 1.9600 + 2.9694*b + 0.4160*b^2 - 0.5324*b^3; 
elseif tsize == 0.1
    crit = 1.6449 + 2.1859*b + 0.3142*b^2 - 0.3427*b^3; 
end

if abs(test(1)) > crit
    reject(1) = 1;
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',1-tsize/2,2*wpeband)
    reject(2) = 1;
else
    reject(2) = 0;
end

%% Results
printmat([test reject], ['Tests for equal predictive accuracy at ',...
    num2str(tsize*100),'% significance level'], 'WCE-B_Fix-b WPE-D_Fix-m', 'Statistic Reject');

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
