clear all;
clc;
% alpha = 0.33; % production function POE
alpha = 0.4; % production function POE
gamma = 0.5; % production function SOE
beta = 0.995; % discount factor
% delta = 0.025; % steady state depreciation rate
% Deltas1 = 2; % steady state difference in recoverability POE vs SOE
% Deltas2 = 1/3; % steady state difference in recoverability government vs SOE
% sp = 0.03; % steady state net spread SOE
delta = 0.035; % steady state depreciation rate, following chang et al. 2019
Deltas1 = 35/25; % steady state difference in recoverability POE vs SOE, according to Geng and Pan (JF forthcoming, the SOE premium is around 100bp)
Deltas2 = 5/25; % steady state difference in recoverability government vs SOE
sp = 0.025; % steady state net spread SOE, according to Geng and Pan (JF forthcoming)
Pis = 1; % steady state gross inflation
kappaf = 1-40^(-1); % decay parameter private
kappab = 1-40^(-1); % decay parameter government
epsips = 11; % elasticity goods SOE
epsipp = 11; % elasticity goods POE
epsiws = 11; % elasticity labor SOE
epsiwp = 11; % elasticity labor POE
phips = 0.75; % price stickiness SOE
phipp = 0.75; % price stickiness POE
phiws = 0.75; % wage stickiness SOE
phiwp = 0.75; % wage stickiness POE
gammaps = 0; % price indexation SOE
gammapp = 0; % price indexation POE
gammaws = 0; % wage indexation SOE
gammawp = 0; % wage indexation POE
% gy = 0.33; % government spending share 
gy = 0.17; % government spending share 
bh = 0; % habit formation 
eta = 1; % inverse Frisch elasticity
ffs = 0; % fraction of SOE debt held by central bank
ffp = 0; % fraction of POE debt held by central bank
levs = 4; % steady state leverage 
sigma = 0.95; % survival probability
% psi = .45; % SOE investment constraint
psi = .2; % SOE investment constraint
bgy = 0.77*4; % steady state debt-gdp ratio

% Monetary policy rule
rhor = 0.8;
phipi = 1.5;
phiy = 0.25;
phistruc = 5;

% investment and utilization adjustment costs
kappai = 2;

% aggregation of output and labor
% rhoL = 0.5; % we should choose rhoL to pin down the sectoral labor
% rhoY = 0.5; % we should choose rhoY to pin down the sectoral output and
% price ratio
sigmaL = 2;
sigmaY = 2;

% shock parameters
rhoA = 0.9;
rhoxi = 0.9;
rhoG = 0.9;
rhoB = 0.98;
rhofs = rhor;
rhofp = rhor;
rhob = rhor;
rhot = 0.98;
rhoD1 = 0.8;
rhoD2 = 0.98;
sA = 0.0065;
sxi = 0.01;
sG = 0.01;
st = 0.04;
sD1 = 0.0025; 
sD2 = 0.00; 
sB = 0.00;

Ldss = 0.25;
Ldps = 1;
pks = 1;

rhof_QTs = rhor;
rhof_QTp = rhor;
rhob_QT = rhor;
fracb = 0;
fracfs = 0;
fracfp = 0;
sr =  0.0025; 
sfs = sr; 
sfp = sr;
sb = sr;

ratio = 1.5; %pp*yp/(ps*ys)
xis = 0.75;

psifs=0;
psifp=0;
psib=0;


save param alpha gamma beta delta Deltas1 Deltas2 sp Pis kappaf kappab epsips epsipp epsiws epsiwp ...
            phips phipp phiws phiwp gammaps gammapp gammaws gammawp gy bh eta ffs ffp levs sigma ...
            sigmaL sigmaY psi bgy rhor phipi phiy kappai rhoA rhoxi ...
            rhoG rhoB rhofs rhofp rhob rhot rhoD1 rhoD2 sr sA sxi sG sfs sfp st sD1 sD2 sb ...
            sB psifs psifp psib Ldss Ldps pks rhof_QTs rhof_QTp rhob_QT fracb  fracfs fracfp ratio xis

param_ss

dynare model_baseline_dydel

