%% MAIN file (Forecast computation based on routines from Michal Rubaszek and Marcin Kolasa, JMCB, 2012, 44(7):1301-1324)
% SW-DSGE-VAR model Forecasts (unconditional and conditional)
% Point and density forecasts
% Variables to forecast:    dy pinfobs robs
% Full Sample:              1947q3 - 2014q1
% Esimated sample:          1966q2 - 2004q4
% Out-of-sample forecasts:  2005q1 - 2014q1

%%
clear all
close all
clc

% settings
XXXX = 2005;
% model name to forecast - dsgevar(p)
%model_string='sw_model_dsgevar2'; % dsgevar(2)
model_string='sw_model_dsgevar3'; % dsgevar(3)
%model_string='sw_model_dsgevar'; % dsgevar(4)
KK = 3; % years horizon of forecast
horizon = KK*4; % quarterly horizon of forecast
datanew = xlsread('usmodel_newdata.xlsx','W7:AD273');


% loop for rolling -sample dsge-var estimation
% optimal lambda estimated (dvar4) - insample fit
% estimating and forecasting based on different lags
for vv = 0 : 0.25 : 6 % Out-of-sample forecasts:  2005q1 - 2014q1
    
    time = [1947.50 : 0.25 : 2014]';

    % start of the data we want to plot 
    start   = find(time == 2003);
    % last point of the observed sample
    endsam  = find(time == XXXX + vv);
    % last point of the forecast sample
    endfor  = find(time == XXXX + vv + KK);

    dc      = datanew(1 : endsam,  1);
    dinve   = datanew(1 : endsam , 2);
    dy      = datanew(1 : endsam , 3);
    labobs  = datanew(1 : endsam , 4);
    pinfobs = datanew(1 : endsam , 5); 
    dw      = datanew(1 : endsam , 6); 
    robs    = datanew(1 : endsam , 7); 
    
    save usmodel_newdata dc dinve dy labobs pinfobs dw robs
    
    % dsge-var model estimation -- mcmc posterior
    model=['dynare ' model_string ' noclearall']; 
    eval(model)
    %copyfile('sw_model_dsgevar/metropolis/sw_model_dsgevar_mh*',['dvar',num2str(options_.dsge_varlag),'/',num2str(XXXX+vv)]);
    %global M_ dr1_test ex0_ it_ oo_recursive_  bayestopt_  estim_params_  fjac  oo_ options_ ys0_  
    %save(['dvar',num2str(options_.dsge_varlag),'/',num2str(XXXX+vv),'/results.mat']);     
    %rmdir('sw_model_dsgevar','s')
    
