close all;
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.831;
chi = 6.43;
ksi = 0.5;
eta = 0.1639; 
sigma = 1/eta;
eps_min = 1-eta;
ppi = 0.00319;
gamma = (1-delta)/(1-ppi);
B = 0.618;
rho_a = 0.8;  % 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.7;
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 = 0.01; %1/100;
% %stdP = 0.01/100;
% stdP = 0.042/100;
% stdO = 30/100; %30/100;
% stdB = 2/100; %2/100;
% % stdO = 0.5/100;
% % stdB = 1/100;
% pphi = 0;

%pampam_policy_unconventional = [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_pampam_policy_unconventional pampam_policy_unconventional;



% 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 = 1.2;

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.*omega_bar;

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. Unconventional Policy: TFP Shock %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

pphi = 0;

save data_steady_state_as_policy steady_state_as_policy;

pam_policy_unconventional = [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_unconventional pam_policy_unconventional;

dynare IRF_as_policy_unconventional_Wen 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.6;
nn =40;

% u_A = zeros(1,nn);
% pphi_index = zeros(1,nn);


for i=1:nn
    
   pphi = a+(b-a)/(nn-1)*(i-1);
   pam_policy_unconventional = [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_unconventional pam_policy_unconventional;
   dynare IRF_as_policy_unconventional_Wen noclearall 
   u_A(i)= oo_.mean(end);
   pphi_A0(i)=a+(b-a)/(nn-1)*(i-1);
   
end

[Uopt_A0,ind_A0]=max(u_A);
disp('Initial Utility, Optimal Utility,  phi')
Opt1_A0 = [u_A0 Uopt_A0 pphi_A0(ind_A0)];
disp('Welfare Gain (Omega) in percentage')
gain_u_A = (exp((1-beta)*(Uopt_A0-u_A0))-1)*100
pphi_A0_optimal = pphi_A0(ind_A0);




% figure(2)
% subplot(2,2,1)
% plot(pphi_A0,u_A,'LineWidth',1.5);
% xlabel('$\phi_\lambda$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Unconventional Policy: TFP Shock','FontSize',11);





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Unconventional Policy: Uncertainty Shock %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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



stdA = 0;
stdP = 0.042/100;
stdO = 0;
stdB = 0;
pphi = 0;
save data_steady_state_as_policy steady_state_as_policy;
pam_policy_unconventional = [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_unconventional pam_policy_unconventional;
dynare IRF_as_policy_unconventional_Wen 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.6;
nn = 40;

% u_P = zeros(1,nn);
% pphi_index = zeros(1,nn);


for i=1:nn;
    
   pphi = a+(b-a)/(nn-1)*(i-1);
   pam_policy_unconventional = [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_unconventional pam_policy_unconventional;
   dynare IRF_as_policy_unconventional_Wen noclearall 
   u_P(i)= oo_.mean(end);
   pphi_P0(i)=a+(b-a)/(nn-1)*(i-1);
   
end

[Uopt_P0,ind_P0]=max(u_P);
disp('Initial Utility, Optimal Utility,  phi')
Opt1_P0 = [u_P0 Uopt_P0 pphi_P0(ind_P0)];
disp('Welfare Gain (Omega) in percentage')
gain_u_P = (exp((1-beta)*(Uopt_P0-u_P0))-1)*100
pphi_P0_optimal = pphi_P0(ind_P0);


% subplot(2,2,2)
% plot(pphi_P0,u_P,'LineWidth',1.5);
% xlabel('$\phi_\lambda$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Unconventional Policy: Uncertainty Shock','FontSize',11);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Unconventional Policy: Resaleability Shock %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

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



pphi = 0;

save data_steady_state_as_policy steady_state_as_policy;

pam_policy_unconventional = [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_unconventional pam_policy_unconventional;

dynare IRF_as_policy_unconventional_Wen 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.6;
nn = 40;

% u_O = zeros(1,nn);
% pphi_index = zeros(1,nn);


for i=1:nn;
    
   pphi = a+(b-a)/(nn-1)*(i-1);
   pam_policy_unconventional = [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_unconventional pam_policy_unconventional;
   dynare IRF_as_policy_unconventional_Wen noclearall 
   u_O(i)= oo_.mean(end);
   pphi_O0(i)=a+(b-a)/(nn-1)*(i-1);
   
end


[Uopt_O0,ind_O0]=max(u_O);
disp('Initial Utility, Optimal Utility,  phi')
Opt1_O0 = [u_O0 Uopt_O0 pphi_O0(ind_O0)];
disp('Welfare Gain (Omega) in percentage')
gain_u_O = (exp((1-beta)*(Uopt_O0-u_O0))-1)*100
pphi_O0_optimal = pphi_O0(ind_O0);


pphi_B0 = pphi_O0;


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


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 4. Unconventional Policy: Fiscal Shock %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 
% % stdA = 1/100;
% % stdP = 0.01/100;
% % stdO = 0.5/100;
% % stdB = 1/100;
% 
% stdA = 0;
% stdP = 0;
% stdO = 0;
% stdB = 2/100; %2/100;
% 
% 
% pphi = 0;
% 
% save data_steady_state_as_policy steady_state_as_policy;
% 
% pam_policy_unconventional = [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_unconventional pam_policy_unconventional;
% 
% dynare IRF_as_policy_unconventional_Wen 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 = 1.2;
% nn = 3;
% 
% % u_B = zeros(1,nn);
% % pphi_index = zeros(1,nn);
% 
% 
% for i=1:nn;
%     
%    pphi = a+(b-a)/(nn-1)*(i-1);
%    pam_policy_unconventional = [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_unconventional pam_policy_unconventional;
%    dynare IRF_as_policy_unconventional_Wen noclearall 
%    u_B(i)= oo_.mean(end);
%    pphi_B0(i)=a+(b-a)/(nn-1)*(i-1);
%    
% end
% 
% [Uopt_B0,ind_B0]=max(u_B);
% disp('Initial Utility, Optimal Utility,  phi')
% Opt1_B0 = [u_B0 Uopt_B0 pphi_B0(ind_B0)];
% disp('Welfare Gain (Omega) in percentage')
% gain_u_B = (exp((1-beta)*(Uopt_B0-u_B0))-1)*100
% pphi_B0_optimal = pphi_B0(ind_B0);

% figure(1)
% subplot(1,3,1)
% plot(pphi_A0,u_A,'LineWidth',1.5);
% xlabel('$\phi_\lambda$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Unconventional Policy: TFP Shock','FontSize',11);
% subplot(1,3,2)
% plot(pphi_P0,u_P,'LineWidth',1.5);
% xlabel('$\phi_\lambda$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Unconventional Policy: Uncertainty Shock','FontSize',11);
% subplot(1,3,3)
% plot(pphi_O0,u_O,'LineWidth',1.5);
% xlabel('$\phi_\lambda$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Unconventional Policy: Resaleability Shock','FontSize',11);


%b = 0.78;

figure(2)
subplot(1,3,1)
plot(pphi_A0,-(u_A-u_A0)./u_A0.*100,'blue','LineWidth',1.5);
xlabel('$\phi$','Interpreter','latex','FontSize',20);
xlim([a 0.35]);
ylabel('$\Delta(\phi)$','Interpreter','latex','FontSize',20);
title('TFP Shock','Interpreter','latex','FontSize',20);
subplot(1,3,2)
plot(pphi_P0,-(u_P-u_P0)./u_P0.*100,'blue','LineWidth',1.5);
xlabel('$\phi$','Interpreter','latex','FontSize',20);
xlim([a b-0.1]);
ylabel('$\Delta(\phi)$','Interpreter','latex','FontSize',20);
title('Safety Shock','Interpreter','latex','FontSize',20);
subplot(1,3,3)
plot(pphi_O0,-(u_O-u_O0)./u_O0.*100,'blue','LineWidth',1.5);
xlabel('$\phi$','Interpreter','latex','FontSize',20);
xlim([a b-0.1]);
ylabel('$\Delta(\phi)$','Interpreter','latex','FontSize',20);
title('Liquidity Shock','Interpreter','latex','FontSize',20);

% print -depsc2 policy_osr_unconventional.eps
 
% subplot(2,2,4)
% plot(pphi_B0,u_B,'LineWidth',1.5);
% xlabel('$\phi_\lambda$','Interpreter','latex','FontSize',11);
% xlim([a b]);
% ylabel('Utility','FontSize',11);
% title('Unconventional Policy: Fiscal Shock','FontSize',11);
% 
% print -depsc2 policy_osr_unconventional.eps


u_no_policy = [u_A0 u_P0 u_O0 u_O0];
save data_u_no_policy u_no_policy;

gain_unconventional = [gain_u_A gain_u_P gain_u_O gain_u_O];

optimal_unconventional = [pphi_A0_optimal pphi_P0_optimal pphi_O0_optimal pphi_O0_optimal];

series_unconventional = [pphi_B0;u_A;u_P;u_O;u_O];

save data_series_unconventional series_unconventional;

save data_optimal_unconventional optimal_unconventional;

gain_and_phi_unconventional = [gain_unconventional;optimal_unconventional];

save data_gain_and_phi_unconventional gain_and_phi_unconventional;

% gain_unconventional =
%     0.7870    1.5106    1.6956    0.7844
% optimal_unconventional =
%     0.8421    0.9474    0.8421    0.8421






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






% gain_unconventional_high_rho = [gain_u_A gain_u_P gain_u_O gain_u_B]
% optimal_unconventional_high_rho = [pphi_A0_optimal pphi_P0_optimal pphi_O0_optimal pphi_B0_optimal]
% 
% 
% series_unconventional_high_rho = [pphi_B0;u_A;u_P;u_O;u_B];
% save data_series_unconventional_high_rho series_unconventional_high_rho;
% 
% gain_and_phi_unconventional_high_rho = [gain_unconventional_high_rho;optimal_unconventional_high_rho];
% 
% save data_gain_and_phi_unconventional_high_rho gain_and_phi_unconventional_high_rho;

% gain_unconventional_high_rho =
%     0.7951    1.8884    2.5531    0.7825
% optimal_unconventional_high_rho =
%     0.8421    0.9474    0.7368    0.8421





