%% BBCDL_figure1and2.m produces Figures 1 and 2 of the paper 
% "Macroeconomic Policy Games" by Martin Bodenstein, Luca Guerrieri, and
% Joe LaBriola.
%%
% The figures compare responses to technology and markup shocks for
% the cooperative and open-loop Nash equilibria of
% the model defined in BBCDLmodelcomp.mod for the decentralized case.
% 
% The files BBCDLmodelcomp_ramsey_c1pid_c2pid.mod and
% BBCDLmodelcomp_nash_c1pid_c2pid.mod can be produced invoking:
% convertmodelfiles("BBCDLmodelcomp","ramsey","c1pid","c2pid")
% and
% convertmodelfiles("BBCDLmodelcomp","nash","c1pid","c2pid")
% 
% Remember to edit the parameters governing the objective of the planner.
% Those parameters sit in the files
% BBCDLmodelcomp_ramsey_c1pid_c2pid_paramfile.m
% for the cooperative solution and in
% BBCDLmodelcomp_nash_c1pid_c2pid_paramfile.m
% for the open loop Nash game.
% 
% For this example, we replace the line that says
% nbeta = 0.99;
% with the line
% nbeta = c1beta;
% so that the planner has the same discount factor as the representative
% agent in country 1 (and in country 2).
%
% 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_ iteration

stem = 'BBCDLmodelcomp';

%convertmodfiles('BBCDLmodelcomp','ramsey','c1pid','c2pid','c1beta',1)

ramsey_modnam = [stem,'_ramsey_c1pid_c2pid'];
nash_modnam = [stem,'_nash_c1pid_c2pid'];

modnam1 = ramsey_modnam;
modnam2 = nash_modnam;
nummod = 2;

% solve models using dynare
for counter = 1:nummod
    iteration = counter;
    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,...
               ppt, ppt,...
               pc, pc);


%set shock
shocknamevector =  char('c1zerr','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),'= 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_ramsey',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),'];']);


    eval(['line_nash',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;

%%
for shockcounter = 1:size(shocknamevector,1)
    shockname = deblank(shocknamevectortex(shockcounter,:));
    
    legendlist = cellstr(char('Cooperative Policy',...
                              'Open-loop Nash'));
    figtitle = [''];

    eval(['figurehandle_shock',num2str(shockcounter),...
        ' = makechart(titlelist,legendlist,figtitle,ylabels,',...
        'line_ramsey',num2str(shockcounter),',',...
        'line_nash',num2str(shockcounter),...
        ');']);
    
    %printpref
end;