% load figures\irf_et_base.mat;
% 
% %% figures
% 
% % Compare irf including or not structure
% figure('Position',[50,70,1600,900]);
% subplot(3,4,1);
% plot(1:40,100*Y_et','-',1:40,100*irf_et_base(:,1)','--','LineWidth',1.5);
% title('Aggregate Output');
% grid on;
% box on;
% subplot(3,4,2);
% plot(1:40,100*I_et','-',1:40,100*irf_et_base(:,2)','--','LineWidth',1.5);
% title('Aggregate Investment');
% grid on;
% box on;
% subplot(3,4,3);
% plot(1:40,100*C_et','-',1:40,100*irf_et_base(:,3)','--','LineWidth',1.5);
% title('Aggregate Consumption');
% grid on;
% box on;
% subplot(3,4,4);
% plot(1:40,400*infl_et','-',1:40,400*irf_et_base(:,4)','--','LineWidth',1.5);
% title('Inflation');
% grid on;
% box on;
% subplot(3,4,5);
% plot(1:40,400*RTR_et','-',1:40,400*irf_et_base(:,5)','--','LineWidth',1.5);
% title('Policy Rate');
% grid on;
% box on;
% subplot(3,4,6);
% plot(1:40,400*rlfs_real_et','-',1:40,400*irf_et_base(:,6)','--','LineWidth',1.5);
% title('Real long yield SOE');
% grid on;
% box on;
% subplot(3,4,7);
% plot(1:40,400*rlfp_real_et','-',1:40,400*irf_et_base(:,7)','--','LineWidth',1.5);
% title('Real long yield POE');
% grid on;
% box on;
% subplot(3,4,8);
% plot(1:40,400*spreadfps_et','-',1:40,400*irf_et_base(:,8)','--','LineWidth',1.5);
% title('Spread between POE and SOE');
% grid on;
% box on;
% subplot(3,4,9);
% plot(1:40,LsLp_et','-',1:40,irf_et_base(:,9)','--','LineWidth',1.5);
% title('Labor ratio L_s/L_p');
% grid on;
% box on;
% subplot(3,4,10);
% plot(1:40,KsKp_et','-',1:40,irf_et_base(:,10)','--','LineWidth',1.5);
% title('Capital ratio K_s/K_p');
% grid on;
% box on;
% subplot(3,4,11);
% plot(1:40,YsYp_et','-',1:40,irf_et_base(:,11)','--','LineWidth',1.5);
% title('GDP ratio GDP_s/GDP_p');
% grid on;
% box on;
% legend('Strucutral Policy','Baseline','Location','east');
% 
% print('-vector','figures/credit_shock_irf_compare_dydel','-depsc');
% savefig('figures/credit_shock_irf_compare_dydel');

%% Welfare Analysis
clear all;
NN = 100;
phistruc_vec = linspace(0,60,NN)';
wel_vec = zeros(NN,1);

% base case
phistruc = 0;
save phist phistruc
dynare model_wel_dydel;
index_w = strmatch('welf', M_.endo_names, 'exact'); 
index_wc = strmatch('welc', M_.endo_names, 'exact');
index_wl = strmatch('well', M_.endo_names, 'exact'); 
index_corr = oo_.dr.inv_order_var(index_w); 
w          = oo_.steady_state(index_w,1);
correction = oo_.dr.ghs2(index_corr,1)/2;
% wel_base = oo_.mean(index_w);
wel_base = w+correction;
index_w = strmatch('welf', M_.endo_names, 'exact'); 
for i = 1:NN
    M_.params(strmatch('phistruc',M_.param_names,'exact'))=phistruc_vec(i);
    [info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, {'welf'});
    index_corr = oo_.dr.inv_order_var(index_w); 
    w          = oo_.steady_state(index_w,1);
    correction = oo_.dr.ghs2(index_corr,1)/2;
	wel = w + correction;
    % wel_vec(i) = oo_.mean(index_w);
    wel_vec(i) = w + correction;
end
wel_vec = 100*(exp((1-beta)*(wel_vec-wel_base))-1);
figure;
plot(phistruc_vec,wel_vec,'LineWidth',1.5);
xlabel('$\phi_{struc}$','Interpreter','latex');
ylabel('Social Welfare');
grid on;
box on;
