%close all;
clc
%clear all;


%pam = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi phi_b phi_Y];
%save data_pam pam;


beta = 0.9848;
delta = 0.025;
alpha = 0.35;
A = 1;
omega_bar = 0.25;
chi = 6.43;
ksi = 0.5;
eta = 0.164; 
sigma = 1/eta;
eps_min = 1-eta;
ppi = 0.005;
gamma = (1-delta)/(1-ppi);
B = 0.618;
rho_a = 0.9;  % persistence coefficient of the AR(1) process the TFP shock
rho_ppi = 0.9;  % persistence coefficient of the AR(1) process the asset-quality shock
rho_omega = 0.9;
rho_b = 0.9;
phi = 0;  % the elasticity of B_{t+1} to private liquidity (L_pri = qk*S, where S denotes trading volume)

stdA = 1/100;
%stdP = 0.01/100;
stdP = 0.042/100;
stdO = 3/100;
stdB = 2/100;
% stdO = 0.5/100;
% stdB = 1/100;
pphi = 0;


pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
save data_pam_policy pam_policy;



% RECURSIVE characterization

NN = 500; % number of grids

%dynare IRF_tfp_2_with_figure noclearall
%qk_ss = exp(oo_.steady_state(1));
qk_ss = 0.6677;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SCENARIO II: INCOMPLETE INFORMATION %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% THOSE ARE REGIONS FOR FIGURES 1 and 2
qk_lb = 0.93*qk_ss;
qk_ub = 1.1*qk_ss;

% THOSE ARE REGIONS FOR FIGURE 3
% qk_lb = 0.618;
% qk_ub = 0.7;

qk = linspace(qk_lb,qk_ub, NN);

eps2 = 1./qk;

theta = ((eps2./eps_min).^(-sigma)).*(1-delta)./(ppi+((eps2./eps_min).^(-sigma)).*(1-delta));

Phi_K = qk./theta;

eps1 = 1./Phi_K;

zeta = (eps1./eps_min).^(-sigma);

qb = beta.*(1+zeta./(sigma-1));

R = 1./qb;

rk = (1-beta.*(1-delta).*((1+(theta.^sigma).*zeta./(sigma-1)).*omega_bar+(1-omega_bar))).*Phi_K.*R-omega_bar.*ppi.*qk;

B_K = delta./(zeta.^(1-1./sigma))-(rk+omega_bar.*ppi.*qk)-(1-delta).*omega_bar.*qk.*(theta.^(sigma-1));

I_K = (B_K+rk+omega_bar.*ppi.*qk).*zeta+(1-delta).*omega_bar.*qk.*(theta.^sigma).*zeta;

Y_K = rk./alpha;

C_K = Y_K-I_K-(1-qb).*B_K;

N = (((1-alpha)./chi).*(Y_K./C_K)).^(1./(1+ksi));

K = ((alpha.*A./rk).^(1./(1-alpha))).*N;

omega = (ppi./(ppi+1-delta)+(1-delta).*(theta.^sigma)./(ppi+1-delta)).*omega_bar;

Y = Y_K.*K;

I = I_K.*K;

C = C_K.*K;

BB = B_K.*K;

% welfare
V = (log(C)-chi.*(N.^(1+ksi))./(1+ksi))./(1-beta);

B_Y = B_K./Y_K;

I_Y = I_K./Y_K;

w = (1-alpha).*Y./N;

%Rk = theta.*(rk+(ppi+1-delta).*qk.*(omega+(1-omega)/theta))./qk;
%Rk = (rk+(1-delta).*qk)./qk;
%Rk = (rk+(ppi+1-delta).*(omega.*qk+(1-omega).*Phi_K))./qk;
Rk = R./theta;

S = (ppi+(1-delta).*(theta.^sigma).*zeta).*K;

spread = Rk-R;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TO SOLVE AND SAVE THE STEADY STATE  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[error,a] = min(abs(BB-B*ones(1,NN)));

