clear

% start from the non-stochastic steady state -- same for all the policies
% define a shock matrix that is 10000 by 3000
% compute welfare based on utility stream for 3000 periods
% average over the 10000 repetitions


% can you choose any initial point?

%% Run model

nperiods = 2;
nreps = 1;

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

modname = 'shadow_BBCDLmodelcomp_ramsey_c1pid_c2pid';

grid_weight = [-.5 -.4 -0.2 -0.1 -0.05 -0.02 0.0001 0.02 0.05 0.1 0.2 0.4 0.5];

ngrid_weight = length(grid_weight);
nshocks = size(shocks_labels,1);

if nperiods >2 
    shock_mat = randn(s1,nreps,nshocks,nperiods);
else
    shock_mat = zeros(nreps,nshocks,nperiods);
end
eval(['dynare ',modname,' noclearall nostrict'])
f1M_ = M_;
f1oo_ = oo_;



%% find a mapping between the positions of the shadow variables and the
% variables in the ramsey portion of the model -- will be used to reset
% predetermined conditions


% first count the number of variables in the shadow part of the model
shadow_pos = strmatch('shadow_',f1M_.endo_names);
n_shadow = length(shadow_pos);

vec2_pos = zeros(n_shadow,1);
for this_var_indx = 1:n_shadow
    this_var = f1M_.endo_names{shadow_pos(this_var_indx)};
    this_var = this_var(8:end);
    vec2_pos(this_var_indx) = find(strcmp(this_var,f1M_.endo_names));
end


f1Welf = zeros(nreps,ngrid_weight);
f2Welf = zeros(nreps,ngrid_weight);
f2Welfcon1 = zeros(nreps,ngrid_weight);

BBCDLmodelcomp_ramsey_c1pid_c2pid_paramfile

for this_weight=1:ngrid_weight
    shadow_c1rn_weight = grid_weight(this_weight);
    
    save extra_parameters shadow_c1rn_weight
    eval([modname,'.driver'])
    f1M_ = M_;
    f1oo_ = oo_;
    
    for this_rep = 1:nreps
        [this_rep shadow_c1rn_weight]
        
        this_shock = reshape(shock_mat(this_rep,:,:),nshocks,nperiods);
        [s1 penultimatepoint1 penultimatepoint2] = mymkdata2(f1M_,f1oo_,nperiods,order,shocks_labels,this_shock);
        f1Welf(this_rep,this_weight) = s1.Welf_irf(end)+s1.Welf_ss+...
                                       (s1.Welfcorrection_irf(end)+...
                                       s1.Welfcorrection_ss)/nbeta;
        
        penultimatepoint1_switch = penultimatepoint1;
        penultimatepoint1_switch(shadow_pos) = penultimatepoint1(vec2_pos);
        
        penultimatepoint2_switch = penultimatepoint2;
        penultimatepoint2_switch(shadow_pos) = penultimatepoint2(vec2_pos);
        
        
        s2 = mymkdata2(f1M_,f1oo_,1,order,shocks_labels,this_shock(:,end),penultimatepoint1_switch,penultimatepoint2_switch);
        f2Welf(this_rep,this_weight) = s2.shadow_Welf_irf(end)+...
                                       s2.shadow_Welf_ss+...
                                       (s2.shadow_Welfcorrection_irf(end)+...
                                       s2.shadow_Welfcorrection_ss)/nbeta;
               
    end
    
    
end


figure

Welf_diff = f1Welf-f2Welf;
plot(grid_weight,Welf_diff)

save shadow_ramsey_results_constinfl Welf_diff grid_weight
