% This file computes the dynamics varying the response of tax and generates 
% 2-D figures.

close all;
clear all;
fclose('all');
clc;

%run once to set up model and workspace
dynare Aggregatestry.mod noclearall;

initial_condition_states = oo_.mean;
initial_policy_functions = oo_.dr;

%define grid
begin_tauB=-0.01;
end_tauB=0.01;
gap_tauB=0.001;
corr_B = begin_tauB:gap_tauB:end_tauB;

results=cell(length(corr_B),1);
  for j = 1:length(corr_B)
        set_param_value('corr_tauaB',corr_B(j));
        info=stoch_simul(var_list_);
        if info==0
            results{j,1}=oo_;
        else
            %model could not be solved
            results{j,1}=NaN;
        end
  end

results_mean = zeros(M_.endo_nbr,length(corr_B));
results_stddev = zeros(M_.endo_nbr,length(corr_B));
results_corr_Y = zeros(M_.endo_nbr,length(corr_B));
y_simult=zeros(M_.endo_nbr,2);
y_sim=zeros(M_.endo_nbr,length(corr_B));
shock_matrix = zeros(1,M_.exo_nbr); %create shock matrix with number of time periods in rows

  for j = 1:length(corr_B)
      if  isstruct(results{j,1})==0
          results_mean(:,j)=NaN;
          results_stddev(:,j) = NaN;
          results_corr_Y(:,j) = NaN;
          y_sim(:,j)=NaN;
      else
results_mean(:,j)=results{j,1}.mean;
results_stddev(:,j) = sqrt(diag(results{j,1}.var));
results_corr_Y(:,j) = results{j,1}.var(:,5)./(sqrt(diag(results{j,1}.var))*sqrt(diag(results{j,1}.var(5,5))));
      end           
  end
