% solve for optimal policy rule
% calculate social welfare in the benchmark case and optimal rules, using second-order approximation 
% it will take hours to run these codes.

clear all

%% set the type of shocks

sig_ap=0;   % =1 or -1 implies that the economy faces a 1% positive or negative specific productivity shocks of POEs, =0 implies that the economy doesn't experience such a shock
sig_a=1;    % =1 or -1 implies that the economy faces a 1% positive or negative productivity shock, =0 implies that the economy doesn't experience this shock
sig_m=0;    % =1 or -1 implies that the economy faces a 1% positive or negative traditional monetary policy shock, =0 implies that the economy doesn't experience this shock
sig_xi=0;   % =1 or -1 implies that the economy faces a 1% positive or negative shadow bank financing efficiency shock, =0 implies that the economy doesn't experience this shock
sig_sigma=0;   % =1 or -1 implies that the economy faces a 1% positive or negative uncertainty shock, =0 implies that the economy doesn't experience such a shock
save shock_types sig_ap sig_a sig_m sig_xi sig_sigma

%% solve for optimal structural monetary policy 
vNU1=[-20:0.1:-10];
vNU2=[0];

for mi=1:1:length(vNU1)
    data1=vNU1(mi);
    for mm=1:1:length(vNU2)
        data2=vNU2(mm);
save info data1 data2
dynare benchmark.mod
totutility_mean=oo_.mean(41);
totutility_ss=oo_.steady_state(41);
var_ch=oo_.var(1,1);var_h=oo_.var(2,2);
ch_ss=oo_.steady_state(1);
vol_y(mi,mm)=sqrt(oo_.var(26,26))/oo_.steady_state(26);var_y(mi,mm)=oo_.var(26,26);var_pi(mi,mm)=oo_.var(8,8);
vol_pi(mi,mm)=sqrt(oo_.var(8,8))/oo_.steady_state(8);var_gdp(mi,mm)=oo_.var(48,48); var_loan(mi,mm)=oo_.var(57,57);
vol_pr(mi,mm)=sqrt(oo_.var(38,38))/oo_.steady_state(38);var_pr(mi,mm)=oo_.var(38,38);
vol_lls(mi,mm)=sqrt(oo_.var(44,44))/oo_.steady_state(44);var_lls(mi,mm)=oo_.var(44,44);
vol_l(mi,mm)=sqrt(oo_.var(57,57))/oo_.steady_state(57);
vol_yys(mi,mm)=sqrt(oo_.var(45,45))/oo_.steady_state(45);var_yys(mi,mm)=oo_.var(45,45);
W_gain(mi,mm)=(exp((totutility_mean-totutility_ss)*(1-BETA))-1)*10000;
end
end

[min_a,index]=max(W_gain,[],1);
opt_sm=vNU1(index);
save OPT_SM opt_sm

%% solve for optimal capital adequacy ratio
clear;
vNU1=[0];
vNU2=[5:0.1:10];
for mi=1:1:length(vNU1)
    data1=vNU1(mi);
    for mm=1:1:length(vNU2)
        data2=vNU2(mm);
save info data1 data2
dynare benchmark.mod
totutility_mean=oo_.mean(41);
totutility_ss=oo_.steady_state(41);
var_ch=oo_.var(1,1);var_h=oo_.var(2,2);
ch_ss=oo_.steady_state(1);
vol_y(mi,mm)=sqrt(oo_.var(26,26))/oo_.steady_state(26);var_y(mi,mm)=oo_.var(26,26);var_pi(mi,mm)=oo_.var(8,8);
vol_pi(mi,mm)=sqrt(oo_.var(8,8))/oo_.steady_state(8);var_gdp(mi,mm)=oo_.var(48,48); var_loan(mi,mm)=oo_.var(57,57);
vol_pr(mi,mm)=sqrt(oo_.var(38,38))/oo_.steady_state(38);var_pr(mi,mm)=oo_.var(38,38);
vol_lls(mi,mm)=sqrt(oo_.var(44,44))/oo_.steady_state(44);var_lls(mi,mm)=oo_.var(44,44);
vol_l(mi,mm)=sqrt(oo_.var(57,57))/oo_.steady_state(57);
vol_yys(mi,mm)=sqrt(oo_.var(45,45))/oo_.steady_state(45);var_yys(mi,mm)=oo_.var(45,45);
W_gain(mi,mm)=(exp((totutility_mean-totutility_ss)*(1-BETA))-1)*10000;
end
end

