%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script uses the forecast errors generated
% by Cmopute_Errors.m to compare Staff forecasts to benchmarks
% using various statistics and tests
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

methods = fieldnames(RT_DATA);

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

    
    variables = fieldnames(RT_DATA.(sprintf(cell2mat(methods(i,1)))));

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

        
        jumps = fieldnames(RT_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))));
            
        for k = 1:size(jumps,1) % Iterate over jump-off delays
            
                                        
            if ~strcmp(methods(i),'Staff')    

                if strcmp(jumps(k),'jump_final') % For last vintage defined as final and benchmark input
                    
                    
                    rmse_ratio.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1)))) = ...
                        RMSE_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))./ ...
                            RMSE_DATA.Staff.(sprintf(cell2mat(variables(j,1)))).jump_final;



                    for l = 1:PARAMS.horizons

                        comp1 = REV1_DATA.Staff.(sprintf(cell2mat(variables(j,1)))).jump_final(1:end,l);

                        comp2 = REV1_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1:end,l);
                        
                        comp1 = (comp1(~isnan(comp1)));
                        
                        comp2 = (comp2(~isnan(comp2)));

                        dm_stats.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1,l) = ...
                            dmtest(comp1,comp2,l-1+2);

                        sup_model_rate.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1,l) = ... % 1 if model beats staff, 0 otherwise
                            mean(abs(comp1)>abs(comp2));
                        
                        hln_staff_model_stats.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1,l) = ...
                            hlntest(comp1,comp2,l-1+2);
                        
                        hln_model_staff_stats.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1,l) = ...
                            hlntest(comp2,comp1,l-1+2);
                        

                    end
                    
                    
                else

                    rmse_ratio.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1)))) = ...
                        RMSE_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))./ ...
                            RMSE_DATA.Staff.(sprintf(cell2mat(variables(j,1)))).jump_0;   



                    for l = 1:PARAMS.horizons

                        comp1 = REV1_DATA.Staff.(sprintf(cell2mat(variables(j,1)))).jump_0(1:end,l);

                        comp2 = REV1_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1:end,l);
                        
                        comp1 = (comp1(~isnan(comp1)));
                        
                        comp2 = (comp2(~isnan(comp2)));

                        dm_stats.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1,l) = ...
                            dmtest(comp1,comp2,l-1+2);

                        sup_model_rate.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1,l) = ... % 1 if model beats staff, 0 otherwise
                            mean(abs(comp1)>abs(comp2));
                        
                        hln_staff_model_stats.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1,l) = ...
                            hlntest(comp1,comp2,l-1+2);
                        
                        hln_model_staff_stats.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(1,l) = ...
                            hlntest(comp2,comp1,l-1+2);

                    end
                
                end

            end

        end

    end

end

clearvars -except PARAMS DATA* FC_DATA RT_DATA FINAL_DATA* REV_DATA* REV*_DATA ERRORS_STATS* *MSE* *RATIO* *L_DATA rmse_ratio dm_stats hln*stats sup_model_rate VAR_* COV_*


