% GKfigure1.m -- Replicates Figure 5 in "Macroeconomic Policy Games" by
% Martin Bodenstein, Luca Guerrieri, and Joe LaBriola
%%
% All the .mod files used in this program were produced (using our toolbox)
% from the file rbcb_monprud.mod which contains the conditions for a
% decentralized equilibrium of the original model.
%
% The file rbcb_monprud_nobias.mod is a slight variant of the
% rbcb_monprud.mod where the biases in the objectives of the regulators
% relative to the utility of the representative agent have been excluded
%
% This program uses three additional .mod files
% 1) rbcb_monprud_nobias_ramsey_infl_bt -- produced from 
%    convertmodfiles('rbcb_monprud_nobias','nash','infl','bt')
% 2) rbcb_monprud_ramsey_infl_bt -- produced from 
%    convertmodfiles('rbcb_monprud','nash','infl','bt')
% 3) rbcb_monprud_nobias_nash_infl_bt -- produced from 
%    convertmodfiles('rbcb_monprud','nash','infl','bt')
%
% notice that to facilitate welfare comparisons across models, after
% processing, we added two variables to rbcb_monprud_ramsey_infl_bt.mod and
% to rbcb_monprud_nash_infl_bt.mod -- these variables are util_nobias and
% welf_nobias.
%
% When processing the Ramsey and Nash files anew, remember to set a value
% for nbeta, the discounnt factor of the planner(s).
% For instance, for rbcb_monprud_nobias_ramsey_infl_bt, the value of nbeta
% is set in rbcb_monprud_nobias_ramsey_infl_bt_paramfile.
% For this model replace the line 
% nbeta = 0.99;
% with the line
% nbeta = betap;
% where betap is the discount factor of the model's representative agent.
%
% Also remember to inspect the shock block of each .mod file. 
% Comparisons of policies will only be accurate if the same exogenous stochastic
% structure is assumed. The routine convertmodelfiles may have the side
% effect of commenting out the shock block of the .mod file -- it can be
% uncommented before running any simulations.
% 
% The convention used here is for the standard deviations of the
% innovations to multiply the innovations in the equations that define of the shocks
% processes. The shock block then needs to set std 1 to be consistent with
% this choice for the shocks that the user wishes to activate.
%
%
%%
clear all
setpathdynare6_5   % Make sure to update the paths in this .m file to reflect 
                 % the local installation of Dynare and of the toolbox that
                 % accompanies this paper.

% Take care of some housekeeping -- if they exist, delete 
% some files that are used to impose changes to the baseline
% parameter choices.
if exist('param_overwrite.mat','file')
    delete param_overwrite.mat
end

if exist('param_overwrite_extra.mat','file')
    delete param_overwrite_extra.mat
end

% define names of dynare model files
stem = 'rbcb_monprud';
nobias_ramsey_modnam = [stem,'_nobias_ramsey_infl_bt'];
ramsey_modnam = [stem,'_ramsey_infl_bt'];
nash_modnam = [stem,'_nash_infl_bt'];

modnam1 = nobias_ramsey_modnam;
modnam2 = ramsey_modnam;
modnam3 = nash_modnam;
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('Investment', 'Consumption ',...
                 'Price of Capital','Output',...
                 'Spread of Loan Rate \newline over Deposit Rate (AR)','Policy Interest Rate (AR)',...
                 'Inflation (AR)', 'Bank Tax');
             
pc = 'Percent Dev. from SS';
ppt = 'PPt. Dev. from SS';

ylabels = char(pc, pc,...
               pc, pc,...
               ppt, ppt,...
               ppt, pc);


%% set shock
shocknamevector = char('errz');       
shocknamevectortex = char('Technology Shock z_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),';'])
        order = 1;
        [zdata] = makeirfsecondorder(M_,oo_,nperiods,shocks,order);
        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_nobias_ramsey',num2str(shockcounter),...
           ' = 100*[in_irf_1_',num2str(shockcounter),'/in_ss_1_',num2str(shockcounter)...
           ', c_irf_1_',num2str(shockcounter),'/c_ss_1_',num2str(shockcounter)...
           ', q_irf_1_',num2str(shockcounter),...
           ', y_irf_1_',num2str(shockcounter),'/y_ss_1_',num2str(shockcounter)...           
           ', 4*spread_irf_1_',num2str(shockcounter),...
           ', 4*rmon_irf_1_',num2str(shockcounter),...
           ', 4*infl_irf_1_',num2str(shockcounter),...           
           ', bt_irf_1_',num2str(shockcounter),'];']);

    eval(['line_ramsey',num2str(shockcounter),...
           ' = 100*[in_irf_2_',num2str(shockcounter),'/in_ss_2_',num2str(shockcounter)...
           ', c_irf_2_',num2str(shockcounter),'/c_ss_2_',num2str(shockcounter)...
           ', q_irf_2_',num2str(shockcounter),...
           ', y_irf_2_',num2str(shockcounter),'/y_ss_2_',num2str(shockcounter)...           
           ', 4*spread_irf_2_',num2str(shockcounter),...
           ', 4*rmon_irf_2_',num2str(shockcounter),...
           ', 4*infl_irf_2_',num2str(shockcounter),...           
           ', bt_irf_2_',num2str(shockcounter),'];']);
       
    eval(['line_nash',num2str(shockcounter),...
           ' = 100*[in_irf_3_',num2str(shockcounter),'/in_ss_3_',num2str(shockcounter)...
           ', c_irf_3_',num2str(shockcounter),'/c_ss_3_',num2str(shockcounter)...
           ', q_irf_3_',num2str(shockcounter),...
           ', y_irf_3_',num2str(shockcounter),'/y_ss_3_',num2str(shockcounter)...           
           ', 4*spread_irf_3_',num2str(shockcounter),...
           ', 4*rmon_irf_3_',num2str(shockcounter),...
           ', 4*infl_irf_3_',num2str(shockcounter),...           
           ', bt_irf_3_',num2str(shockcounter),'];']);       
       
end

for shockcounter = 1:size(shocknamevector,1)
    shockname = deblank(shocknamevectortex(shockcounter,:));
    
    legendlist = cellstr(char('Cooperation, no bias',...
                              'Cooperation, biased objectives',...
                              'Open-loop Nash'));
                  
    figtitle = '';

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