[min_a,index]=max(W_gain,[],2);
opt_car=vNU2(index);
save OPT_CAR opt_car

%% solve for optimal capital adequacy ratio and optimal structural monetary policy jointly
clear;
vNU1=[-15:0.1:-10];
vNU2=[-35:0.1:-33];
for mi=1:1:length(vNU1)
    data1=vNU1(mi);
    for mm=1:1:length(vNU2)
        data2=vNU2(mm);
save info data1 data2
dynare benchmark1.mod
totutility_mean=oo_.mean(41);
totutility_ss=oo_.steady_state(41);
var_ch=oo_.var(1,1);var_h=oo_.var(2,2);
ch_ss=oo_.steady_state(1);
vol_y(mi,mm)=sqrt(oo_.var(26,26))/oo_.steady_state(26);var_y(mi,mm)=oo_.var(26,26);var_pi(mi,mm)=oo_.var(8,8);
vol_pi(mi,mm)=sqrt(oo_.var(8,8))/oo_.steady_state(8);var_gdp(mi,mm)=oo_.var(48,48); var_loan(mi,mm)=oo_.var(57,57);
vol_pr(mi,mm)=sqrt(oo_.var(38,38))/oo_.steady_state(38);var_pr(mi,mm)=oo_.var(38,38);
vol_lls(mi,mm)=sqrt(oo_.var(44,44))/oo_.steady_state(44);var_lls(mi,mm)=oo_.var(44,44);
vol_l(mi,mm)=sqrt(oo_.var(57,57))/oo_.steady_state(57);
vol_yys(mi,mm)=sqrt(oo_.var(45,45))/oo_.steady_state(45);var_yys(mi,mm)=oo_.var(45,45);
W_gain(mi,mm)=(exp((totutility_mean-totutility_ss)*(1-BETA))-1)*10000;
end
end

maxValue = max(W_gain(:));
[indexi, indexj] = find(W_gain == maxValue);
opt_sm1=vNU1(indexi);
opt_car1=vNU2(indexj);
save OPT_JOINT opt_car1 opt_sm1


%% set the type of shocks

sig_ap=0;   % =1 or -1 implies that the economy faces a 1% positive or negative specific productivity shocks of POEs, =0 implies that the economy doesn't experience such a shock
sig_a=0;    % =1 or -1 implies that the economy faces a 1% positive or negative productivity shock, =0 implies that the economy doesn't experience this shock
sig_m=0;    % =1 or -1 implies that the economy faces a 1% positive or negative traditional monetary policy shock, =0 implies that the economy doesn't experience this shock
sig_xi=1;   % =1 or -1 implies that the economy faces a 1% positive or negative shadow bank financing efficiency shock, =0 implies that the economy doesn't experience this shock
sig_sigma=0;   % =1 or -1 implies that the economy faces a 1% positive or negative uncertainty shock, =0 implies that the economy doesn't experience such a shock
save shock_types sig_ap sig_a sig_m sig_xi sig_sigma

%% solve for optimal structural monetary policy 
clear all

vNU1=[-38.1:0.1:-37];
vNU2=[0];

for mi=1:1:length(vNU1)
    data1=vNU1(mi);
    for mm=1:1:length(vNU2)
        data2=vNU2(mm);
save info data1 data2
dynare benchmark.mod
totutility_mean=oo_.mean(41);
totutility_ss=oo_.steady_state(41);
var_ch=oo_.var(1,1);var_h=oo_.var(2,2);
ch_ss=oo_.steady_state(1);
vol_y(mi,mm)=sqrt(oo_.var(26,26))/oo_.steady_state(26);var_y(mi,mm)=oo_.var(26,26);var_pi(mi,mm)=oo_.var(8,8);
vol_pi(mi,mm)=sqrt(oo_.var(8,8))/oo_.steady_state(8);var_gdp(mi,mm)=oo_.var(48,48); var_loan(mi,mm)=oo_.var(57,57);
vol_pr(mi,mm)=sqrt(oo_.var(38,38))/oo_.steady_state(38);var_pr(mi,mm)=oo_.var(38,38);
vol_lls(mi,mm)=sqrt(oo_.var(44,44))/oo_.steady_state(44);var_lls(mi,mm)=oo_.var(44,44);
vol_l(mi,mm)=sqrt(oo_.var(57,57))/oo_.steady_state(57);
vol_yys(mi,mm)=sqrt(oo_.var(45,45))/oo_.steady_state(45);var_yys(mi,mm)=oo_.var(45,45);
W_gain(mi,mm)=(exp((totutility_mean-totutility_ss)*(1-BETA))-1)*10000;
end
end

