clear all;

% parameters

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;

% 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];
save data_pam_haoning pam

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);
l_ppi_ss = R_ss*beta*(1-theta_ss^(1/eta))*(eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta))*(1-delta)*omega_bar-R_ss*beta*(1+(eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta)))*omega_bar*ppi*qk_ss/Phi_k_ss;
l_omega_ss = R_ss*beta*((eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta)))*(1-delta)*(1-omega_bar);

disp("B/Y is:");
disp(BY_ss);
disp("I/Y is:");
disp(IY_ss);
disp("mean of liquidity preimum is:");
disp(100*l_omega_ss);
disp("mean of safety premium is:");
disp(100*l_ppi_ss);
disp("R^B is:");
disp(R_ss);
disp("N is:");
disp(N_ss);
