% This program replicates tables A.2, A.3, and A.4 in Online Appendix of
% "A Static Capital Buffer is Hard to Beat" by Matthew Canzoneri,
% Behzad Diba, Luca Guerrieri and Arsenii Mishin.
%
% This program computes the results (parameter values of the shock processes,
% variance decomposition, and matching moments) using a SMM procedure under
% an alternative calibration that considers TFP and ISP shocks (there
% is no volatility shock compared with the main framework)


clear all

setpath_windows

global M_ oo_ options_
global overwrite overwrite_param_names

%% set some switches

load_from_disk = 0;     % set to 1 to load initial guess from disk

if ~load_from_disk
%% run reference model
nperiods = 5000;
maxiter = 20;

irfshock = char('eps_tau','eps_a','eps_isp');
randn('seed',1)
shockssequence = randn(nperiods,3);

% draw a data sample from the calibrated model with optimal rule
modopt = 'model_calibration';


[zdatalinear, zdatass, ooopt_, Mopt_ ] = ...
    solve_no_constraint(modopt,...
    shockssequence,irfshock,nperiods);



% record some statistics
gamma_var = 100*zdatalinear(:,find(strcmp('gamma_var',Mopt_.endo_names)));
sigmap_c_opt = Mopt_.params(strcmp('sigmap_c',Mopt_.param_names));

std_gamma_var_opt = std(gamma_var);

welf_opt = mean(zdatalinear(:,strcmp('welf',Mopt_.endo_names))+...
    zdatass(strcmp('welf',Mopt_.endo_names)));
welfc_opt = mean(zdatalinear(:,strcmp('welfc',Mopt_.endo_names))+...
    zdatass(strcmp('welfc',Mopt_.endo_names)));

gamma_ss_opt = zdatass(strcmp('gamma_var',Mopt_.endo_names));

% select baseline calibrated model with static buffer
% model setup
mod00 = 'model_safe_opt_all_shocks';
mod11 = 'model_mixed_all_shocks';
mod10 = 'model_risky_all_shocks';
mod01 = 'model_spare_all_shocks';

rule_param_names = char('gammap_increment');


constraint1 = 'chi2_upper_bar+chi2_upper_bar_ss<0';    % If the Lagrange multiplier on the
% no-shorting constraints for risky
% loans is slack, it means that I want
% to make risky loans.
% Accordingly, go into 11 or 10

constraint_relax1 = '((l+l_ss)<0)|(l_upper_bar<0)'; % If it turns out that the quantity of
% risky loans is negative,  00
% notice that l is set to 0 in the 10
% file, therefore the first condition
% is irrelevant when transition back
% from 10, but it might be relevant
% when transition back from 11.


constraint2 = '((chi2_lower_bar<0)&(l_upper_bar>0))'; % If the Lagrange multiplier on the
% no-shorting constraints for risky
% loans  and on safe loans are both slack, it means that I want
% to make risky and safe loans. Go to
% 01 or 11, but since if this
% condition is met constraint1 will
% also be active, we can only go into
% 11. Hence, 01 is spare.

constraint_relax2 ='((l+l_ss)<0)|(l_upper_bar<0)';  % return to 10 or 00, but
% notice that
% constraint_relax1 is identical,
% therefore, if either
% condition is violated
% we will go back to 00

% draw a data sample from the calibrated model with a baseline static
% buffer
shockval = 0;
[~, zdatapiecewise, zdatass, oo00_,  M00_, violvecbool,...
        M10_,M01_,M11_,oo10_,oo01_,oo11_,...
        options00_,options10_,options01_,options11_] = ...
        solve_two_constraints_add_violvecbool(mod00,mod10,mod01,mod11,...
        constraint1, constraint2,...
        constraint_relax1, constraint_relax2,...
        shockval,'eps_a',nperiods,0,maxiter);

