%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script uses the forecast errors generated
% by Compute_Errors.m to evaluate forecast bias for
% each variable and forecast horizons, controlling for recessions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

methods = fieldnames(RT_DATA);

for m = 1:size(methods,1) % Iterate over forecast methods

    
    variables=fieldnames(REV_DATA.(sprintf(cell2mat(methods(m,1)))));

    for i = 1:size(variables,1) % Iterate over variables

        
        jumps = fieldnames(REV_DATA.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))));        
        
        for j = 1:size(jumps,1) % Iterate over jump-off delays        

            vintages = fieldnames(REV_DATA.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))));
            
            ERRORS_DATA.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1))))(:,:) = NaN(size(vintages,1),PARAMS.horizons);
            
            for k = 1:size(vintages,1) % Iterate over vintages

                
                max_horizon = size(REV_DATA.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))). ...
                    (sprintf(cell2mat(jumps(j,1)))).(sprintf(cell2mat(vintages(k,1)))),2);

                for l = 1:max_horizon % Iterate over horizons

                    % Organise data

                    ERRORS_DATA.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1))))(k,l) = ...
                        REV_DATA.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))).(sprintf(cell2mat(vintages(k,1))))(1,l);
                    

                    
                    DUMMY_VECTOR(k,l) = PARAMS.recession_dummy.(sprintf(cell2mat(vintages(k,1))))(l,1);
                    

                end
                
            end


            %Evaluate bias

            for l = 1:PARAMS.horizons % Iterate over horizons

                % Regress the forecast errors on a dummy and calculate the
                % Newey-West standard errors with lags = horizon - 1

                Y = ERRORS_DATA.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1))))(:,l);
                Y=Y(1:length(Y)-sum(isnan(Y)),1)
                X = [cell2mat(DUMMY_VECTOR(1:length(Y),l))];
                
                               
                Mdl = fitlm(X,Y);
                [WhiteCov,WhiteSe,coeff] = hac(Mdl,'bandwidth',1);

                % Calculate significance of bias and store statistics

                df = Mdl.DFE;

                ERRORS_STATS.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))). ...
                    REG_COEFF_CONST_control(1,l) = coeff(1);
                
                ERRORS_STATS.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))). ...
                    REG_SE_control(1,l) = WhiteSe(1);
                
                ERRORS_STATS.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))). ...
                    REG_TSTAT_control(1,l) = coeff(1)/WhiteSe(1);
                
                ERRORS_STATS.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))). ...
                    REG_PVAL_control(1,l) = 2*(1-tcdf(abs(coeff(1)/WhiteSe(1)),df));
                
                ERRORS_STATS.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))). ...
                    REG_PVAL_SIGN_control(1,l) = (1-2*(1-tcdf(abs(coeff(1)/WhiteSe(1)),df)))*-sign(coeff(1));
                
                ERRORS_STATS.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))). ...
                    REG_CRIT90_control(1:2,l) = tinv([0.050  0.950],df)'*WhiteSe(1);
                
                ERRORS_STATS.(sprintf(cell2mat(methods(m,1)))).(sprintf(cell2mat(variables(i,1)))).(sprintf(cell2mat(jumps(j,1)))). ...
                    REG_CRIT95_control(1:2,l) = tinv([0.025  0.975],df)'*WhiteSe(1);

            end
            
        end

    end

end

clearvars -except DATA_level PARAMS FC_DATA RT_DATA FINAL_DATA* REV_DATA* REV*_DATA *MSE* ERRORS_STATS* *RATIO* *L_DATA


