%% BBCDL_table1 produces Table 1 of the paper "Macroeconomic Policy Games" by 
% Martin Bodenstein, Luca Guerrieri, and Joe LaBriola.
% results reported in the paper are stored in BBCDLwelfarecomparison_1.mat

clear
setpathdynare6_5;

stem = 'BBCDLmodelcomp';

%ramsey model -- responses independent of instruments
ramsey_modnam = [stem,'_ramsey_c1pid_c2pid'];

%instrument choices 
instr1 = char({'c1pid','c1dcore','c1y','c1ngdp','c1dnxf'});
instr2 = char({'c2pid','c2dcore','c2y','c2ngdp','c2dnxf'});

%% Run Ramsey model to obtain initial conditions
nperiods = 1000;
nreps = 10000;

order = 2;
shocks_labels = char('c1zerr','c1taupvarerr','c2zerr','c2taupvarerr');
s1 = RandStream('mt19937ar','seed',2);

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

eval(['dynare ',eval(['ramsey_modnam']),' noclearall nostrict'])
f1M_ = M_;
f1oo_ = oo_;

%% open all .mod files 
% for ind1 = 1:size(instr1,1)
%     for ind2 = 1:size(instr2,1)
%         instrcomb = ['_',deblank(instr1(ind1,:)),'_',deblank(instr2(ind2,:))];
%         nash_modnam = [stem,'_nash',instrcomb];
%         open([nash_modnam,'.mod'])
%     end
% end

%% Run Nash models 
for ind1 = 1:size(instr1,1)
    for ind2 = 1:size(instr2,1)
        instrcomb = ['_',deblank(instr1(ind1,:)),'_',deblank(instr2(ind2,:))];
        nash_modnam = [stem,'_nash',instrcomb];
        data_nash_modnam = ['data_',nash_modnam];

%         if exist([data_nash_modnam,'.mat']) == 1
%             load([data_nash_modnam,'.mat'])
%             eval(['f',instrcomb,'_M_ = M_;']);
%             eval(['f',instrcomb,'_oo_ = oo_;']);
%             eval(['vec',instrcomb,'_pos = vec2_pos;']);
%         else   
            eval(['dynare ',eval(['nash_modnam']),' noclearall nostrict']);
            eval(['f',instrcomb,'_M_ = M_;']);
            eval(['f',instrcomb,'_oo_ = oo_;']);
            vec2_pos = zeros(M_.endo_nbr,1);
            for this_var = 1:M_.endo_nbr
                % find a mapping from the declaration order of the first model to the
                % second model
                varname = deblank(char(M_.endo_names(this_var,:)));
                if regexp(varname,'mult')
                    % if mult is in the name of the variable, assume that variable is a
                    % multiplier
                    vec2_pos(this_var) = strmatch(varname(3:end),f1M_.endo_names,'exact');
                else
                    vec2_pos(this_var) = strmatch(varname,f1M_.endo_names,'exact');
                end
            end        
            eval(['vec',instrcomb,'_pos = vec2_pos;']);
            save([data_nash_modnam],'instrcomb','M_','oo_','vec2_pos');        
%         end
    end
end

% now deal with the lagrange multipliers -- 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


%%
f1Welf = zeros(nreps,1);
f2Welf = zeros(nreps,1);
f2Welf1 = zeros(nreps,1);
f2Welf2 = zeros(nreps,1);
f2Welfcon1 = zeros(nreps,1);
for this_rep = 1:nreps
    this_rep
    
    % Ramsey
    this_shock = reshape(shock_mat(this_rep,:,:),nshocks,nperiods);
    [s1 penultimatepoint1 penultimatepoint2] = mymkdata2(f1M_,f1oo_,nperiods,order,shocks_labels,this_shock);
    f1Welf(this_rep) = s1.Welf_irf(end)+s1.Welf_ss;
     
    % Nash cases
    for ind1 = 1:size(instr1,1)
        for ind2 = 1:size(instr2,1)
            instrcomb = ['_',deblank(instr1(ind1,:)),'_',deblank(instr2(ind2,:))];
            eval(['f2M_ = f',instrcomb,'_M_;']);
            eval(['f2oo_ = f',instrcomb,'_oo_;']);
            eval(['vec2_pos = vec',instrcomb,'_pos ;']);
            
            s2 = mymkdata2(f2M_,f2oo_,1,order,shocks_labels,this_shock(:,end),penultimatepoint1(vec2_pos),penultimatepoint2(vec2_pos));
            f2Welf(this_rep) = s2.Welf_irf(end)+s2.Welf_ss;
            f2Welf1(this_rep) = s2.Welf1_irf(end)+s2.Welf1_ss;
            f2Welf2(this_rep) = s2.Welf2_irf(end)+s2.Welf2_ss;
            f2Welfcon1(this_rep) = s2.Welfcon1_irf(end)+s2.Welfcon1_ss;
            
            eval(['f',instrcomb,'_Welf(',num2str(this_rep),',1) = f2Welf(',num2str(this_rep),',1);']);
            eval(['f',instrcomb,'_Welf1(',num2str(this_rep),',1) = f2Welf1(',num2str(this_rep),',1);']);
            eval(['f',instrcomb,'_Welf2(',num2str(this_rep),',1) = f2Welf2(',num2str(this_rep),',1);']);
        end
    end
end

eval([eval('ramsey_modnam'),'_paramfile'])

%%
for ind1 = 1:size(instr1,1)
    for ind2 = 1:size(instr2,1)
        instrcomb = ['_',deblank(instr1(ind1,:)),'_',deblank(instr2(ind2,:))];
        
        eval(['f2Welf = f',instrcomb,'_Welf;']) 
        tau2to1 = 100*mean(((f1Welf-f2Welf)./omega_welf1./f2Welfcon1+1).^(1/(1-c1sigma))-1);

        tau2to1_matrix(ind1,ind2) = tau2to1;
        eval(['tau',instrcomb,'_to1 = tau2to1;']);    
        
        eval(['f2Welf1 = f',instrcomb,'_Welf1;'])
        Welf1_matrix(ind1,ind2) = mean(f2Welf1); 
        eval(['f2Welf2 = f',instrcomb,'_Welf2;'])
        Welf2_matrix(ind1,ind2) = mean(f2Welf2); 
        
    end
end

save('BBCDLwelfarecomparison','tau2to1_matrix','instr1','instr2','Welf1_matrix','Welf2_matrix')