%init_theta = 0.078000000000000-0.072148591648013;
init_theta = 0.005851408352;
% [dist, zdatapiecewise, zdatass] = ...
%         distance_function_rules(init_theta,rule_param_names,...
%         M00_,M10_,M01_,M11_,oo00_,...
%         options00_,constraint1, constraint2,...
%         constraint_relax1, constraint_relax2);
% 
% 
% % record some statistics
% welf_baseline = mean(zdatapiecewise(:,strcmp('welf',M00_.endo_names))+...
%     zdatass(strcmp('welf',M00_.endo_names)));
% welfc_baseline = mean(zdatapiecewise(:,strcmp('welfc',M00_.endo_names))+...
%     zdatass(strcmp('welfc',M00_.endo_names)));
% 
% non_defaulted = zdatapiecewise(:,strcmp('non_defaulted',M00_.endo_names));
% non_defaulted_ss = zdatass(strcmp('non_defaulted',M00_.endo_names));
% average_failure_rate_baseline = 400*(1 - mean(non_defaulted_ss+non_defaulted));
% 
% gamma_ss_baseline = zdatass(strcmp('gamma_var',Mopt_.endo_names));
% buffer_baseline = gamma_ss_baseline-gamma_ss_opt;
% 
betap= M00_.params(strcmp('betap',M00_.param_names));
sigmap_c = M00_.params(strcmp('sigmap_c',M00_.param_names));
% tau_eq_baseline = 10000 * (1 - (1 - (welf_opt-(welf_baseline))/ ...
%                          (welfc_opt+1/((1-betap)*(1-sigmap_c))))^(1/(1-sigmap_c)) );
% 
% %% now run the simple rules
% %mod00 = 'model_safe_opt_all_shocks';
% 
% %initialize structure that will hold the estimation results
% %opt = 0.001679897308350
buffer_grid = 0.004:0.0005:0.012;

%buffer_grid = 0.00704;
nrules = length(buffer_grid);

average_failure_rate_vec = nan(nrules,1);
welf_vec = nan(nrules,1);
welfc_vec = nan(nrules,1);
buffer_vec = nan(nrules,1);
tau_eq_vec = nan(nrules,1);
gamma_vec = nan(nrules,1);
row_labels_struc{nrules} = {};


for set_rule = 1:nrules
    result_file_name = 'my_buffer_grid';
    this_buffer = buffer_grid(set_rule);

    rule_label = ['Buffer ', num2str(this_buffer)];
    
 
    init_theta = this_buffer;

    % call distance function
    [dist, zdatapiecewise, zdatass] = ...
        distance_function_rules(init_theta,rule_param_names,...
        M00_,M10_,M01_,M11_,oo00_,...
        options00_,constraint1, constraint2,...
        constraint_relax1, constraint_relax2);
    
    
    % record coefficients
    buffer_vec(set_rule) = this_buffer;
    
    row_labels_struc{set_rule} = rule_label;

    % mean_deposit_change = 
    if ~isnan(zdatapiecewise)

    % annualized average failure rate
    non_defaulted = zdatapiecewise(:,strcmp('non_defaulted',M00_.endo_names));
    non_defaulted_ss = zdatass(strcmp('non_defaulted',M00_.endo_names));
    average_failure_rate = 400*(1 - mean(non_defaulted_ss+non_defaulted));
    average_failure_rate_vec(set_rule)=average_failure_rate;
    
    % capital requirement
    gamma_vec(set_rule) = zdatass(strcmp('gamma_var',M00_.endo_names));
        
    % welfare
    welf = mean(zdatapiecewise(:,strcmp('welf',M00_.endo_names))+...
           zdatass(strcmp('welf',M00_.endo_names)));
    welf_vec(set_rule) = welf;
    
    welfc = mean(zdatapiecewise(:,strcmp('welfc',M00_.endo_names))+...
           zdatass(strcmp('welfc',M00_.endo_names)));
    welfc_vec(set_rule) = welfc;

    betap= M00_.params(strcmp('betap',M00_.param_names));
    
    % consumption equivalent variation (in basis points)
    tau_eq = 10000 * (1 - (1 - (welf_opt-(welf))/ ...
                       (welfc_opt+1/((1-betap)*(1-sigmap_c))))^(1/(1-sigmap_c)) );
    tau_eq_vec(set_rule) = tau_eq;
    end
