clear all; clc


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

% other parameters, useless here
rho_a = 0.7728;
rho_ppi = 0.7728;
rho_omega = 0.7728;
rho_b = 0.9000;
phi = 0;  % the elasticity of B_{t+1} to private liquidity (L_pri = qk*S, where S denotes trading volume)
stdA = 0.0127;
stdP = stdA;
stdO = stdA;
stdB = 0.0100;
pphi = 0;

pam = [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];

% phi: conveitonal policy
% pphi: unconventional policy

% save data_pam pam;





phi_conventional = 0.3; % leaning againg wind on the supply of public liquidity
pam_conventional_ex = [beta delta alpha A omega_bar chi ksi eta sigma eps_min ppi gamma B rho_a rho_ppi rho_omega rho_b phi_conventional stdA stdP stdO stdB pphi];

% phi: conveitonal policy
% pphi: unconventional policy

% save data_pam_conventional_ex pam_conventional_ex;


pphi_unconventional = 0.3; % leaning againg wind on the supply of public liquidity
pam_unconventional_ex = [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_unconventional];

% phi: conveitonal policy
% pphi: unconventional policy

% save data_pam_unconventional_ex pam_unconventional_ex;




NN = 150; % number of grids

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The Figure on Phi_K versus q_K  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

q_K_ub = 1;
q_K_lb = 0.5;
% capital price
q_K = linspace(q_K_lb,q_K_ub,NN);

eps2 = 1./q_K;

% theta_M: proportion of good-quality capital in markets

theta_M = (((eps2./eps_min).^(-1/eta)).*(1-delta))./((((eps2./eps_min).^(-1/eta)).*(1-delta))+ppi);

% eps1: eps-single-star

Phi_KK = q_K./theta_M;


% figure(1)
% plot(q_K,Phi_KK,'-b',q_K,q_K,'--k','linewidth',1.5);
% ylabel('$\Phi_{t}^{K}:\mbox{Shadow Value of High-Quality Assets}$','Interpreter','latex','FontSize',12);
% xlabel('$q_{t}^{K}:\mbox{Asset Price}$','Interpreter','latex','FontSize',12);
% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h_legend = legend('Information Asymmetry','Benchmark');
% legend('Information Asymmetry','Benchmark');
% %ylim([0.3*min(Phi_K) max(Phi_K)]);
% ylim([0.5 1.2]);
% %xlim([min(q_K) max(q_K)]);
% xlim([q_K_lb q_K_ub]);

% print -depsc2 phik_qk_1.eps



% RECURSIVE characterization



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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SCENARIO I:  COMPLETE   INFORMATION %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

qk_lb_ci = 0.75*qk_ss;
qk_ub_ci = 1.12*qk_ss;

qk_ci = linspace(qk_lb_ci,qk_ub_ci, NN);
% NN = 150; % number of grids

eps1_ci = 1./qk_ci; 

zeta_ci = (eps1_ci./eps_min).^(-sigma);

qb_ci = beta.*(1+zeta_ci./(sigma-1));

R_ci = 1./qb_ci;

rk_ci = ((1./beta-(1-delta).*(1-omega_bar)).*R_ci-omega_bar.*(1-delta)).*qk_ci;

B_K_ci = delta./(zeta_ci.^(1-1./sigma))-rk_ci-(1-delta).*omega_bar.*qk_ci;

I_K_ci = (B_K_ci+rk_ci+(1-delta).*omega_bar.*qk_ci).*zeta_ci;

Y_K_ci = rk_ci./alpha;

C_K_ci = Y_K_ci-I_K_ci-(1-qb_ci).*B_K_ci;

N_ci = (((1-alpha)./chi).*(Y_K_ci./C_K_ci)).^(1./(1+ksi));

K_ci = ((alpha.*A./rk_ci).^(1./(1-alpha))).*N_ci;

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

Y_ci = Y_K_ci.*K_ci;

I_ci = I_K_ci.*K_ci;

C_ci = C_K_ci.*K_ci;

BB_ci = B_K_ci.*K_ci;

B_Y_ci = B_K_ci./Y_K_ci;

I_Y_ci = I_K_ci./Y_K_ci;

