%% BBCDL_figure4.m produces Figure 4 of the paper 
%%
% The figure compare responses to a markup shock for the open-loop Nash 
% equilibrium for different instruments.
% 
% The files BBCDLmodelcomp_nash_c1pid_c2pid.mod 
% BBCDLmodelcomp_nash_c1dnxf_c2dnxf';
% BBCDLmodelcomp_nash_c1y_c2y';
% are created by invoking 
% convertmodelfiles("BBCDLmodelcomp","nash","instrument1","instrument2")
% where instrument1 and instrument2 are the desired instrument combination.
% 
% Remember to edit the parameters governing the objective of the planner.
% Those parameters sit in the files
% BBCDLmodelcomp_nash_c1pid_c2pid_paramfile.m
% BBCDLmodelcomp_nash_c1dnxf_c2dnxf_paramfile.m
% BBCDLmodelcomp_nash_c1y_c2y_paramfile.m
%
%
% Also remember to check out the shock block of each model file.
% The shock blocks are commented out by convertmodelfiles.
% They need to be uncommented to obtain the correct second-solution.
% Remember that the convention used here is that the innovations are
% defined as Standard Normal (std 1) in the shock block, but are multiplied
% by their standard deviations in the equations that define the shock
% processes in the .mod files.
% NB: This convention is not necessary for the toolbox to work correctly.
%%

clear all
setpathdynare6_5

global paramlabels_ paramvec_ M_ oo_

stem = 'BBCDLmodelcomp';

modnam1 = [stem,'_nash_c1pid_c2pid'];
modnam2 = [stem,'_nash_c1dnxf_c2dnxf'];
modnam3 = [stem,'_nash_c1y_c2y'];
nummod = 3;

% solve models using dynare
for counter = 1:nummod    
    eval(['dynare ',eval(['modnam',num2str(counter)]),' noclearall nostrict'])
    eval(['M',num2str(counter),'_ = M_;'])
    eval(['oo',num2str(counter),'_ = oo_;'])
    eval(['ys',num2str(counter),'= oo_.dr.ys;'])
    eval(['nvars',num2str(counter),'= size(M_.endo_names,1);'])
    eval(['nshocks',num2str(counter),'= size(M_.exo_names,1);'])
end

%% generate model IRFs
nperiods = 10;

titlelist = char('Output, home','Output, foreign',...
                 'Output price inflation, home','Output price inflation, foreign',...
                 'Relative consumption','Real exchange rate');

pc = 'Percent Dev. from SS';
ppt = 'PPt. Dev. from SS';             
             
ylabels = char(pc, pc,...
               ppt, ppt,...
               pc, pc);


%set shock
shocknamevector =  char('c1taupvarerr');             
shocknamevectortex =  char('Technology Shock z_t',...
                           'Cost Push Shock \tau_{p,t}');  

for shockcounter =1:size(shocknamevector,1)
    shockname = deblank(shocknamevector(shockcounter,:));
    shocksize = 1;
    
    for counter = 1:nummod
        nshocks = eval(['nshocks',num2str(counter)]);
        shocks  = zeros(1,nshocks);
        shocks(1,strmatch(shockname,eval(['M',num2str(counter),'_.exo_names']),'exact')) = shocksize;
        eval(['shocks',num2str(counter),'= shocks;']);
    end
    
    for counter = 1:nummod
        eval(['M_ = M',num2str(counter),'_;'])
        eval(['oo_ = oo',num2str(counter),'_;'])
        eval(['nvars = nvars',num2str(counter),';'])
        eval(['shocks = shocks',num2str(counter),';'])
%        [zdata] = makeirfsecondorder(M_,oo_,nperiods,shocks,2);
        [zdata] = makeirfsecondorder(M_,oo_,nperiods,shocks,1);
        endo_names_reordered = M_.endo_names(oo_.dr.order_var,:);
        for indxi = 1:nvars
