% GKfigure3.m -- third figure for the GK model
% Replicates Figure 7 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.
%
%%
clear

setpathdynare6_5;     % make sure to update to reflect the local installation 
                    % Dynare and of the toolbox that accompanies the paper

load_results = 0;   % set to 0 to generate the results anew; 
                    % set to 1 to load from disk


paramfilename = 'rbcb_monprud_paramfile';

nperiods = 2;  % set to 2 to start from the nonstochastic steady state
               % set to a large number i.e., 1000, to start from points 
               % sampled from the ergodic space of the model
               
nreps = 1;     % if sampling from the erogic space, set to a large number
               % i.e., 1000, to obtain results for an average of the points
               % drawn from the ergodic space

order = 2;
shocks_labels = char('errz');  % shocks active in the simulation
                               % make sure that the shocks here 
                               % match the shocks active in the shocks
                               % section of the Dynare model file.
                               % The convention used by this program is
                               % that the shock processes are defined in
                               % the .mod file so that the innovation to
                               % each process is multiplied by its standard 
                               % deviation -- for the programs to work
                               % the shock section needs to declare all
                               % innovations as Standard Normal 
                               % -- set std 1.
                               
s1 = RandStream('mt19937ar','seed',3);

nshocks = size(shocks_labels,1);
shock_mat = randn(s1,nreps,nshocks,nperiods);

if nperiods <3
    shock_mat = zeros(nreps,nshocks,nperiods);
end



if load_results
    load welfare_results_baseline
else
    
    eval(paramfilename)
    
    % Run model
    my_mpr_range =[0 0.05 0.1 0.15 .2 .25 0.5 .4 0.6 .8 1];  
    my_cp_range = [0 0.1 .2 .4 .8 1]; 
    
    tau2to1 = zeros(length(my_mpr_range),length(my_cp_range));
    tau3to1 = zeros(length(my_mpr_range),length(my_cp_range));
    
    
    if exist('param_overwrite.mat','file')
        delete param_overwrite.mat
    end
    
    % NB: this model is out of the loop that, further down, will 
    % recompute the model solution
    dynare rbcb_monprud_nobias_ramsey_infl_bt  noclearall nostrict
    f1M_ = M_;
    f1oo_ = oo_;
    
    dynare rbcb_monprud_ramsey_infl_bt noclearall nostrict       
    f2M_ = M_;
    f2oo_ = oo_;
            
    dynare rbcb_monprud_nash_infl_bt noclearall nostrict
    f3M_ = M_;
    f3oo_ = oo_;
    
    
    % calculate vector for mapping the positions of the endogenous variables of model 2 into
    % those of model 1
    vec_2to1_pos = zeros(f1M_.endo_nbr,1);
    for this_var = 1:f2M_.endo_nbr
    % find a mapping from the declaration order of the first model to the
    % second model
    vec_2to1_pos(this_var) = strmatch(deblank(f1M_.endo_names{this_var}),f2M_.endo_names,'exact');
    end
    
    % calculate vector for mapping the positions of the endogenous variables of model 3 into
    % those of model 2
    
    % In this case, one has to deal with the lagrange multipliers for the Nash problem 
    % -- Remember to make sure that the
    % declaration of the objective functions in the original .mod file happens
    % at the end of the file, right before the declaration of the exogenous
    % block of equations
    
    vec_2to3_pos = zeros(f2M_.endo_nbr,1);
    for this_var = 1:f3M_.endo_nbr
        % find a mapping from the declaration order of the first model to the
        % second model
        varname = deblank(f3M_.endo_names{this_var});
        
        % deal with the fact that the Nash model has 2 Welfcorrection
        % variables -- both are mapped into the same Welfcorrection 
        % from teh Ramsey problem
        if regexp(varname,'Welfcorrection')
            varname = 'Welfcorrection';
        end
        
        if regexp(varname,'mult')
            % if mult is in the name of the variable, assume that variable is a
            % multiplier
            vec_2to3_pos(this_var) = strmatch(varname(3:end),f2M_.endo_names,'exact');
        else
            vec_2to3_pos(this_var) = strmatch(varname,f2M_.endo_names,'exact');
        end
    end

    
    % now lookp through the values chosen for the bias parameters
    my_mpr_index = 0;
    for my_mpr_loop = my_mpr_range
        my_mpr_index = my_mpr_index+1;
        
        my_cp_index = 0;
        for my_cp_loop = my_cp_range
            my_cp_index=my_cp_index+1;
            
            alt1p = my_mpr_loop;
            alt2p = my_cp_loop;
            save param_overwrite alt1p alt2p
                      
            % recompute the decision rules with updated parameters
            rbcb_monprud_ramsey_infl_bt.driver
            f2M_ = M_;
            f2oo_ = oo_;
            
            rbcb_monprud_nash_infl_bt.driver       
            f3M_ = M_;
            f3oo_ = oo_;
            
            f1Welf = zeros(nreps,1);
            f2Welf = zeros(nreps,1);
            f3Welf = zeros(nreps,1);
            
            for this_rep = 1:nreps
                [this_rep alt1p alt2p]
                
                % select the shock from the shock matrix
                this_shock = reshape(shock_mat(this_rep,:,:),nshocks,nperiods);
                
                % sample from the ergodic distribution of the second model
                % -- in this case biased cooperative ramsey
                [s2 penultimatepoint1 penultimatepoint2] = mymkdata2(f2M_,f2oo_,nperiods,order,shocks_labels,this_shock);
                f2Welf(this_rep) = s2.Welf_nobias_irf(end)+s2.Welf_nobias_ss;
                
                % ask how much would I pay to switch to unbiased policies
                s1 = mymkdata2(f1M_,f1oo_,1,order,shocks_labels,this_shock(:,end),penultimatepoint1(vec_2to1_pos),penultimatepoint2(vec_2to1_pos));
                f1Welf(this_rep) = s1.Welf_irf(end)+s1.Welf_ss;
                
                % ask how much I would have to get paid to switch to NASH
                % policies with a bias
                s3 = mymkdata2(f3M_,f3oo_,1,order,shocks_labels,this_shock(:,end),penultimatepoint1(vec_2to3_pos),penultimatepoint2(vec_2to3_pos));
                f3Welf(this_rep) = s3.Welf_nobias_irf(end)+s1.Welf_nobias_ss;
                
                
            end
            tau2to1(my_mpr_index,my_cp_index) = 100*mean(exp((1-betap)*(f1Welf-f2Welf))-1);
            tau3to1(my_mpr_index,my_cp_index) = 100*mean(exp((1-betap)*(f1Welf-f3Welf))-1);
            
        end
    end
    
    save welfare_results_baseline my_mpr_range tau2to1 tau3to1