qk_ss = qk(a);
Phi_k_ss = Phi_K(a);
qb_ss = qb(a);
B_ss = BB(a);
C_ss = C(a);
I_ss = I(a);
K_ss = K(a);
Y_ss = Y(a);
N_ss = N(a);
R_ss = R(a);
rk_ss = rk(a);
zeta_ss = zeta(a);
omega_ss = omega(a);
eps1_ss = eps1(a);
eps2_ss = eps2(a);
theta_ss = theta(a);
BY_ss = B_Y(a);
IY_ss = I_Y(a);
w_ss = w(a);
Lambda_ss = 1/C_ss;
Rk_ss = Rk(a);
spread_ss = spread(a);

% as: as denotes for adverse selection
steady_state_as_policy = [qk_ss Phi_k_ss eps1_ss eps2_ss theta_ss qb_ss R_ss rk_ss K_ss N_ss Y_ss w_ss I_ss C_ss omega_ss BY_ss Rk_ss];
save data_steady_state_as_policy steady_state_as_policy;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Taylor Rule: TFP Shock %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% stdA = 1/100;
% stdP = 0.01/100;
% stdO = 0.5/100;
% stdB = 1/100;

% stdA = 1/100;
% stdP = 0.042/100;
% stdO = 30/100;
% stdB = 2/100;


stdA = 1/100;
stdP = 0;
stdO = 0;
stdB = 0;

phi = 0;



pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
save data_pam_policy pam_policy;

dynare IRF_as_policy noclearall   
u_A0 = oo_.mean(end);    % Household Utility (V in dynare code)

a = 0;
% = -0.22 if response to L rather than L(-1)
b = 0.3;
nn = 40;

u_A = zeros(1,nn);
phi_index = zeros(1,nn);

for i=1:nn;
    
   phi = a+(b-a)/(nn-1)*(i-1);
   pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
   save data_pam_policy pam_policy;
   dynare IRF_as_policy noclearall 
   u_A(i)= oo_.mean(end);
   phi_index(i)=a+(b-a)/(nn-1)*(i-1);
   
end

[Uopt,ind]=max(u_A);
disp('Initial Utility, Optimal Utility,  phi')
Opt1 = [u_A0 Uopt phi_index(ind)];
disp('Welfare Gain (Omega) in percentage')
gain_u_A = (exp((1-beta)*(Uopt-u_A0))-1)*100
phi_A0_optimal = phi_index(ind);


% subplot(2,2,1)
% plot(phi_index,u_A,'LineWidth',1.5);
% xlabel('$\phi$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Conventional Liquidity Policy: TFP Shock','FontSize',11);







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Taylor Rule: Uncertainty Shock %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% stdA = 1/100;
% stdP = 0.01/100;
% stdO = 0.5/100;
% stdB = 1/100;

% stdA = 1/100;
% stdP = 0.042/100;
% stdO = 30/100;
% stdB = 2/100;

stdA = 0;
stdP = 0.042/100;
stdO = 0;
stdB = 0;

phi = 0;

save data_steady_state_as_policy steady_state_as_policy;

pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
save data_pam_policy pam_policy;

dynare IRF_as_policy noclearall   
u_P0 = oo_.mean(end);    % Household Utility (V in dynare code)

a = 0;
% = -0.22 if response to L rather than L(-1)
b = 0.3;
nn = 40;

u_P = zeros(1,nn);
phi_index = zeros(1,nn);

for i=1:nn;
    
   phi = a+(b-a)/(nn-1)*(i-1);
   pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
   save data_pam_policy pam_policy;
   dynare IRF_as_policy noclearall 
   u_P(i)= oo_.mean(end);
   phi_index(i)=a+(b-a)/(nn-1)*(i-1);
   
end

[Uopt,ind]=max(u_P);
disp('Initial Utility, Optimal Utility,  phi')
Opt1 = [u_P0 Uopt phi_index(ind)]
disp('Welfare Gain (Omega) in percentage')
gain_u_P = (exp((1-beta)*(Uopt-u_P0))-1)*100
phi_P0_optimal = phi_index(ind);

