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

r_e_opt = zdataopt(:,strcmp('r_e',Mopt_.endo_names));
r_l_opt = zdataopt(:,strcmp('r_l',Mopt_.endo_names));
ytot_opt = zdataopt(:,strcmp('ytot',Mopt_.endo_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 = 7; % 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 = 7%1:nrules

    if set_rule == 1
        result_file_name = 'my_credit2gdp';
        rule_label = 'Credit-to-GDP Ratio Gap';
        rule_param_names = char('ltot2ytot_gammap','gammap_increment');

        param_lb = [-1
                    -.01];

        param_ub = [1
                    1];
    elseif set_rule == 2
        result_file_name = 'my_gdp';
        rule_label = 'GDP Gap';
        rule_param_names = char('ytot_gammap','gammap_increment');
        param_lb = [-1
                    -.01];

        param_ub = [1
                    1];
    elseif set_rule == 3
        result_file_name = 'my_expected_banking_spread';
        rule_label = 'Expected Banking Spread';
        rule_param_names = char('banking_spread_gammap','gammap_increment');
        param_lb = [-1
                    -.01];

        param_ub = [1
                    1];

    elseif set_rule == 4
        result_file_name = 'my_buffer';
        rule_label = 'Optimal Buffer';
        rule_param_names = char('gammap_increment');

        param_lb = [-.01];

        param_ub = [0.01];
    
    elseif set_rule == 5
        result_file_name = 'my_all_shocks';
        rule_label = 'All Shocks \& Lagged $\gamma$';
        rule_param_names = char('tau_gammap',...
            'eps_tau_gammap',...
            'a_gammap',...
            'eps_a_gammap',...
            'isp_gammap',...
            'eps_isp_gammap',...
            'rho_gammap',...
            'gammap_increment');

    elseif set_rule == 6
        result_file_name = 'my_all_state_variables_and_shocks';
        rule_label = 'All Shocks \& State Variables';
        rule_param_names = char('tau_gammap',...
            'eps_tau_gammap',...
            'a_gammap',...
            'eps_a_gammap',...
            'isp_gammap',...
            'eps_isp_gammap',...
            'k_gammap',...
            'q_gammap',...
            'rho_gammap',...
            'c_gammap',...
            'ginvest_gammap',...
            'l_gammap',...
            'gammap_increment');


    elseif set_rule == 7
        result_file_name = 'my_credit2gdp_no_optimization';
        rule_label = 'Credit-to-GDP Ratio Gap';
        rule_param_names = char('ltot2ytot_gammap','gammap_increment');

        param_lb = [-1.5
                    -1.5];

        param_ub = [1.5
                    1.5];
    end


    %% This part does some setup
    if ~load_from_disk
        if set_rule ==1 % rule that responds to the credit-to-GDP ratio
            ltot2ytot = zdataopt(:,strcmp('ltot2ytot',Mopt_.endo_names));

            olsb=regress(gamma_var_opt,[0*ltot2ytot+1 ltot2ytot]);
            %ltot2ytot_gammap = olsb(2);
            ltot2ytot_gammap = -0.000034989247322; %-0.00001;
            gammap_increment = 0.008424481600523; %0.01;

            init_theta = [ltot2ytot_gammap gammap_increment];

        elseif set_rule == 2

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

            olsb=regress(gamma_var_opt,[0*ytot+1 ytot]);
            ytot_gammap = 0.000985751819440; %olsb(2);

            gammap_increment = 0.007908363342285; %0.02;

            init_theta = [ytot_gammap gammap_increment];


        elseif set_rule == 3
            banking_spread = zdataopt(:,strcmp('banking_spread',Mopt_.endo_names));

            olsb=regress(gamma_var_opt,[0*banking_spread+1 banking_spread]);
            banking_spread_gammap = 0.051260653138161; %0.05;
            gammap_increment = 0.008474121093750; %0.01;

            init_theta = [banking_spread_gammap gammap_increment];

        elseif set_rule == 4 %simple buffer
            %gammap_increment = 0.0011;
            gammap_increment = 0.008427991999156; % buffer_baseline;
            init_theta = [gammap_increment];
            

        elseif set_rule == 5
            tau_var = zdataopt(:,strcmp('tau_var',Mopt_.endo_names));
            eps_tau = shockssequence(:,3);

            a = zdataopt(:,strcmp('a',Mopt_.endo_names));
            eps_a = shockssequence(:,1);

            isp = zdataopt(:,strcmp('isp',Mopt_.endo_names));
            eps_isp = shockssequence(:,2);

            [olsb,~,~,~,stats] = regress(gamma_var_opt,...
                [0*gamma_var_opt+1 ...
                tau_var eps_tau ...
                a eps_a ...
                isp eps_isp ...
                [0; gamma_var_opt(1:end-1)]]);

            tau_gammap = -0.747289505809556; %olsb(2);
            eps_tau_gammap = 0.004475424571238; %olsb(3);
            a_gammap = 0.000173779795356; %olsb(4);
            eps_a_gammap = -0.000025319227128; %olsb(5);
            isp_gammap = -0.023239800082801; %olsb(6);
            eps_isp_gammap = 0.000646118700751; %olsb(7);
            rho_gammap = 0.979846862688460; %olsb(8);
            gammap_increment = 0.000097073568982; %0.0005;

            init_theta = [tau_gammap eps_tau_gammap...
                a_gammap   eps_a_gammap...
                isp_gammap eps_isp_gammap ...
                rho_gammap ...
                gammap_increment];


        elseif set_rule == 6

            tau_var = zdataopt(:,strcmp('tau_var',Mopt_.endo_names));
            eps_tau = shockssequence(:,3);

            a = zdataopt(:,strcmp('a',Mopt_.endo_names));
            eps_a = shockssequence(:,1);

            isp = zdataopt(:,strcmp('isp',Mopt_.endo_names));
            eps_isp = shockssequence(:,2);

            k = zdataopt(:,strcmp('k',Mopt_.endo_names));
            q=zdataopt(:,strcmp('q',Mopt_.endo_names));
            d=zdataopt(:,strcmp('d',Mopt_.endo_names));

            [olsb,~,~,~,stats] = regress(gamma_var_opt,...
                [0*gamma_var_opt+1 ...
                tau_var eps_tau ...
                a eps_a ...
                isp eps_isp ...
                [0; k(1:end-1)] ...
                [0; q(1:end-1)]  ...
                [0; gamma_var_opt(1:end-1)]  ...
                ]);

            tau_gammap = -0.746518911631293; %-0.748106402481215; %olsb(2);
            eps_tau_gammap = 0.004470021885493; %0.004473277213358; %olsb(3);
            a_gammap = 0.000171146597192; %0.000167889801743; %olsb(4);
            eps_a_gammap = -0.000025784144348; %-0.000025748994180; %olsb(5);
            isp_gammap = -0.023232199767078; %-0.023210846538468; %olsb(6);
            eps_isp_gammap = 0.000644803875834; %0.000644072004898; %olsb(7);
            k_gammap = -0.000001106992801; %0; %olsb(8);
            q_gammap = 0.000001035546905; %0; %olsb(9);
            rho_gammap = 0.979812569490889; %0.980502858133236; %olsb(10);

            % The regression r2 is one indicating that the
            % no-switching rule can be expressed without including all of the state
            % variables.
            % Add the remaining state variables for the safe model to allow the
            % optimized rule to diverge from the no-switching rule if desirable

            c_gammap = 0.000004974320751; %0;
            ginvest_gammap = 0.000018326937731; %0;
            l_gammap = -0.000000445093424; %0;
            gammap_increment = 0.000097674227665; %0.000104610385389; %0.0002;

            init_theta = [tau_gammap eps_tau_gammap...
                a_gammap   eps_a_gammap...
                isp_gammap eps_isp_gammap ...
                k_gammap ...
                q_gammap ...
                rho_gammap ...
                c_gammap ...
                ginvest_gammap ...
                l_gammap ...
                gammap_increment];

            % init_theta=[-0.219125136233880
            %   0.003394210888858
            %  -0.000066667089555
            %  -0.000014653273946
            %   0.002646808122715
            %   0.000407132579206
            %   0.000017798523535
            %  -0.000377491417361
            %   0.795059633866177
            %   0
            %   0
            %   0
            %   0.000007791317646];


        elseif set_rule ==7 % rule that responds to the credit-to-GDP ratio
            ltot2ytot = zdataopt(:,strcmp('ltot2ytot',Mopt_.endo_names));

            olsb=regress(gamma_var_opt,[0*ltot2ytot+1 ltot2ytot]);
            ltot2ytot_gammap = 0.001;
            gammap_increment = buffer_baseline; %0.01;

            init_theta = [ltot2ytot_gammap gammap_increment];

        
        end
    end

    irfshock = char('eps_a');
    shockval = -0.5;
    nperiods = 25;
    % solve the model
    % eval(['dynare ',modnam,' noclearall']
    [ zdatalinear_, zdatapiecewise1_, 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,irfshock,nperiods,0,maxiter);

    [zdatapiecewise2_, zdatass_, ~,  ~, violvecbool_last_shock] = ...
        solve_two_constraints_add_violvecbool_no_parse(M00_,M10_,M01_,M11_,oo00_,...
        options00_,constraint1, constraint2,...
        constraint_relax1, constraint_relax2,...
        shockval,irfshock,nperiods,0,maxiter);

    if max(max(abs(zdatapiecewise1_-zdatapiecewise2_)))>10^(-15)
        error('Check no_parse')
    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);  
            


    mydist = @(init_theta) distance_function_rules(init_theta,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
    thetaunc = init_theta;
    if overwrite_results
        exit_flag = 0;
        rep = 0
        while exit_flag<1 
            rep = rep+1;
            
            if set_rule < 5
            [thetaunc, fval, exit_flag] = patternsearch(mydist,init_theta,[],[],[],[],param_lb,param_ub,[],options_patternsearch);
            
            init_theta=thetaunc;
            distance_param = thetaunc;
            eval(['save ',result_file_name,' distance_param fval exit_flag rep'])
            
            else
            [thetaunc,fval,exit_flag] = ...
                fminsearch('distance_function_rules',...
                init_theta,optimization_options,rule_param_names,...
                M00_,M10_,M01_,M11_,oo00_,...
                options00_,constraint1, constraint2,...
                constraint_relax1, constraint_relax2);
            init_theta = thetaunc;
            distance_param = thetaunc;
            eval(['save ',result_file_name,' distance_param fval exit_flag rep'])
            end
        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';
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 4 4 4 4 6 6];
row_labels = char(row_labels_struc);
maketable(results_table,texfilename,column_labels,row_labels,sideways_switch,caption,column_precision)

end