model_names = {'results_rec_for_xmas',   'XMAS';
                'results_rec_for_mas',         'MAS';
                };


var_names={...
 %%%        var_name        long_name                           transformation
 %%%                                                            0=> nothing
 %%%                                                            1=> mov_sum(4)
 %%%                                                            2=> multx4
           'piZ_obs',       'Core annual Inflation',            1;          
           'gam_YR_obs',    'Non min. GDP annual growth',       1;
           'gam_YCo_obs',   'Min. GDP annual growth',           1;
           'gam_I_obs',     'Tot. investment annual growth', 	1;
           
           'gam_C_obs',     'Priv. cons. annual growth ',	1;           
           'gam_CG_obs',    'Gov. cons. annual growth ', 	1;
           'gam_HN_obs',    'Total hours annual growth ',     	1; 
           'gam_WN_obs',    'Nom. wage annual inflation ',   	1;
           
           'R_obs',         'Nominal interest rate',            2;     
           'rer_obs',       'Real exchange rate',               0;         
           'stb_obs',       'Trade balance/GDP',                0; 
           'xi_obs',        'Risk premia ',                     0;           
          
           'piA_obs',       'Food annual inflation',          	1; 
           'gam_ICo_obs',  	'Min. investment annual growth',  	1;
           'gam_N_obs',     'Employement annual growth',     	1; 
           
           'piO_obs',       'Energy annual inflation',        	1;            
           'gam_IG_obs', 	'Gov. investment annual growth ',  	1;
           'gam_H_obs',   	'Hours per worker annual growth ',  1; 

           }

model_vars=[18 12];
       
%%% COMPUTE FORECAST AND ERRORS  
drop=3;
for kk=1:2  
    modelname=model_names{kk,1};
    eval(['load ' modelname])
    
    modelname=model_names{kk,2};
    
    for nn=1:model_vars(kk)      
        namevar1=var_names{nn,1}; %% name of variable nn
        var1=oo_.SmoothedVariables.(namevar1)+oo_.SmoothedMeasurementErrors.(namevar1); %%actual value of variable nn
        %var1=var1(1+drop);
        rec_forecast=oo_.RecursiveForecast.Mean.(namevar1)(1+drop:end,:)'; %% recursive forecast foregoing first "dropped" obs
        nobs=numel(var1); %% number of observation in data
        nfor=size(rec_forecast,1);%%length of forecasts
        nrec=size(rec_forecast,2);%% number of recursive forecasts
        all_series=NaN(nobs+nfor,1+nrec); %%initialize matrix that will contain all relevant variables
        all_series(1:numel(var1),1)=var1; %%first column is actual data

        
        for ii=1:size(rec_forecast,2) %% "ii" is the number of recursive forecasts available
            inipos=nobs-nrec+ii+1; %% first forecast is done with "nobs-nrec+ii" observations
            endpos=inipos+nfor-1; %% last forecasted date is "nfor-1" after the first one
            all_series(inipos:endpos,ii+1)= rec_forecast(:,ii); %% fill all_series with the forecast
            all_series(1:inipos-1,ii+1)= var1(1:inipos-1); %% fill all_series with the actual data

        end
        
        %%% transform variables
        
        if var_names{nn,3}==0
            all_series4=all_series(4:end,:);
        end
        if var_names{nn,3}==1
            all_series4=all_series(1:nobs+nfor-3,:)...
                        +all_series(2:nobs+nfor-2,:)...
                        +all_series(3:nobs+nfor-1,:)...
                        +all_series(4:nobs+nfor,:);
        end

        if var_names{nn,3}==2
            all_series4=4* all_series(4:end,:);
        end  
        
        error_sq2=NaN(size(all_series4,1),size(all_series4,2)-1); 
           
           
         for ii=1:size(rec_forecast,2) %% "ii" is the number of recursive forecasts available          
            error_sq2(:,ii)=(all_series4(:,ii+1)-all_series4(:,1)).^2;
         end
        
        error_sq2_aux=NaN(nfor,nrec);
        
        for ii=1:nrec
            aux=error_sq2(:,ii);
            aux=aux((isnan(aux)+(aux==0))==0); %% eliminate zeros and NaNs
            error_sq2_aux(1:numel(aux),ii)=aux;
        end
          
        error_sq2=error_sq2_aux';

        for ii=1:nfor
            aux1=error_sq2(:,ii);
            error_sq_exnan=aux1(1-isnan(aux1)==1);
            
            MSE_bs = bootstrp(1000,@mean,error_sq_exnan).^0.5;

            MSE1(ii)=mean(error_sq_exnan)^0.5;
            MSE1_std(ii)=std(MSE_bs);
            
            MSE1_025(ii)=prctile(MSE_bs,2.5);
            MSE1_05(ii)=prctile(MSE_bs,05);
            MSE1_10(ii)=prctile(MSE_bs,10);
            MSE1_90(ii)=prctile(MSE_bs,90);
            MSE1_95(ii)=prctile(MSE_bs,95);  
            MSE1_975(ii)=prctile(MSE_bs,97.5); 
                       
            SD2(ii)=nanstd(all_series4(nobs-nrec+ii+1-3:end,1));
        end

        MSE.(var_names{nn,1}).(modelname)=MSE1;
        MSE_std.(var_names{nn,1}).(modelname)=MSE1_std;
        
        MSE_025.(var_names{nn,1}).(modelname)=MSE1_025;
        MSE_05.(var_names{nn,1}).(modelname)=MSE1_05;
        MSE_10.(var_names{nn,1}).(modelname)=MSE1_10;
        MSE_90.(var_names{nn,1}).(modelname)=MSE1_90;
        MSE_95.(var_names{nn,1}).(modelname)=MSE1_95;
        MSE_975.(var_names{nn,1}).(modelname)=MSE1_975;
        SD.(var_names{nn,1}).(modelname)=SD2;
        
        
        if var_names{nn,3}==0
            all_series4=all_series(4:end,:);
        end
        if var_names{nn,3}==1
            all_series4=all_series(1:nobs+nfor-3,:)...
                        +all_series(2:nobs+nfor-2,:)...
                        +all_series(3:nobs+nfor-1,:)...
                        +all_series(4:nobs+nfor,:);
        end

        if var_names{nn,3}==2
            all_series4=4* all_series(4:end,:);
        end  



        Forecasts.(var_names{nn,1}).(modelname)=all_series4;

    end
end


%%%%%%%%%%%% START PLOTING %%%%%%%%%%%

modelname1=model_names{1,2};
modelname2=model_names{2,2};

dates1=(2002.25:0.25:2019.25)'; 
dates2=cell(size(dates1));

  