end

%% export results to LaTex table 
% add results for baseline calibrated buffer


% append results for the optimal rule
% buffer_vec(end+1) = 0;
% average_failure_rate_vec(end+1) = 0;
% welf_vec(end+1) = welf_opt;
% welfc_vec(end+1) = welfc_opt;
% tau_eq_vec(end+1) = 0;
% row_labels_struc{end+1} ='Optimal Rule';

save grid_buffers buffer_vec average_failure_rate_vec welfc_vec welf_vec tau_eq_vec row_labels_struc

else
    load grid_buffers.mat

end

figure
plot(buffer_vec,welf_vec)

%%
figure
subplot(4,1,1)
welf_vec = welf_vec(1:end-1);
buffer_vec = buffer_vec(1:end-1);
[~,max_pos] = max(welf_vec);
buffer_opt = buffer_vec(max_pos);

welf_vec = spline(buffer_vec(~isnan(welf_vec)),welf_vec(~isnan(welf_vec)), buffer_vec);
plot(10000*buffer_vec,welf_vec,LineWidth=1.5)
hold on; hhh=vline(10000*buffer_opt,'k--','Optimal buffer');
hhh.LineWidth=1.5;
title('Welfare')
ylabel('Level')
%xlabel('Capital buffer, basis points')

subplot(4,1,2)
welfc_vec = welfc_vec(1:end-1);
welfd_vec = welf_vec - welfc_vec;
welfd_vec = spline(buffer_vec(~isnan(welfd_vec)),welfd_vec(~isnan(welfd_vec)), buffer_vec);
plot(10000*buffer_vec,welfd_vec,LineWidth=1.5)
hold on; hhh=vline(10000*buffer_opt,'k--');
hhh.LineWidth=1.5;

title('Welfare Contribution of Deposits')
ylabel('Level')
%xlabel('Capital buffer, basis points')

subplot(4,1,3)
welfc_vec = spline(buffer_vec(~isnan(welfc_vec)),welfc_vec(~isnan(welfc_vec)), buffer_vec);
plot(10000*buffer_vec,welfc_vec,LineWidth=1.5)
hold on; hhh=vline(10000*buffer_opt,'k--');
hhh.LineWidth=1.5;

title('Welfare Contribution of Consumption')
ylabel('Level')
%xlabel('Capital buffer, basis points')


subplot(4,1,4)
average_failure_rate_vec = average_failure_rate_vec(1:end-1);
average_failure_rate_vec = spline(buffer_vec(~isnan(average_failure_rate_vec)),average_failure_rate_vec(~isnan(average_failure_rate_vec)), buffer_vec);
plot(10000*buffer_vec,average_failure_rate_vec,LineWidth=1.5)
hold on; hhh=vline(10000*buffer_opt,'k--');
hhh.LineWidth=1.5;

title('Average Bank Failure Rate (annualized)')
ylabel('Percent')
xlabel('Capital buffer, basis points')

printpref

print -depsc2 figure_buffer.eps
% make latex table --- Continue editing here
% results_table = [coefs_vec',buffer_vec',std_gamma_var_vec',average_failure_rate_vec',welf_vec', tau_eq_vec'];
% 
% texfilename = 'simple_rules_results';
% column_labels = char('Slope','Buffer','Std.Dev. \gamma','Ave. Fail. Rate','Welfare','Cons. Equiv. Variation');
% 
% sideways_switch = 0;
% caption = 'Evaluating Simple Rules';
% column_precision = [6 2 2 2 3 3];
% row_labels = char(row_labels_struc);
% maketable(results_table,texfilename,column_labels,row_labels,sideways_switch,caption,column_precision)