%% Find point forecast (forecast horizon=12)

    close all
    global options_ M_
    % loads estimation results, mh settings, observables
    load([M_.fname,'_results.mat']);   
    load([M_.fname,'/metropolis/',M_.fname,'_mh_history_0.mat']);

    % Initializes variables 
    cond_value           = robs(end);               % SPF nowcast conditional on interest rate
    %SPF_nowcast          = [dy(end) dc(end) dinve(end) labobs(end) pinfobs(end) robs(end)];
    cond_options.controlled_varexo='eta_r';
    cond_options.replic  = 1;               % nr of draws for future shocks
    cond_options.periods = horizon;              % forecast horizon
    cond_options.Ydata   = Ydata; 
    cond_options.Xdata   = Xdata;
    
        % reads MH settings (for draws of theta)
        npar = estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn + estim_params_.np;
        draw_.mh_nblck      = options_.mh_nblck;
        draw_.NumberOfDraws = sum(record.MhDraws(:,1))-floor(options_.mh_drop*sum(record.MhDraws(:,1)));
        draw_.fname         = M_.fname;
        draw_.FirstLine     = record.KeepedDraws.FirstLine; 
        draw_.FirstMhFile   = record.KeepedDraws.FirstMhFile; 
        draw_.MAX_nruns     = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
        
        % Forecast is calculated for every k-th step
        k = 20;
        
        % some parameters
        h  = cond_options.periods+1;     % max forecast horizon
        n  = draw_.mh_nblck;
        m  = floor(draw_.NumberOfDraws/k);
        
        % initializes matices in which we will save forecast draw
        cf.dlGDP    = NaN(h,n*m);  f.dlGDP    = NaN(h,n*m); spf_cf.dlGDP    = NaN(h,n*m); 
        cf.dlCONS   = NaN(h,n*m);  f.dlCONS   = NaN(h,n*m); spf_cf.dlCONS   = NaN(h,n*m);
        cf.dlINV    = NaN(h,n*m);  f.dlINV    = NaN(h,n*m); spf_cf.dlINV    = NaN(h,n*m);
        cf.lHOURS   = NaN(h,n*m);  f.lHOURS   = NaN(h,n*m); spf_cf.lHOURS   = NaN(h,n*m);
        cf.dlP      = NaN(h,n*m);  f.dlP      = NaN(h,n*m); spf_cf.dlP      = NaN(h,n*m);
        cf.FEDFUNDS = NaN(h,n*m);  f.FEDFUNDS = NaN(h,n*m); spf_cf.FEDFUNDS = NaN(h,n*m);
        
        % LOOP for MH draws
       for ChainNumber = 1:n
            for DrawNumber  = 1:m
                
                % draw from posterior distribution of theta
                cond_options.parameter_set = ParamPosteriorDraw(draw_, ChainNumber, DrawNumber*k);
                
                % dsge-var forecasts (conditional and unconditional)
                forecasts = Dvars_forecast(cond_value, [], cond_options);
                
                % collects forecasts for all draws
                cnt = DrawNumber+(ChainNumber-1)*m;
                
                cf.dlGDP(:, cnt)        = forecasts.condIR.dy;       cf.dlCONS(:, cnt)   = forecasts.condIR.dc;
                cf.dlINV(:, cnt)        = forecasts.condIR.dinve;    cf.lHOURS(:, cnt)   = forecasts.condIR.labobs;
                cf.dlP(:, cnt)          = forecasts.condIR.pinfobs;  cf.FEDFUNDS(:, cnt) = forecasts.condIR.robs;
                
                f.dlGDP(:, cnt)         = forecasts.uncond.dy;        f.dlCONS(:, cnt)   = forecasts.uncond.dc;
                f.dlINV(:, cnt)         = forecasts.uncond.dinve;     f.lHOURS(:, cnt)   = forecasts.uncond.labobs;
                f.dlP(:, cnt)           = forecasts.uncond.pinfobs;   f.FEDFUNDS(:, cnt) = forecasts.uncond.robs;
        
                %spf_cf.dlGDP(:, cnt)    = forecasts.condALL.dlGDP;    spf_cf.dlCONS(:, cnt)   = forecasts.condALL.dlCONS;
                %spf_cf.dlINV(:, cnt)    = forecasts.condALL.dlINV;    spf_cf.lHOURS(:, cnt)   = forecasts.condALL.lHOURS;
                %spf_cf.dlP(:, cnt)      = forecasts.condALL.dlP;      spf_cf.FEDFUNDS(:, cnt) = forecasts.condALL.FEDFUNDS;
                
                % saves forecast draws
                %save(['dvar',num2str(options_.dsge_varlag),'_',num2str(XXXX+vv),'_forecasts.mat'], 'f', 'cf');
            end
       end
       
       % save dc dinve dy labobs pinfobs dw robs
       actual_data       = [datanew(start:endfor,3),datanew(start:endfor,5),datanew(start:endfor,7)];
       actual_data_block = [datanew(start:endsam,3),datanew(start:endsam,5),datanew(start:endsam,7)];
    
       % unconditional point forecast for dy, inf and int (12 period).
       f.dlGDP(1,:)=[]; % first row end of sample
       f.dlP(1,:)=[];
       f.FEDFUNDS(1,:)=[];
       
       mean_forecast_dy  = mean(f.dlGDP,2);
       mean_forecast_inf = mean(f.dlP,2);
       mean_forecast_int = mean(f.FEDFUNDS,2);
       
       % transformation (annualization) ??
       %if isequal(var_name,'FEDFUNDS')
       %draws    = 4*draws';         
       %else
       %draws     = ((exp(draws'/100)).^4-1)*100;
       %end
       
       % point forecast for plot
       mean_forcasts = [ actual_data_block;...
       mean_forecast_dy,...
       mean_forecast_inf,...
       mean_forecast_int]; 
       
       %% Computes confidence band for the forecast
       yf = [mean_forecast_dy,mean_forecast_inf,mean_forecast_int];
       nstatic = M_.nstatic;
       nspred = M_.nspred;
       dr = oo_.dr;
       nc = size(dr.ghx,2);
       inv_order_var = oo_.dr.inv_order_var;
       [A,B] = kalman_transition_matrix(dr,nstatic+(1:nspred),1:nc,M_.exo_nbr);
       
       var_list = options_.varlist;
       if size(var_list,1) == 0
          %var_list = M_.endo_names(1:M_.orig_endo_nbr,:);
          error (['Specifiy forecast variables in the .mod']);
       end
       nvar = size(var_list,1);
       ivar=zeros(nvar,1);
       for i=1:nvar
           i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
           if isempty(i_tmp)
              disp(var_list(i,:));
              error (['One of the variable specified does not exist']);
           else
              ivar(i) = i_tmp;
           end
       end
       
       ghx1 = dr.ghx(inv_order_var(ivar),:);
       ghu1 = dr.ghu(inv_order_var(ivar),:);
       
       sigma_u = B*M_.Sigma_e*B';
       sigma_u1 = ghu1*M_.Sigma_e*ghu1';
       sigma_y = 0;

       for i=1:horizon
           sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
           var_yf(i,:) = diag(sigma_y1)';
           if i == horizon
              break
           end
           sigma_u = A*sigma_u*A';
           sigma_y = sigma_y+sigma_u;
       end
       
       fact = norminv((1-options_.conf_sig)/2,0,1);
       
       int_width = zeros(horizon,M_.endo_nbr);
       for i=1:nvar
           int_width(:,i) = -fact*sqrt(var_yf(:,i));
       end       
       
       var_name = strvcat('GDP','Inflation','Interest');
       
       maximum_lag = M_.maximum_lag;
       
       for i=1:nvar
       eval(['forbands.HPDinf.' var_name(i,:) '= yf(:,' int2str(i) ')-' ...
          ' int_width(:,' int2str(i) ');']);
       eval(['forbands.HPDsup.' var_name(i,:) '= yf(:,' int2str(i) ')+' ...
          ' int_width(:,' int2str(i) ');']);
       end
       
       % Confidence bands for plot
       HPDinf_forcasts = [ actual_data_block;...
       forbands.HPDinf.GDP,...
       forbands.HPDinf.Inflation,...
       forbands.HPDinf.Interest];  

       HPDsup_forcasts = [ actual_data_block;...
       forbands.HPDsup.GDP,...
       forbands.HPDsup.Inflation,...
       forbands.HPDsup.Interest]; 
       
       %% plot
       % Fan charts
       strout = ['dvar',num2str(options_.dsge_varlag),'_forecast_output_' num2str(XXXX+vv)];
       mkdir(strout);
       
       
       for var = 1 : size(actual_data,2)
           figure('Name',['SW Forecast of ' deblank(var_name(var,:))]);
           h = area([time(start:endfor)],[HPDinf_forcasts(:,var), HPDsup_forcasts(:,var) - HPDinf_forcasts(:,var)]);%,'FaceColor',[.85 .85 .85]);
           set(h(2),'FaceColor',[.95 .95 .95])
           set(h(1),'FaceColor',[1 1 1])
           hold on;
           plot(time(start:endfor),mean_forcasts(:,var),'--r','Linewidth',2);
           hold on;
           plot(time(start:endfor),actual_data(:,var),'k');
           hold on
           axis tight
           title(var_name(var,:),'fontsize',16);    
           ext  = 'eps';
           type = 'psc2';
           saveas(gcf,['.\' strout '\DVAR', num2str(options_.dsge_varlag), '_' num2str(XXXX + vv) ...
                  '_' deblank(var_name(var,:)) '.' num2str(ext)],num2str(type))
           ext='fig';
           type='fig';
           saveas(gcf,['.\' strout '\DVAR', num2str(options_.dsge_varlag), '_' num2str(XXXX + vv) ...
                  '_' deblank(var_name(var,:)) '.' num2str(ext)],num2str(type))
       end
       
       % forecast errors with transformation (annualization)
       close all
       actual = [datanew(endsam+1:endfor,3),datanew(endsam+1:endfor,5),datanew(endsam+1:endfor,7)];
       forecast =  [mean_forecast_dy, mean_forecast_inf,mean_forecast_int]; 
       forerror(:,:,4*vv+1)  = actual - forecast;
       forsqerr(:,:,4*vv+1)  = (actual - forecast).^2;
       
end

%% Calculate ME, MAE and RMSE  
close all
absfor    = abs(forerror);
forME     = squeeze(mean(forerror,3));
forMAE    = squeeze(mean(absfor,3));
forRMSE   = squeeze(mean(forsqerr,3).^(0.5));

save RMSEdvar forME forMAE forRMSE

%% Plots RMSEs from DSGE and BVAR forecasts - need rmse.mat files (optional)
load rmse % load results from dsge and bvar forecast rmse
 for var = 1 : size(actual_data,2)
        figure('Name',['RMSE of ' var_name(var,:)])
        plot(rmse_sw(:,var),'b','Linewidth',2); hold on;
        plot(rmse_bvar(:,var),'r','Linewidth',2); hold on;
        plot(forRMSE(:,var),'k','Linewidth',2); hold on;
        axis tight
        ext  = 'eps';
        type = 'psc2';
        saveas(gcf,['rmse_' deblank(var_name(var,:)) '.'   num2str(ext)],num2str(type))
        ext='fig';
        type='fig';
        saveas(gcf,['rmse_' deblank(var_name(var,:)) '.'   num2str(ext)],num2str(type))
 end