%            eval([deblank(M_.endo_names(indxi,:)),'_ss_',num2str(counter),'_',num2str(shockcounter),' = oo_.dr.ys(indxi);']);
%            eval([deblank(endo_names_reordered(indxi,:)),'_irf_',num2str(counter),'_',num2str(shockcounter),'=zdata(:,indxi);']);
eval([deblank(M_.endo_names{indxi}),'_ss_',num2str(counter),'_',num2str(shockcounter),'= zdata.',deblank(M_.endo_names{indxi}),'_ss;']);
eval([deblank(endo_names_reordered{indxi}),'_irf_',num2str(counter),'_',num2str(shockcounter),'= zdata.',deblank(endo_names_reordered{indxi}),'_irf;']);
        end
    end
    
    eval(['line_nash_1',num2str(shockcounter),...
           ' = 100*[c1y_irf_1_',num2str(shockcounter),...
           ', c2y_irf_1_',num2str(shockcounter),...   
           ', 4*c1pid_irf_1_',num2str(shockcounter),...
           ', 4*c2pid_irf_1_',num2str(shockcounter),...           
           ', c1c_irf_1_',num2str(shockcounter),'-c2c_irf_1_',num2str(shockcounter),...
           ', c1rer_irf_1_',num2str(shockcounter),'];']);

    if nummod > 1
    eval(['line_nash_2',num2str(shockcounter),...
           ' = 100*[c1y_irf_2_',num2str(shockcounter),...
           ', c2y_irf_2_',num2str(shockcounter),...
           ', 4*c1pid_irf_2_',num2str(shockcounter),...
           ', 4*c2pid_irf_2_',num2str(shockcounter),...           
           ', c1c_irf_2_',num2str(shockcounter),'-c2c_irf_2_',num2str(shockcounter),...
           ', c1rer_irf_2_',num2str(shockcounter),'];']);     
    end
       
    if nummod > 2
    eval(['line_nash_3',num2str(shockcounter),...
           ' = 100*[c1y_irf_3_',num2str(shockcounter),...
           ', c2y_irf_3_',num2str(shockcounter),...
           ', 4*c1pid_irf_3_',num2str(shockcounter),...
           ', 4*c2pid_irf_3_',num2str(shockcounter),...           
           ', c1c_irf_3_',num2str(shockcounter),'-c2c_irf_3_',num2str(shockcounter),...
           ', c1rer_irf_3_',num2str(shockcounter),'];']);     
    end   
end;

%%
for shockcounter = 1:size(shocknamevector,1)
    shockname = deblank(shocknamevectortex(shockcounter,:));
    
    legendlist = cellstr(char('Output price inflation',...
                              'Nom. ex. rate (growth)',...  
                              'Real output'));
    figtitle = [''];
    if nummod == 1
        eval(['figurehandle_shock',num2str(shockcounter),...
            ' = makechart(titlelist,legendlist,figtitle,ylabels,',...
            'line_nash_1',num2str(shockcounter),...
            ');']);
    elseif nummod == 2
        eval(['figurehandle_shock',num2str(shockcounter),...
            ' = makechart(titlelist,legendlist,figtitle,ylabels,',...
            'line_nash_1',num2str(shockcounter),',',...
            'line_nash_2',num2str(shockcounter),...
            ');']);
    elseif nummod == 3
        eval(['figurehandle_shock',num2str(shockcounter),...
            ' = makechart(titlelist,legendlist,figtitle,ylabels,',...
            'line_nash_1',num2str(shockcounter),',',...
            'line_nash_2',num2str(shockcounter),',',...
            'line_nash_3',num2str(shockcounter),...
            ');']);
    end
    
%     printpref
%     directory = pwd;
%     cd ../..
%     directory1 = cd;
%     cd(directory);
% 
%     print('-dpsc',[directory1,'/paper/figurespaper/BBCDLfigure4',deblank(shocknamevector(shockcounter,:))]); 
end;