end


%% Make figure 
figure

subplot(2,1,1)

spline_range = 0:0.02:my_mpr_range(end);
spline_tau = spline(my_mpr_range,tau2to1(:,1),spline_range);
plot(spline_range,spline_tau,'k-','LineWidth',1.5); hold on
spline_tau = spline(my_mpr_range,tau2to1(:,2),spline_range);
plot(spline_range,spline_tau,'b--','LineWidth',1.5);
spline_tau = spline(my_mpr_range,tau2to1(:,3),spline_range);
plot(spline_range,spline_tau,'r-.','LineWidth',1.5);
spline_tau = spline(my_mpr_range,tau2to1(:,4),spline_range);
plot(spline_range,spline_tau,'g:','LineWidth',1.5);
spline_tau = spline(my_mpr_range,tau2to1(:,5),spline_range);
plot(spline_range,spline_tau,'c.','LineWidth',1.5);


xlim([0 spline_range(end)])
xlabel('Altruism parameter for the Central Bank, \omega_{cb}')
ylabel('Percent of Consumption')
title('Welfare Cost of Stabilization Biases with Cooperative Policies')

subplot(2,1,2)
spline_range = 0:0.02:my_mpr_range(end);
spline_tau = spline(my_mpr_range,tau3to1(:,1),spline_range);
plot(spline_range,spline_tau,'k-','LineWidth',1.5); hold on
spline_tau = spline(my_mpr_range,tau3to1(:,2),spline_range);
plot(spline_range,spline_tau,'b--','LineWidth',1.5);
spline_tau = spline(my_mpr_range,tau3to1(:,3),spline_range);
plot(spline_range,spline_tau,'r-.','LineWidth',1.5);
spline_tau = spline(my_mpr_range,tau3to1(:,4),spline_range);
plot(spline_range,spline_tau,'g:','LineWidth',1.5);
spline_tau = spline(my_mpr_range,tau3to1(:,5),spline_range);
plot(spline_range,spline_tau,'c.','LineWidth',1.5);


xlim([0 spline_range(end)])
xlabel('Altruism parameter for the Central Bank, \omega_{cb}')
ylabel('Percent of Consumption')
title('Welfare Cost of Stabilization Biases with Open-loop Nash Policies')

legend('No altruism for MPR \omega_{mpr} = 0','MPR Altruism  \omega_{mpr} = 0.1','MPR Altruism \omega_{mpr} = 0.2','MPR Altruism \omega_{mpr} = 0.4','MPR Altruism \omega_{mpr} = 0.8','MPR Altruism \omega_{mpr} = 1','location','best')
printpref