[min_a,index]=max(W_gain,[],1);
opt_sm=vNU1(index);
save OPT_SM_sb opt_sm

%% solve for optimal capital adequacy ratio
clear;
vNU1=[0];
vNU2=[3:0.1:5.5];
for mi=1:1:length(vNU1)
    data1=vNU1(mi);
    for mm=1:1:length(vNU2)
        data2=vNU2(mm);
save info data1 data2
dynare benchmark.mod
totutility_mean=oo_.mean(41);
totutility_ss=oo_.steady_state(41);
var_ch=oo_.var(1,1);var_h=oo_.var(2,2);
ch_ss=oo_.steady_state(1);
vol_y(mi,mm)=sqrt(oo_.var(26,26))/oo_.steady_state(26);var_y(mi,mm)=oo_.var(26,26);var_pi(mi,mm)=oo_.var(8,8);
vol_pi(mi,mm)=sqrt(oo_.var(8,8))/oo_.steady_state(8);var_gdp(mi,mm)=oo_.var(48,48); var_loan(mi,mm)=oo_.var(57,57);
vol_pr(mi,mm)=sqrt(oo_.var(38,38))/oo_.steady_state(38);var_pr(mi,mm)=oo_.var(38,38);
vol_lls(mi,mm)=sqrt(oo_.var(44,44))/oo_.steady_state(44);var_lls(mi,mm)=oo_.var(44,44);
vol_l(mi,mm)=sqrt(oo_.var(57,57))/oo_.steady_state(57);
vol_yys(mi,mm)=sqrt(oo_.var(45,45))/oo_.steady_state(45);var_yys(mi,mm)=oo_.var(45,45);
W_gain(mi,mm)=(exp((totutility_mean-totutility_ss)*(1-BETA))-1)*10000;
end
end

[min_a,index]=max(W_gain,[],2);
opt_car=vNU2(index);
save OPT_CAR_sb opt_car

%% solve for optimal capital adequacy ratio and optimal structural monetary policy jointly
clear;
vNU1=[-55.5:.1:-55];
vNU2=[1.8:0.1:2];
for mi=1:1:length(vNU1)
    data1=vNU1(mi);
    for mm=1:1:length(vNU2)
        data2=vNU2(mm);
save info data1 data2
dynare benchmark.mod
totutility_mean=oo_.mean(41);
totutility_ss=oo_.steady_state(41);
var_ch=oo_.var(1,1);var_h=oo_.var(2,2);
ch_ss=oo_.steady_state(1);
vol_y(mi,mm)=sqrt(oo_.var(26,26))/oo_.steady_state(26);var_y(mi,mm)=oo_.var(26,26);var_pi(mi,mm)=oo_.var(8,8);
vol_pi(mi,mm)=sqrt(oo_.var(8,8))/oo_.steady_state(8);var_gdp(mi,mm)=oo_.var(48,48); var_loan(mi,mm)=oo_.var(57,57);
vol_pr(mi,mm)=sqrt(oo_.var(38,38))/oo_.steady_state(38);var_pr(mi,mm)=oo_.var(38,38);
vol_lls(mi,mm)=sqrt(oo_.var(44,44))/oo_.steady_state(44);var_lls(mi,mm)=oo_.var(44,44);
vol_l(mi,mm)=sqrt(oo_.var(57,57))/oo_.steady_state(57);
vol_yys(mi,mm)=sqrt(oo_.var(45,45))/oo_.steady_state(45);var_yys(mi,mm)=oo_.var(45,45);
W_gain(mi,mm)=(exp((totutility_mean-totutility_ss)*(1-BETA))-1)*10000;
end
end

maxValue = max(W_gain(:));
[indexi, indexj] = find(W_gain == maxValue);
opt_sm1=vNU1(indexi);
opt_car1=vNU2(indexj);
save OPT_JOINT_sb opt_car1 opt_sm1