% 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

plot_irf_switch = 0;    % switch to plot IRFs
overwrite_results = 0;  % set to 1 to overwrite the results used for the paper

load_from_disk = 0;     % set to 1 to load initial guess from disk
% Here the values of initial guesses match the values from the disk
% Commented lines show the initial values that were used to derive the
% results

make_latex_table = 1;   % set to 1 to export latex table with all rules for the paper

%% run reference model
nperiods = 5000;
maxiter = 50;

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

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

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

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

std_gamma_var_opt = std(100*gamma_var_opt);

welf_opt = mean(zdataopt(:,strcmp('welf',Mopt_.endo_names))+...
    zdatass(strcmp('welf',Mopt_.endo_names)));
welfc_opt = mean(zdataopt(:,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');
init_theta = 0;

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 = 1;
[~, zdatapiecewise, zdatass, oo00_,  M00_, violvecbool,...
        M10_,M01_,M11_,oo10_,oo01_,oo11_,...
        options00_,options10_,options01_,options11_] = ...
        solve_two_constraints_add_violvecbool('model_safe_all_shocks',mod10,mod01,mod11,...
        constraint1, constraint2,...
        constraint_relax1, constraint_relax2,...
        shockssequence,irfshock,nperiods,0,maxiter);

% 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';

nrules = 1; % number of rules

%initialize structure that will hold the estimation results
coefs_struct{nrules} = {};
coefs_name_struct{nrules} = {};
rule_label_struct{nrules} = {};
row_labels_struc{nrules} = {};

for set_rule = 1:nrules
    
    if set_rule == 1

        result_file_name = 'credit2gdp';
        rule_label = 'Credit-to-GDP Ratio Gap';
        rule_param_names = char('gammap_increment','ltot2ytot_gammap');

        param_lb = [-.01];

        param_ub = [0.2];

    end


    %% This part does some setup
    if ~load_from_disk
        if set_rule == 1

            ltot2ytot = zdataopt(:,strcmp('ltot2ytot',Mopt_.endo_names));

            ltot2ytot_gammap = 0.25; %-0.00001;
            gammap_increment = 0.06850299835205079; %0.1;
            % for coefficient 0.0025: 0.01116333007812500; %0.1; %0.01;

            init_theta = [gammap_increment ltot2ytot_gammap];
        
        end
    end



    %% optimization
    if load_from_disk
        % load starting guess from disk
        eval(['load ',result_file_name]);
        init_theta = distance_param;
    end

    % test the 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);

    dist
    [model_momm,...
          mean_negative_equity_at_default,...
          average_failure_rate,...
          average_gdp_loss]=make_moments(M00_,zdatapiecewise,zdatass);  
            

    init_theta1 = init_theta(1);
    mydist = @(init_theta1) distance_function_rules2(init_theta1,...
        init_theta(2:end),...
        rule_param_names,...
        M00_,M10_,M01_,M11_,oo00_,...
        options00_,constraint1, constraint2,...
        constraint_relax1, constraint_relax2);

    % options for the numerical minimizer
    %optimization_options = optimset('display','iter','MaxIter',150,'MaxFunEvals',1e10,'TolFun',10e-5,'TolX',1e-6);
    optimization_options = optimset('display','iter','MaxIter',450,'MaxFunEvals',1e10,'TolFun',10e-5,'TolX',1e-5);

    options_patternsearch = optimoptions('patternsearch');
    options_patternsearch = optimoptions(options_patternsearch,'display','iter','MaxIterations',100,'MeshTolerance',1e-5,'StepTolerance',1e-5);


    % iterate through guesses -- set a number for the rep loop; fminsearch has
    % better convergence properties if we set a low maxfunevals (above) and
    % call it many times --- we have reset the loop to 1 iteration after
    % ascertaining convergence
    thetaunc1 = init_theta1;
    thetaunc = init_theta;
    if overwrite_results
        exit_flag = 0;
        rep = 0
        while exit_flag<1 
            rep = rep+1;
            
            [thetaunc1, fval, exit_flag] = patternsearch(mydist,init_theta1,[],[],[],[],param_lb,param_ub,[],options_patternsearch);
            
            init_theta=[thetaunc1 init_theta(2:end)];
            thetaunc = init_theta;
            distance_param = thetaunc;
            eval(['save ',result_file_name,' distance_param fval exit_flag rep'])
            
            
        end

        [dist, zdatapiecewise, zdatass] = ...
            distance_function_rules(init_theta,rule_param_names,...
            M00_,M10_,M01_,M11_,oo00_,...
            options00_,constraint1, constraint2,...
            constraint_relax1, constraint_relax2);


    end

    % record coefficients
    buffer_vec(set_rule) = 100*(thetaunc(strcmp(cellstr(rule_param_names),'gammap_increment')));
    coefs_sruct{set_rule} = thetaunc(~strcmp(cellstr(rule_param_names),'gammap_increment'));
    if length(thetaunc(~strcmp(cellstr(rule_param_names),'gammap_increment')))== 1
        coefs_vec(set_rule) = thetaunc(~strcmp(cellstr(rule_param_names),'gammap_increment'));
    else
        coefs_vec(set_rule) = nan;
    end
    
 

    coefs_names_struct{set_rule} = rule_param_names(~strcmp(cellstr(rule_param_names),'gammap_increment'),:);

    row_labels_struc{set_rule} = rule_label;

    % standard deviation of capital requirement
    gamma_var = zdatapiecewise(:,strcmp('gamma_var',M00_.endo_names));
    mean_gamma_var_vec(set_rule) = 100*(mean(gamma_var)+zdatass(strcmp('gamma_var',M00_.endo_names)));
    std_gamma_var_vec(set_rule) = std(100*gamma_var);

    % average GDP loss given excess risk
    excess_risk_switch = zdatapiecewise(:,strcmp('excess_risk_switch',M00_.endo_names));
    gdp = zdatapiecewise(:,strcmp('ytot',M00_.endo_names));
    gdp_ss = zdatass(strcmp('ytot',M00_.endo_names));
    gdp_cycle = 100*bpass(gdp/gdp_ss,6,32);
    average_gdp_loss = sum(gdp_cycle(excess_risk_switch>0.00001))/sum(excess_risk_switch>0.00001);
    
    % 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;

    % welfare
    welf = mean(zdatapiecewise(:,strcmp('welf',M00_.endo_names))+...
           zdatass(strcmp('welf',M00_.endo_names)));
    welf_vec(set_rule) = welf;

    betap= M00_.params(strcmp('betap',M00_.param_names));
    
    % consumption equivalent variation (in basis points)
    tau_eq = 100 * (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

%% export results to LaTex table 
if make_latex_table

% append results for the optimal rule
coefs_vec(end+1) = nan;
coefs_names_struct{end+1} = {};
coefs_struct{end+1} = {};
buffer_vec(end+1) = 0;
mean_gamma_var_vec(end+1) = 100*(mean(gamma_var_opt)+gamma_ss_opt);
std_gamma_var_vec(end+1) = std_gamma_var_opt;
average_failure_rate_vec(end+1) = 0;
welf_vec(end+1) = welf_opt;
tau_eq_vec(end+1) = 0;
row_labels_struc{end+1} ='Optimal Rule';


% make latex table
results_table = [coefs_vec',buffer_vec',std_gamma_var_vec',mean_gamma_var_vec',average_failure_rate_vec',tau_eq_vec',welf_vec'];

texfilename = 'simple_rules_results_fixed_coefficients';
column_labels = char('Slope','Buffer','Std.Dev. $\gamma$','Ave. $\gamma$','Ave. Fail. Rate','Cons. Equiv. Variation','Welfare');

sideways_switch = 1;
caption = 'Evaluating Simple Rules';
column_precision = [6 3 3 3 3 3 4];
row_labels = char(row_labels_struc);
maketable(results_table,texfilename,column_labels,row_labels,sideways_switch,caption,column_precision)



% % make latex table
% results_table = [coefs_vec',buffer_vec',average_failure_rate_vec',tau_eq_vec',welf_vec'];
% 
% texfilename = 'simple_rules_results_fixed_coefficients';
% column_labels = char('Slope','Buffer','Ave. Fail. Rate','Cons. Equiv. Variation','Welfare');
% 
% sideways_switch = 1;
% caption = 'Evaluating Simple Rules';
% column_precision = [6 3 3 3 4];
% row_labels = char(row_labels_struc);
% maketable(results_table,texfilename,column_labels,row_labels,sideways_switch,caption,column_precision)

end