% subplot(2,2,2)
% plot(phi_index,u_P,'LineWidth',1.5);
% xlabel('$\phi$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Conventional Liquidity Policy: Uncertainty Shock','FontSize',11);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Taylor Rule: Resaleability Shock %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% stdA = 1/100;
% stdP = 0.01/100;
% stdO = 0.5/100;
% stdB = 1/100;

% stdA = 1/100;
% stdP = 0.042/100;
% stdO = 30/100;
% stdB = 2/100;

stdA = 0;
stdP = 0;
stdO = 3/100;
stdB = 0;

phi = 0;

save data_steady_state_as_policy steady_state_as_policy;

pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
save data_pam_policy pam_policy;

dynare IRF_as_policy noclearall   
u_O0 = oo_.mean(end);    % Household Utility (V in dynare code)

a = 0;
% = -0.22 if response to L rather than L(-1)
b = 0.3;
nn = 40;

u_O = zeros(1,nn);
phi_index = zeros(1,nn);

for i=1:nn;
    
   phi = a+(b-a)/(nn-1)*(i-1);
   pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
   save data_pam_policy pam_policy;
   dynare IRF_as_policy noclearall 
   u_O(i)= oo_.mean(end);
   phi_index(i)=a+(b-a)/(nn-1)*(i-1);
   
end

[Uopt,ind]=max(u_O);
disp('Initial Utility, Optimal Utility,  phi')
Opt1 = [u_O0 Uopt phi_index(ind)];
disp('Welfare Gain (Omega) in percentage')
gain_u_O = (exp((1-beta)*(Uopt-u_O0))-1)*100
phi_O0_optimal = phi_index(ind);

% subplot(2,2,3)
% plot(phi_index,u_O,'LineWidth',1.5);
% xlabel('$\phi$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Conventional Liquidity Policy: Resaleability Shock','FontSize',11);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4. Taylor Rule: Fiscal Shock %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% stdA = 1/100;
% stdP = 0.01/100;
% stdO = 0.5/100;
% stdB = 1/100;

% stdA = 1/100;
% stdP = 0.042/100;
% stdO = 30/100;
% stdB = 2/100;

stdA = 0;
stdP = 0;
stdO = 0;
stdB = 2/100;

phi = 0;

save data_steady_state_as_policy steady_state_as_policy;

pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
save data_pam_policy pam_policy;

dynare IRF_as_policy noclearall   
u_B0 = oo_.mean(end);    % Household Utility (V in dynare code)

a = 0;
% = -0.22 if response to L rather than L(-1)
b = 0.3;
nn = 40;

u_B = zeros(1,nn);
phi_index = zeros(1,nn);

for i=1:nn;
    
   phi = a+(b-a)/(nn-1)*(i-1);
   pam_policy = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi stdA stdP stdO stdB pphi];
   save data_pam_policy pam_policy;
   dynare IRF_as_policy noclearall 
   u_B(i)= oo_.mean(end);
   phi_index(i)=a+(b-a)/(nn-1)*(i-1);
   
end

[Uopt,ind]=max(u_B);
disp('Initial Utility, Optimal Utility,  phi')
Opt1 = [u_B0 Uopt phi_index(ind)];
disp('Welfare Gain (Omega) in percentage')
gain_u_B = (exp((1-beta)*(Uopt-u_B0))-1)*100
phi_B0_optimal = phi_index(ind);

phi_B0 = phi_index;

% subplot(2,2,4)
% plot(phi_index,u_B,'LineWidth',1.5);
% xlabel('$\phi$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Conventional Liquidity Policy: Fiscal Shock','FontSize',11);

% 
% figure(1)
% subplot(2,2,1)
% plot(phi_index,u_A,'LineWidth',1.5);
% xlabel('$\phi$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Conventional Liquidity Policy: TFP Shock','FontSize',11);
% 
% subplot(2,2,2)
% plot(phi_index,u_P,'LineWidth',1.5);
% xlabel('$\phi$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Conventional Liquidity Policy: Uncertainty Shock','FontSize',11);
% 
% subplot(2,2,3)
% plot(phi_index,u_O,'LineWidth',1.5);
% xlabel('$\phi$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Conventional Liquidity Policy: Resaleability Shock','FontSize',11);
% 
% subplot(2,2,4)
% plot(phi_index,u_B,'LineWidth',1.5);
% xlabel('$\phi$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Conventional Liquidity Policy: Fiscal Shock','FontSize',11);


