clc; clear;

addpath('ADD_PATH_HERE\Workfiles');
addpath('ADD_PATH_HERE\Scripts');
addpath('ADD_PATH_HERE\Scripts\RossiSekhposyan');

cd 'ADD_PATH_HERE'

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data

load RT_data
%%

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
            
            
            vintages = fieldnames(FINAL_DATA.(sprintf(cell2mat(variables(j,1)))));
            
            FORECAST_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1)))) = ...
                nan(size(vintages,1),PARAMS.horizons);
            

            
            for l = 1:size(vintages,1) % Iterate over vintages
                
                
                horizons = size(FINAL_DATA.(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(vintages(l,1)))),2);
                
                for h = 1:horizons % Iterate over horizons  
                                          
                        FORECAST_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(l,h) = ...
                            RT_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1)))).(sprintf(cell2mat(vintages(l,1))))(1,h);
                         
                        REALIZED_DATA.(sprintf(cell2mat(variables(j,1))))(l,h) = ...
                            FINAL_DATA.(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(vintages(l,1))))(1,h);
                        
                        ERROR_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(l,h) = ...
                            REALIZED_DATA.(sprintf(cell2mat(variables(j,1))))(l,h)-FORECAST_DATA.(sprintf(cell2mat(methods(i,1)))).(sprintf(cell2mat(variables(j,1)))).(sprintf(cell2mat(jumps(k,1))))(l,h);
                            
                end
            end
            
            REALIZED_DATA.(sprintf(cell2mat(variables(j,1))))(REALIZED_DATA.(sprintf(cell2mat(variables(j,1))))==0) = NaN;
            
        end     
    end
end

%% Loop over all variables and horizons to calculate statistics

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

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

    horizons = 9;

    for h = 1:horizons % Iterate over horizons. The algorithm can't compile for the nowcast due to colinearity.
                
        j
        h
        
        m = 20; % Five years
        forecasterror = ERROR_DATA.Staff.(sprintf(cell2mat(variables(j,1)))).jump_0(1:end,h);
        x1 = ones(size(forecasterror)); % Use only a constant (no additional predictors)
        x2 = [];
                
        results = FluctuationRationalityRecNW(m,forecasterror,x1,x2);
        
        tv_bias.(sprintf(cell2mat(variables(j,1)))).(strcat('h',num2str(h))).coef = results.coef;
        tv_bias.(sprintf(cell2mat(variables(j,1)))).(strcat('h',num2str(h))).wald = results.wseries;
        tv_bias.(sprintf(cell2mat(variables(j,1)))).(strcat('h',num2str(h))).cv_5 = results.cv;
        
    end

end


%% Loop over all variables and other horizons to produce charts

variables = fieldnames(RT_DATA.(sprintf('VAR')));

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

    horizons = 9;
    
    figure
    pbaspect([1 2 1])
    num = 1;
    
    for m = 2:horizons % Iterate over horizons  
                    
            subplot(4,2,num) 
            
            num = num+1;
            
            h=m-1;
            a = 1;
            time = datetime(1982,06,30):calquarters(1):datetime(2015,12,31);
            x=m:size(time,2)-(PARAMS.horizons-m);
            critline=ones(1,size(x,2))*tv_bias.(sprintf(cell2mat(variables(j,1)))).(strcat('h',num2str(m))).cv_5;
            startdate = datenum(time(m));
            enddate = datenum(time(x(end)));
            lenght = size(x,2);
            dates = linspace(startdate,enddate,lenght);

            WALD = plot(dates,[nan(19,1);tv_bias.(sprintf(cell2mat(variables(j,1)))).(strcat('h',num2str(m))).wald],'-k','Color',[0.6 0.6 0.6],'LineWidth',1.5);
            hold on;
            CRIT = plot(dates,critline,'-k','LineWidth',0.75);
            hold on;
            if m==1
                title({'nowcast'},'Fontsize',10);
            else
                title({char(strcat('horizon ',' {\itt}','+',num2str(m-1)))},'Fontsize',10);
            end
            ylabel('Wald stat.','Fontsize',10);
            
            datetick('x','yyyy','keeplimits');

            set(gca,'XLim',[datenum(1982,03,30) datenum(2015,12,30)]);
              
    end

        saveas(gcf,char(strcat('Figure_Appendix_E_',variables(j))),'epsc');   

end