%% GKfigure1.m produces Figures 4

clear all
setpathdynare6_5

addpath ../BBCDL_model_2ndorder


LQ_ramsey_modnam = 'LQBBCDL';
ramsey_modnam = 'BBCDLmodelcomp_ramsey_c1pid_c2pid';

modnam1 = LQ_ramsey_modnam;
modnam2 = ramsey_modnam;
nummod = 2;

% solve models using dynare
for counter = 1:nummod      
    if counter == 2
        cd ../BBCDL_model_2ndorder
    end        
    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);'])

    if counter == 2
       cd ../LQ_BBCDL_model 
    end
end

% generate model IRFs
nperiods = 10;

titlelist = cellstr(char(...
              'Output gap, home country',...
              'Output gap, foreign country',...
              'Inflation, home country',...
              'Inflation, foreign country',...
              'Relative consumption gap',...
              'Real exchange rate gap'));             

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

%set shock
shocknamevector =  char('c1uerr');             
shocknamevectortex =  char('Cost Push Shock');  

for shockcounter =1:size(shocknamevector,1)
    shockname = deblank(shocknamevector(shockcounter,:));
    shocksize = 1;
    
    for counter = 1:nummod
        if counter == 2
            shockname = char('c1taupvarerr');
        end        
        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_ramseyCDL',num2str(shockcounter),...
           ' = 100*[c1y_irf_1_',num2str(shockcounter),...
           ', c2y_irf_1_',num2str(shockcounter),...
           ', 4*c1pi_irf_1_',num2str(shockcounter),...
           ', 4*c2pi_irf_1_',num2str(shockcounter),...           
           ', c1relc_irf_1_',num2str(shockcounter),...
           ', c1q_irf_1_',num2str(shockcounter),'];']);

    eval(['line_ramsey',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('LQ solution',...
                              'Toolbox solution'));
    figtitle = ['Optimal Policy under Cooperation'];

    eval(['figurehandle_shock',num2str(shockcounter),...
        ' = makechart(titlelist,legendlist,figtitle,ylabels,',...
        'line_ramseyCDL',num2str(shockcounter),',',...
        'line_ramsey',num2str(shockcounter),...
        ');']);
    
%     directory = pwd;
%     cd ..
%     directory1 = cd;
%     cd(directory);
% 
%     print('-dpsc',[directory1,'/paper/figurespaper/LQtoolboxcomparison']); 
end;    
    
  
