% Plot and save impulse responses
function doing_plots_a(model_names,Impulse, Response, plot_options)   

    plot_options.marks      = {'-k','-.k',':k','--k'};
    plot_options.horizon    = 20;
    plot_options.grid       = 0;
    plot_options.latex      = 1;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot_options.model_names= model_names(:,1)';
legend_names            = model_names(:,2)';

shocks                  = Impulse(:,1)';
shocks_latex            = Impulse(:,2)';
shocks_scale            = Impulse(:,3)';
shocks_long_name        = Impulse(:,4)';

v_sel                   = Response(:,1)';
v_name                  = Response(:,2)';
plot_options.v_do_cusum = cell2mat(Response(:,3)');
plot_options.v_adj      = cell2mat(Response(:,4)');
plot_options.v_div_ss   = cell2mat(Response(:,5)');

plot_options.v_do_cusum = [plot_options.v_do_cusum 	0];  	% shock doesn't sum over
plot_options.v_div_ss   = [plot_options.v_div_ss 	1];    	% shock divided by its ss
plot_options.v_adj_orig = [plot_options.v_adj  		1]*100;	% no adjustment for shock

for ii=1:length(plot_options.model_names)
    load([plot_options.model_names{ii} '.mat'])
    eval(['data.M_.' plot_options.model_names{ii} '=M_;'])
    eval(['data.oo_.' plot_options.model_names{ii} '=oo_;'])
end


for j=1:size(shocks,2)
    %disp(['plotting IRFs for ' shocks{j} ' (' num2str(j) '/'  num2str(size(shocks,2)) ')']);
    plot_options.u_sel=shocks(j);
    plot_options.u_title=shocks_long_name(j);
    if size(plot_options.model_names,2)>1
        for i=2:size(plot_options.model_names,2)
            plot_options.u_sel=[plot_options.u_sel shocks(j)];
        end
    end
    plot_options.v_sel  = [v_sel shocks(j)];
    plot_options.u_name = shocks_latex(j);
    plot_options.v_name = [ v_name shocks_latex(j)];
    plot_options.v_adj  = shocks_scale{j} * plot_options.v_adj_orig;
    

    plots_for_dynare(plot_options,data)
    h=legend(legend_names);
    set(h,'Position',[0.4 0 0.23 0.03],...
        'Orientation','Horizontal',...
        'Fontsize',8,...
        'Interpreter','Latex');
    legend('boxoff')
    
    
end


function plots_for_dynare(plot_options,data)

% This function creates a graph that combines the impulse responses of
% different models and/or different shocks for teh same model, generates by
% Dynare. If the responses come from different models, it is assumed that
% all the models have the same names for variables and shocks in Dynare.
% The input plot_options is a structure that has to include the following:
%
%   plot_options.model_names: A set of strings containing the mat-file names
%       (without includeing the .mat extension) were the oucomes from dynare 
%       have been stored. The function will create a line in the graph for 
%       each string in model_names.
%   plot_options.marks: A set of strings containg the mark-tipe desired for
%       the line corresponding to each line from model_names.
%   plot_options.u_sel: A set of strings of the same dimension as model_names
%       indicating the name (as declared in dynare) of the exogenous variable whose shock 
%       generates the responses to be ploted.
%   plot_options.v_sel: A set of strings containing the name of the
%       variables (as declared in dynare) whose responses will be displayed.
%       The variables toe be plated will be the same for all the alements in
%       model_name. 
%  plot_options.v_adj: A line vector with the same number of elements as
%       strings in v_sel, indicating the scale adjustment that will be
%       applied for each variable in v_sel.
%  plot_options.v_do_cusum: A line vector with the same number of elements as
%       strings in v_sel, indicating with a 1 if the cummulative impulse
%       response for that variable is required and a 0 if not.
%  plot_options.v_div_ss: A line vector with the same number of elements as
%       strings in v_sel, indicating with a 1 if that variable should be
%       divided by its steady state value and a 0 if not.
%  plot_options.horizon: a numebr indicating the horizon of teh responses to
%       be included in the graphs.
%  plot_options.n_col: a number indicating the number of columns in each the graph.
%  plot_options.n_row: a number indicating the number of rows in each the graph.
%  plot_options.grid: equal to 1 if a grid in each graph is desired, 0 otherwise.
%  plot_options.latex: equal to 1 if a latex format is wanted for titles, 0
%       otherwise. If 1 is selected, the next two have to be included.
%  plot_options.u_name (OPTIONAL): A set of strings of the dimension as u_sel 
%       with the LaTex name for each shock to be included in the title of each graph. 
%       If not included, the name in u_sel is used in the title.
%  plot_options.v_name (OPTIONAL): A set of strings of the dimension as v_sel 
%       with the LaTex name for each variable to be included in the title of each graph. 
%       If not included, the name in v_sel is used in the title.
%
% The function uses the sub-funtion loc defined at the end.

tt=plot_options.horizon;
ir_all=nan(tt,length(plot_options.v_sel));

figure('NumberTitle', 'off','Name',char(plot_options.u_title(1)),'units','normalized','outerposition',[0 0 1 1])
set(gcf,'WindowStyle','docked');


for jj=1:length(plot_options.v_sel)
    subplot(plot_options.n_row,plot_options.n_col,jj);
    va = plot_options.v_sel{jj};
    for ii=1:length(plot_options.model_names)
        eval(['M_=data.M_.' plot_options.model_names{ii} ';'])
        eval(['oo_=data.oo_.' plot_options.model_names{ii} ';'])
        sh = plot_options.u_sel{ii};
        if isempty(loc(M_.endo_names,va))==0
            if plot_options.v_div_ss(jj)
                div_ss=abs(oo_.steady_state(loc(M_.endo_names,va)));
            else
                div_ss=1;
            end
            if isfield(oo_.irfs,[va '_eps_' sh])
                if plot_options.v_do_cusum(jj)==0
                    ir = getfield(oo_.irfs,[va '_eps_' sh])...
                        *plot_options.v_adj(jj)/div_ss;
                elseif plot_options.v_do_cusum(jj)==1
                    ir = cumsum(getfield(oo_.irfs,[va '_eps_' sh]))...
                        *plot_options.v_adj(jj)/div_ss;
                end
            else
                ir = zeros(1,plot_options.horizon);
            end
            ir_all(:,ii)=ir(1:tt);
        else
            ir_all(:,ii)=0;
        end
    end
    for ii=1:length(plot_options.model_names)
        hold on;
        plot(1:tt,ir_all(:,ii),[plot_options.marks{ii}],'LineWidth',1.25);
        
    end
    if plot_options.grid==1
        grid on;
    end
    if plot_options.latex==1
        va_name=plot_options.v_name{jj};
        sh_name=plot_options.u_name{1};
    else
        va_name=va;
        sh_name=sh;
    end
    title(['$' sh_name ' \Rightarrow ' va_name '$'],'interpreter','latex','FontSize',12);
    set(gca,'FontSize',8,'FontName','Times')
    xlim([1 tt]);
    yl=ylim;
    
    lim_=0.01;
    
    if max(abs(yl))<lim_
    ylim([min(-lim_,yl(1)) ,   max(lim_,yl(2))])
    end
    
    hold off;

end

