% This program replicates Figure 1 in "A Static Capital Buffer is 
% Hard to Beat" by Matthew Canzoneri, Behzad Diba, Luca Guerrieri 
% and Arsenii Mishin.
%
% The program produces Figure 1, corresponding to the empirical and model 
% relationships between the ratio of non-financial credit to GDP and GDP two
% years ahead using data for the United States from the first quarter
% of 1980 through the fourth quarter of 2024.

% Empirical validation of the model

%% Empirics

clear all
close all


setpath_windows

%global M_ oo_
global overwrite overwrite_param_names 



data_table = readtable('Commercial_banks.xlsx','Sheet',"quarterly");
data_table = data_table(year(data_table.observation_date)>1979 & year(data_table.observation_date)<2025,:);

data_table_end = readtable('Commercial_banks.xlsx','Sheet',"Quarterly, End of Quarter");
data_table_end = data_table_end(year(data_table_end.observation_date)>1979 & year(data_table_end.observation_date)<2025,:);

data_table.CRDQUSAPABIS = data_table_end.CRDQUSAPABIS;



% Selecting sample horizons:
% from 01-Jan-1947 to 01-Jul-2022
% Possible Months (quarterly data): Jan, Apr, Jul, Oct

n_obs = length(data_table.observation_date);

[gdp_cyclical,gdp_trend] = hpfast(log(data_table.A191RA3Q086SBEA), 1600);

% credit_real = data_table.total_credit./data_table.gdp_price_deflator;
% [credit_trend,credit_cyclical] = hpfilter(log(credit_real), 400000);

% Notice that in defining the credit-to-GDP ratio, we multiply the credit
% measure by 1000 to express it in $ billions, the same units as GDP.
credit2gdp = 1000*data_table.CRDQUSAPABIS./(4*data_table.NGDPSAXDCUSQ);

[credit2gdp_cyclical,credit2gdp_trend] = hpfast(credit2gdp, 400000);

%% Model
mod00 = 'model_safe_all_shocks';
mod11 = 'model_mixed_all_shocks';
mod10 = 'model_risky_all_shocks';
mod01 = 'model_spare_all_shocks';

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

maxiter = 20;


nperiods = round(n_obs*4/3); %including pre-sample that consists of 1/4 of observations
n_simul = 1000;
irf_shock = char('eps_a','eps_isp','eps_tau');

randn('seed',1);
nshocks = size(irf_shock,1);

f0_ytot_vec = zeros(n_simul,n_obs);
f0_ltot_vec = zeros(n_simul,n_obs);
f0_ltot2ytot_vec = zeros(n_simul,n_obs);

shockssequence = [0.0001,0.0001,0.0001];

[ ~, 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,...
    shockssequence,irf_shock,nperiods,0,maxiter);


for j=1:n_simul
j
shockssequence = randn(nperiods,nshocks);

[zdatapiecewise, zdatass, oo00_,  M00_, violvecbool_last_shock] = ...
  solve_two_constraints_add_violvecbool_no_parse(M00_,M10_,M01_,M11_,oo00_,...
    options00_,constraint1, constraint2,...
    constraint_relax1, constraint_relax2,...
    shockssequence,irf_shock,nperiods,0,maxiter);


ytot_pos = strmatch('ytot',M00_.endo_names,'exact');
ltot_pos = strmatch('ltot',M00_.endo_names,'exact');
ltot2ytot_pos = strmatch('ltot2ytot',M00_.endo_names,'exact');
            

f0_ytot_data = zdatapiecewise(nperiods-n_obs+1:end,ytot_pos);
f0_ltot_data = zdatapiecewise(nperiods-n_obs+1:end,ltot_pos);
f0_ltot2ytot_data = zdatapiecewise(nperiods-n_obs+1:end,ltot2ytot_pos);

f0_ytot_ss = zdatass(ytot_pos);
f0_ltot_ss = zdatass(ltot_pos);
f0_ltot2ytot_ss = zdatass(ltot2ytot_pos);
        
[f0_gdp_cyclical,f0_gdp_trend] = hpfast(f0_ytot_data/f0_ytot_ss, 1600);
[f0_ltot_cyclical,f0_ltot_trend] = hpfast(f0_ltot_data/f0_ltot_ss, 400000);
[f0_ltot2ytot_cyclical,f0_ltot2ytot_trend] = hpfast(f0_ltot2ytot_data/4, 400000);

f0_ytot_vec(j,:) = f0_gdp_cyclical';
f0_ltot_vec(j,:) = f0_ltot_cyclical';
f0_ltot2ytot_vec(j,:) = f0_ltot2ytot_cyclical';


end


%% Results

shift = 8;
CI_level = 0.90;

%% Relationship between Credit-to-GDP ratio and shift quarters ahead GDP
% Model

[r_lb, r_up, r_emp, xconf, yconf, credit_shifted_cyclical_emp, gdp_emp_predicted] = ...
    function_shift_credit_conf_int(f0_ltot2ytot_vec, f0_ytot_vec, credit2gdp_cyclical, ...
    gdp_cyclical, n_simul, shift, CI_level);
%%
figure
p = fill(100*xconf,100*yconf,'red','DisplayName',"90% Confidence Interval (Model-Generated)");
p.FaceColor = [1 0.8 0.8];      
p.EdgeColor = 'none';

hold on
scatter(100*credit2gdp_cyclical(1:end-8),100*gdp_cyclical(9:end),'DisplayName','Observed Data')

plot(100*credit_shifted_cyclical_emp,100*gdp_emp_predicted,'LineWidth',2, 'DisplayName',"Regression Line (Fitted to Observed Data)");
xlabel('Non-Financial Credit-to-GDP Ratio, Perc. Point Dev. From HP Trend ($\lambda = 400,000$)', 'interpreter', 'latex')
ylabel('Two-Year-Ahead Real GDP, Perc. Point Dev. From HP Trend ($\lambda = 1,600$)', 'interpreter', 'latex')
text(-5, -5.0, ['Empirical Correlation = ' num2str(r_emp,2)])
text(-5, -5.5, ['90% Confidence Interval of Model Correlation = [' num2str(r_lb,3) ', ' num2str(r_up,2) ']'])
hold off
legend('Location','southeast')
legend show
xlim([-10 15])

grid on

print -depsc2 figure_1_credit_gdp_compare_new.eps

