%zhang and JPE paper model with capiital and assets same question as firm_vol and taxonprofits mod file
%but unlike taxonprofits file, there is a rep. capiital producer to piin down
%aggregate K and endogenize pii, needs fsolve helper. Model solved, replace the helper equation in terms of original parameters not deep ones (solved in SnIB diary)
% calibrating the FD parameter - data equivalent? - so go back to taxonprofits.mod file - intro investment adjustment costs to rectify the negative corr (Y,I)?


clc;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

var ke, kh, K, ae, ah, le, lh, ls, ye, yh, Y, Ae, Ah, pii, ce, ch, C, w, r, mu1, mu2 I, d, z, MUe, MUh; 
predetermined_variables K;
varexo e eb;
parameters cai_0, cai_1, a1, a2, a3, a4, Ey, Eke, Ekh, psi $\psi$, rhob $\rho_b$, gamma $\gamma$, beta $\beta$, sigma $\sigma$, lambda $\lambda$, delta $delta$, alpha $\alpha$, rho $\rho$, eta $\eta$; 

beta=0.99;
gamma=0.96;
sigma=0.006;
sigmab=0.007;
alpha=0.33;
delta=0.075;
rho=0.95;
rhob=0.95;
rhobh=0.95;
psi=1;
cai_0=0.73;
cai_1=0.1;
rhoj=0.95;
eta=2.75;
lambda=1.25;

//Computation of the Steady States

Eke=lambda*alpha*beta/((beta-gamma)+lambda*beta*(((1/beta)-1)+delta));
Ekh=alpha*beta/(((1/beta)-1)+delta);
a1=((1-alpha)+(((1/beta)-1)*(lambda-1)/lambda)*Eke)/(Eke-Ekh);
a2=((1-alpha)+(((1/beta)-1)*(lambda-1)/lambda)*Eke)*Ekh/(Eke-Ekh);
a3=(1-delta*Ekh)*Eke/(Eke-Ekh);
a4=(1-delta*Ekh)/(Eke-Ekh);
Ey=((Ekh^alpha)^((1+eta)/(eta*(1-alpha)))*((1-alpha)/psi))^(1/(1+eta));

model;

//F.O.C consumptions
MUe=1/ce;
MUh=1/ch;
mu2*ce=0;

//F.O.C. assets
MUe=gamma*(1+r(+1))*MUe(+1)+mu1(+1)*lambda;
beta*MUh(+1)*(1+r(+1))=MUh;

//F.O.C capiital
alpha*yh/kh=(r+delta);                %entrepreneur
MUe*(alpha*ye/ke-(delta+r))-mu1=0;    %listed firms
r=((K(+1)-(1-delta)*K)/cai_0)^(1/(1-cai_1))-((K(+1)-(1-delta)*K)/cai_0)^(cai_1/(1-cai_1))*((1-delta)/(cai_0*(1-cai_1)))*K;

//F.O.C labor 
((1-alpha)*ye)/le=w;
(1-alpha)*yh/lh=w;
w/ch=psi*(ls^(1/eta));

//Investments
K(+1)=(1-delta)*K+cai_0*(I/K)^(1-cai_1);
/////

//Budget Constraints
ce+ae(+1)=ye-w*le-(delta+r)*ke+(1+r)*ae;
ch+ah(+1)=d+(1+r)*ah+w*ls;
d=yh-w*lh-(r+delta)*kh;

//Borrowing Constraints
ke=lambda*ae;

//Output
ye=z*Ae*(ke^alpha)*le^(1-alpha);
yh=z*Ah*(kh^alpha)*lh^(1-alpha);

//Aggregate
Ah=1;
ae+ah=K;
K=pii*Y;
ke+kh=K;
le+lh=ls;
ce+ch=C;
ye+yh=Y;

log(z) = (rho*log(z(-1)))+(e);
end;

steady_state_model;
z=1;
r=(1/beta)-1;
Ah=1;
Ae=((Ekh/Eke)^alpha)*Ah;
K=((cai_0*(1-cai_1)*r/(1-delta))*(1/((delta/cai_0)^(1/(1-cai_1))-(delta/cai_0)^(cai_1/(1-cai_1)))))^(1-cai_1);
pii=firm_vol2_2_steady_state_helper(cai_0, cai_1, beta, delta, eta, alpha, psi, gamma, lambda);
Y=K/pii;
ye=(pii-Ekh)*Y/(Eke-Ekh);
yh=Y-ye;
ch=((1-alpha)+r*((lambda-1)/lambda)*Eke)*ye+(1-delta*Ekh)*yh;
%ch=(((Ah*Ekh^alpha)^((1+eta)/(1-alpha))*(1/Y))^(1/eta))*(1-alpha)/psi;
w=(((psi*ch)^eta)*(1-alpha)*Y)^(1/(1+eta));
%w=(Ae*Eke^alpha)^(1/(1-alpha))*(1-alpha);
le=(1-alpha)*ye/w;
lh=(1-alpha)*yh/w;
ls=le+lh;
ke=Eke*ye;
kh=Ekh*yh;
ae=ke/lambda;
ah=Eke*ye+Ekh*yh-ae;
I=(delta/cai_0)^(1/(1-cai_1))*K^((2-cai_1)/(1-cai_1));
d=yh-w*lh-(r+delta)*kh;
ce=ye-w*le-(delta+r)*ke+r*ae;
MUe=1/ce;
MUh=1/ch;
mu1=(beta-gamma)*MUe/(beta*lambda);
mu2=0;
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4. Simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

shocks;
var e= 0.006^2;
%periods 1:9;
%values 0.1;

%var ebe=0.002^2;
%var ebh=(0.004)^2;
%var eie = 0.005^2;
%var e, eie = cai*0.004*0.005;
%var eie = 0.005^2;
%var eih = 0.005^2;
%var eie, eih = cai*0.005*0.005;
end;

resid(1);
steady;
%solve algo=1;
check;

stoch_simul(order=2, hp_filter=1600, periods=10000, drop=200, irf=0);
%write_latex_dynamic_model ;

%simul (periods=2100);
// Some Results

statistic1 = 100*sqrt(diag(oo_.var(1:26,1:26)))./oo_.mean(1:26);
dyntable('Relative standard deviations',strvcat('VARIABLE','REL. S.D.'),M_.endo_names(1:26,:),statistic1,10,8,4);