w_ci = (1-alpha).*Y_ci./N_ci;

Rk_ci = (rk_ci+(1-delta).*qk_ci)./qk_ci;


[error_ci,a_ci] = min(abs(BB_ci-B*ones(1,NN)));

qk_ss_ci = qk_ci(a_ci);
qb_ss_ci = qb_ci(a_ci);
B_ss_ci = BB_ci(a_ci);
C_ss_ci = C_ci(a_ci);
I_ss_ci = I_ci(a_ci);
K_ss_ci = K_ci(a_ci);
Y_ss_ci = Y_ci(a_ci);
N_ss_ci = N_ci(a_ci);
R_ss_ci = R_ci(a_ci);
rk_ss_ci = rk_ci(a_ci);
zeta_ss_ci = zeta_ci(a_ci);
%omega_ss_ci = omega_ci(a_ci);
eps1_ss_ci = eps1_ci(a_ci);
%theta_ss_ci = theta_ci(a_ci);
BY_ss_ci = B_Y_ci(a_ci);
IY_ss_ci = I_Y_ci(a_ci);
w_ss_ci = w_ci(a_ci);
Lambda_ss_ci = 1/C_ss_ci;
Rk_ss_ci = Rk_ci(a_ci);
omega_ss_ci = omega_bar;
V_ss_ci = log(C_ss_ci)-chi*(N_ss_ci^(1+ksi))/(1+ksi);


% ci: ci denotes for complete infomration
steady_state_ci = [qk_ss_ci eps1_ss_ci qb_ss_ci R_ss_ci rk_ss_ci K_ss_ci N_ss_ci Y_ss_ci w_ss_ci I_ss_ci C_ss_ci BY_ss_ci Rk_ss_ci omega_ss_ci];
% save data_steady_state_ci steady_state_ci;


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



% THE INTERVAL FOR FIGURE 2: the demand for public liquidity
qk_lb = 0.68*qk_ss;
qk_ub = 1.2*qk_ss;


% THOSE ARE REGIONS FOR FIGURES 3: comparative statics
% qk_lb = 0.92*qk_ss;
% qk_ub = 1.2*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;
Rk = (rk+omega_bar.*ppi.*qk)./qk+(1+((eps2./eps_min).^(-1*sigma))./(sigma-1)).*R.*(1-delta).*omega_bar.*Phi_K./qk+R.*(1-delta).*(1-omega_bar).*Phi_K./qk;

spread = Rk-R;

spread_liquidity = (rk+(1-delta).*qk)./qk+(1-delta).*(1-omega_bar).*(1./qb-1)-R;

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);
LP_liquidity_proportion = spread_liquidity(a)./spread(a);
V_ss = log(C_ss)-chi*(N_ss^(1+ksi))/(1+ksi);


% % as: as denotes for adverse selection
% steady_state_as = [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 steady_state_as;


%%%%%%%%%%%%
% FIGURES  %
%%%%%%%%%%%%

figure(2)
subplot(1,2,1)
plot(BB_ci,qk_ci,'-.r',BB,qk,'--b',B_ss.*ones(NN,1),linspace(qk_lb,0.7,NN),'-k','linewidth',1.5);
xlabel('Demand for Public Liquidity (B)');
ylabel('Asset Price (q_k)')
xlim([0 max(BB).*1.028]);
legend('Complete Information','Adverse Selection');
%ylim([qk_lb qk_ub]);
%ylim([qk_lb 0.7]);
ylim([0.5 0.7]);

BB_ci_as_low_omega = [BB_ci;R_ci;BB;R];
save data_BB_ci_as_low_omega_haoning BB_ci_as_low_omega

subplot(1,2,2)
plot(BB_ci,R_ci,'-.r',BB,R,'--b',B_ss.*ones(NN,1),linspace(0.985,1.012,NN),'-k','linewidth',1.5);
xlabel('Demand for Public Liquidity (B)');
ylabel('Interst Rate (R^B)');
xlim([0 max(BB).*1.028]);
legend('Complete Information','Adverse Selection');
ylim([0.99 1/beta]);
ylim([0.985 1.012]);

% print -depsc2 D_public_liquidity_qk.eps