%print -depsc2 hello.eps



% gain_conventional = [gain_u_A gain_u_P gain_u_O gain_u_B]
% optimal_conventional = [pphi_A0_optimal pphi_P0_optimal pphi_O0_optimal pphi_B0_optimal]
% 
% series_conventional = [pphi_B0;u_A;u_P;u_O;u_B];
% 
% save data_series_conventional series_conventional;
% 
% save data_optimal_conventional optimal_conventional;
% 
% gain_and_phi_conventional = [gain_conventional;optimal_conventional];
% 
% save data_gain_and_phi_conventional gain_and_phi_conventional;


% rho = 0.9
% gain_conventional =
%     0.1781    0.3017    0.0675    0.0857
% optimal_conventional =
%    -0.1737   -0.1737   -0.0811   -0.1274



u_no_policy = [u_A0 u_P0 u_O0 u_B0];
save data_u_no_policy_other u_no_policy;


gain_conventional = [gain_u_A gain_u_P gain_u_O gain_u_B];

optimal_conventional = [phi_A0_optimal phi_P0_optimal phi_O0_optimal phi_B0_optimal];

series_conventional = [phi_B0;u_A;u_P;u_O;u_B];

save data_series_conventional_other series_conventional;

% save data_optimal_conventional_other optimal_conventional;
% 
% gain_and_phi_conventional = [gain_conventional;optimal_conventional];
% 
% save data_gain_and_phi_conventional_other gain_and_phi_conventional;
% 
% 
% gain_conventional_high_rho = [gain_u_A gain_u_P gain_u_O gain_u_B]
% optimal_conventional_high_rho = [phi_A0_optimal phi_P0_optimal phi_O0_optimal phi_B0_optimal]
% 
% 
% series_conventional_high_rho = [phi_index;u_A;u_P;u_O;u_B];
% save data_series_conventional_high_rho_other series_conventional_high_rho;
% 
% gain_and_phi_conventional_high_rho = [gain_conventional_high_rho;optimal_conventional_high_rho];
% 
% save data_gain_and_phi_conventional_high_rho gain_and_phi_conventional_high_rho;

% high rho: rho = 0.95
% gain_u_B =
%     0.0900
% gain_conventional_high_rho =
%     0.1900    0.4402    0.1591    0.0900
% optimal_conventional_high_rho =
%    -0.1968   -0.1737   -0.1042   -0.1274


load data_series_conventional_other;
phi_index = series_conventional(1,:);
u_A = series_conventional(2,:);
u_P = series_conventional(3,:);
u_O = series_conventional(4,:);
u_B = series_conventional(5,:);


% Plot with initial calibration
figure(1)
subplot(1,3,1)
plot(phi_index,-(u_A-u_A(1))./u_A(1)*100,'blue','LineWidth',1.5);
xlabel('$\phi$','Interpreter','latex','FontSize',11);
xlim([0 0.3]);
ylabel('$\Delta(\phi)$','Interpreter','latex','FontSize',11);
title('TFP Shock','FontSize',11);
grid on 

subplot(1,3,2)
plot(phi_index,-(u_P-u_P(1))./u_P(1)*100,'blue','LineWidth',1.5);
xlabel('$\phi$','Interpreter','latex','FontSize',11);
xlim([0 0.3]);
ylabel('$\Delta(\phi)$','Interpreter','latex','FontSize',11);
title('Safety Shock','FontSize',11);
grid on 

subplot(1,3,3)
plot(phi_index,-(u_O-u_O(1))./u_O(1)*100,'blue','LineWidth',1.5);
xlabel('$\phi$','Interpreter','latex','FontSize',11);
xlim([0 0.3]);
ylabel('$\Delta(\phi)$','Interpreter','latex','FontSize',11);
title('Liquidity Shock','FontSize',11);
grid on 

   

