clear;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%OPTIONS TO BE CHANGED BY THE USER%%%%%%%%%%%%%%%%
RESULT_names={'NK_SW_results','GK_equity_SW_results', 'BGG_normal_dist_SW_results','BGG_uniform_dist_SW_results' }; %%EXACT NAMES OF THE RESULTS FILES TO LOAD
lin_vs_nonlin=[0,0,0,0]; %=0 IF THE CORRESPONDING MODEL IN RESULT_NAMES IS NON-LINEAR, 1 IF IT IS LINEAR
VAR_IRFs_nonlin = {'Y', 'C', 'I', 'H', 'W', 'R',  'Rn', 'Q','PIE', 'RK', 'spread', 'N'}; %% NAMES OF ENDOGENOUS VARIABLES OF INTEREST FOR THE IRFS (NON-LINEAR MODELS) - NEED TO HAVE EXACTLY SAME NAME AS IN THE .MOD FILE
VAR_IRFs_linear= {'Y','C','I','H','W','R','Rn','Q','PIE'};
NAME_SHOCKS={'_epsA','_epsAtrend','_epsG','_epsMC','_epsMPS','_epsMRSS', '_epsRPS', '_epsIS' };%% NAMES OF EXOGENOUS SHOCKS OF INTEREST FOR THE IRFS - NEED TO HAVE EXACTLY SAME NAME AS IN THE .MOD FILE
names = char('Output', 'Consumption', 'Investment', 'Hours worked', 'Real wage', 'Real interest rate',  'Nominal interest rate', 'Tobin Q', 'Inflation','Rental Rate',...
'spread', 'Net Worth'    ); % NAMES OF THE EDOGENOUS VARIABLES FOR THE GRAPHS. YOU CAN SPECIFY ANY NAME YOU WANT BUT NEED TO HAVE THE SAME NUMBER OF ELEMENT OF VAR_IRFs
Model_names=char('NK ','GK', 'BGG (normal)','BGG (uniform)' ); %NAMES OF THE MODEL VARIANTS TO APPEAR IN THE LEGEND OF THE GRAPHS
%Model_names=char('NK (No Habit)','NK (External Habit)','NK (Internal Habit)'); %NAMES OF THE MODEL VARIANTS TO APPEAR IN THE LEGEND OF THE GRAPHS
Rows_figure=3; % NUMBER OF ROWS IN EACH PLOT
Column_figure=4; % NUMBER OF COLUMNSS IN EACH PLOT
%Rows_figure=2; % NUMBER OF ROWS IN EACH PLOT
%Column_figure=5; % NUMBER OF COLUMNSS IN EACH PLOT
irf_horizon=40; %LENGTH OF THE IRFS GENERATED BY THE .MOD FILE
%%%%%%%%%%%%%%%%%%%%%%%END OF OPTIONS%%%%%%%%%%%%%%%%%%%%%%%%%%



%% DO NOT CHANGE ANYTHING BELOW.
Num_models=length(RESULT_names); %%NUMBER OF MODELS TO COMPARE
NUM_SHOCKS=length(NAME_SHOCKS);  %%NUMBER OF SHOCKS TO COMPARE
NUM_VAR=length(VAR_IRFs_nonlin);%%NUMBER OF VARIABLES TO COMPARE

irfs_matrix=zeros(irf_horizon,NUM_VAR,NUM_SHOCKS,Num_models);
for mm=1:Num_models;
%load Model 
eval(['load ' RESULT_names{mm} ';']);
if lin_vs_nonlin(mm)~=0
    VAR_IRFs=VAR_IRFs_linear;
else
    VAR_IRFs=VAR_IRFs_nonlin;
end
for xx=1:NUM_SHOCKS;
for jj=1:NUM_VAR;
%Rename the IRFs for each variable of interest 
genrate_irf_names=[VAR_IRFs{jj},NAME_SHOCKS{xx},'=','oo_.irfs.',VAR_IRFs{jj},NAME_SHOCKS{xx},';'];
evalin('base', genrate_irf_names);
%generate_irf_matrix=['irf_',num2str(mm),NAME_SHOCKS{xx},'(jj,:)=',VAR_IRFs{jj},NAME_SHOCKS{xx}];
%evalin('base', generate_irf_matrix)
irfs_matrix(:,jj,xx,mm)=eval([VAR_IRFs{jj},NAME_SHOCKS{xx}]);
end
end

end



%%PLOT


color = {'-k','--r','.-g','+b','*y'};%color and line style for the plots


for xx=1:NUM_SHOCKS;

%Options for the plot
h=figure('Position', [600, 0, 1000, 900]);
axes ('position', [0, 0, 1, 1]);

%Figure 
%Loop over the number of endogenous variables to plot
F1=figure(xx);
set(F1, 'numbertitle','off');
set(F1, 'name', ['Impulse response functions to',NAME_SHOCKS{xx}]);
for jj = 1:length(VAR_IRFs);
    for mm=1:Num_models;
                subplot(Rows_figure,Column_figure,jj), plot(irfs_matrix(:,jj,xx,mm),color{mm},'LineWidth',2); hold on;
                    xlabel('Quarters');
                    ylabel('% dev from SS');
                    grid on
                    title(names(jj,:),'FontSize',10)

axis tight;                   
    end;
end;
legend(Model_names)
end;