
%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

var Y I K Bl R N w C Ce Ch Le Lh Eql ql qk psi lambdaa varphi vz vq theta lambdaz lambdaq gz gq ggamma mue mub muh;

varexo sun evarphi ea epsi ez eq evq etheta evz;

parameters lambdak lambdai IY delta beta zeta qlLeY qlLhY IK KY Nbar Lbar rhovz alpha Rbar thetabar ggammabar gammae gammah lambdaqbar Omega rhoa rhopsi rhoq rhovq rhotheta rhovarphi rhoz sigmaa sigmavarphi sigmaz sigmavz sigmaq sigmavq sigmatheta sigmapsi sigmasun;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

zeta=0.05;
ggammabar=1.0042211;
lambdaqbar=1.012126;  
gammah=0.4976;
gammae=0.6584;
Omega=0.1753;

Lbar=1; % arbitrary
alpha=0.3; % the average labor income share is 70%
Rbar=1.01; % the anual interest rate is 4%
qlLeY = 2.6; % the average land-output ratio is 0.65 at the annual frequency
qlLhY = 5.8011; % the average housing output ratio is 1.45 at the annual frequency
KY=4.6194; % the capital-output ratio is on average 1.15 at the annual frequency
IK=0.2093/4; % the investment-capital ratio is on average 0.209 at the annual frequency  
Nbar=1/4; % average market hours is 25% of time endowment
thetabar=0.75;  % the average nonfarm and nonfiancial businesses' loan-asset ratio is 0.75 at the annual frequency

rhoa=0.9055;
rhoz=0.4263;
rhovz=0.0095;
rhoq=0.5620;
rhovq=0.2949;
rhovarphi=0.9997;
rhopsi=0.9829;
rhotheta=0.9804;
sigmaa=0.1013;
sigmaz=0.0042;
sigmavz=0.0037;
sigmaq=0.0042;
sigmavq=0.0029;
sigmapsi=0.0073;
sigmatheta=0.0112;
sigmavarphi=0.0462;
sigmasun=0.0462;

lambdak=ggammabar*lambdaqbar;
lambdai=lambdak;
IY=IK*KY;
delta=1-(1-IK)*ggammabar*lambdaqbar;
beta=solvebeta(alpha,delta,ggammabar,lambdak,thetabar,Rbar,zeta,qlLeY,KY);

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model;
# lambdaabar=ggammabar/(beta*Rbar)-1;
# phi =qlLeY*(1-beta-zeta*thetabar*ggammabar*(ggammabar-beta*Rbar)/(Rbar*(ggammabar-beta*(1-zeta))))/(beta*alpha);
# lambdazbar=(ggammabar^(1-(1-phi)*alpha))/lambdaqbar^((1-phi)*alpha);
# BY=zeta*thetabar*(ggammabar*qlLeY+KY/lambdaqbar)/(1-(1-zeta)/ggammabar);
# CeY = alpha-IY-BY*(1/ggammabar-1/Rbar);
# ChY = 1-CeY-IY;
# varphibar=(qlLhY/ChY)*ggammabar*(1-ggammabar/Rbar)*(1-gammah/Rbar)/(ggammabar-gammah);
# psibar=((1-alpha)*ggammabar*(1-gammah/Rbar)*(1/ChY))/(Nbar*(ggammabar-gammah));

1=beta*(muh(+1)/muh)*(1/ggamma(+1))*(1+lambdaa(+1))*R(+1);
1=(beta*(mue(+1)/mue)/ggamma(+1)+mub/mue)*R(+1)-beta*(mub(+1)/mue)/ggamma(+1)*R(+2)*(1-zeta);
alpha*Y=Ce+I+ql*(Le-Le(-1))+R*Bl(-1)/ggamma-Bl;
R(+1)*(Bl-(1-zeta)*Bl(-1)/ggamma)=zeta*theta*(Eql*ggamma(+1)*Le+qk(+1)*K/gq(+1));
qk=beta*mue(+1)/mue*(alpha*(1-phi)*Y(+1)/K+qk(+1)/(gq(+1)*ggamma(+1))*(1-delta))+mub/mue*zeta*theta*qk(+1)/gq(+1);
ql=beta*mue(+1)/mue*(alpha*phi*Y(+1)/Le+Eql)+mub/mue*zeta*theta*Eql*ggamma(+1);

ggamma=(gz*(gq^((1-phi)*alpha)))^(1/(1-(1-phi)*alpha));
muh=1/(Ch-gammah*Ch(-1)/ggamma)-beta*gammah*(1+lambdaa(+1))/(Ch(+1)*ggamma(+1)-gammah*Ch);
w=psi/muh;
ql=beta*(muh(+1)/muh)*(1+lambdaa(+1))*Eql+varphi/(muh*Lh);
mue=1/(Ce-gammae*Ce(-1)/ggamma)-beta*gammae/(Ce(+1)*ggamma(+1)-gammae*Ce);
w=(1-alpha)*Y/N;
1=qk*(1-Omega/2*(I/I(-1)*gq*ggamma-lambdai)^2-Omega*(I/I(-1)*gq*ggamma-lambdai)*I/I(-1)*gq*ggamma)+beta*Omega*mue(+1)/mue/gq(+1)/ggamma(+1)*qk(+1)*(I(+1)/I*gq(+1)*ggamma(+1)-lambdai)*(I(+1)/I*gq(+1)*ggamma(+1))^2;
Y=((gz*gq)^(-(1-phi)*alpha/(1-(1-phi)*alpha)))*(((Le(-1))^phi*(K(-1))^(1-phi))^alpha)*N^(1-alpha);
K=(1-delta)*K(-1)/gq/ggamma+(1-Omega/2*(I/I(-1)*gq*ggamma-lambdai)^2)*I;
Y=Ch+Ce+I;
Lbar=Lh+Le;
C=Ch+Ce;
gz=lambdaz*vz/vz(-1);
gq=lambdaq*vq/vq(-1);
log(lambdaa)=(1-rhoa)*log(lambdaabar)+rhoa*log(lambdaa(-1))+sigmaa*ea;
log(varphi)=(1-rhovarphi)*log(varphibar)+rhovarphi*log(varphi(-1))+sigmavarphi*evarphi;
log(psi)=(1-rhopsi)*log(psibar)+rhopsi*log(psi(-1))+sigmapsi*epsi;
log(lambdaz)=(1-rhoz)*log(lambdazbar)+rhoz*log(lambdaz(-1))+sigmaz*ez;
log(vz)=rhovz*log(vz(-1))+sigmavz*evz;
log(lambdaq)=(1-rhoq)*log(lambdaqbar)+rhoq*log(lambdaq(-1))+sigmaq*eq;
log(vq)=rhovq*log(vq(-1))+sigmavq*evq;
log(theta)=(1-rhotheta)*log(thetabar)+rhotheta*log(theta(-1))+sigmatheta*etheta;
log(ql)-log(steady_state(ql))=log(Eql(-1))-log(steady_state(Eql))+sigmasun*sun;
end;

%----------------------------------------------------------------
% 4. Steasy State
%----------------------------------------------------------------

steady_state_model;
ggamma=ggammabar;
lambdaq=lambdaqbar;
R=Rbar;
lambdakss=ggamma*lambdaq;
IYss=IK*KY;
deltapara=1-(1-IK)*lambdakss;
theta=thetabar;
betapara=beta;
lambdaa=ggammabar/(betapara*R)-1;
phipara=qlLeY*(1-betapara-zeta*theta*ggamma*(ggamma-betapara*R)/(R*(ggamma-betapara*(1-zeta))))/(betapara*alpha);
lambdaz=(ggamma^(1-(1-phipara)*alpha))/lambdaq^((1-phipara)*alpha);
gz=lambdaz;
gq=lambdaq;
vz=1;
vq=1;
BYss=zeta*theta*(ggamma*qlLeY+KY/lambdaq)/(1-(1-zeta)/ggamma);
CeYss= alpha-IK*KY-BYss*(1/ggamma-1/R);
ChYss=1-CeYss-IK*KY;
N=Nbar;
psi=((1-alpha)*ggamma*(1-gammah/R)*(1/ChYss))/(N*(ggamma-gammah));
varphi=(qlLhY/ChYss)*ggamma*(1-ggamma/R)*(1-gammah/R)/(ggamma-gammah);
LhLe=qlLhY/qlLeY;
Le=Lbar/(1+LhLe);
Lh=Lbar-Le;
Y=((gz*gq)^(-(1-phipara)*alpha/(1-(1-phipara)*alpha))*(Le^phipara*KY^(1-phipara))^alpha*N^(1-alpha))^(1/(1-(1-phipara)*alpha));
I=IK*KY*Y;
K=KY*Y;
B=BYss*Y;
Bl=B/R;
Ce=CeYss*Y;
Ch=ChYss*Y;
C=Ce+Ch;
ql=qlLeY*Y/Le;
Eql=ql;
w=(1-alpha)*Y/N;
mue=(ggamma-betapara*gammae)/(Ce*ggamma-gammae*Ce);
mubmue=(ggamma-betapara*R)/(R*(ggamma-betapara*(1-zeta)));
mub=mubmue*mue;
muh=(ggamma-betapara*gammah*(1+lambdaa))/(Ch*ggamma-gammah*Ch);
qk=(betapara*alpha*(1-phipara)*Y/K)/(1-betapara*(1-deltapara)/(gq*ggamma)-mubmue*zeta*theta/gq);
end;

%----------------------------------------------------------------
% 5. Computation
%----------------------------------------------------------------

steady;

check;

shocks;
var sun=1;
%var evarphi=1;
%var ea=1;
%var epsi=1;
%var ez=1;
%var eq=1;
%var evq=1;
%var evz=1;
%var etheta=1;
end;

stoch_simul(loglinear,irf=200,order=1,periods=10000);
datatomfile('simudata',[]);
%rplot  Y I K Bl R N w C Ce Ch Le Lh Eql ql qk mue mub muh;

SS=oo_.steady_state;
MEAN=oo_.mean;
SD=sqrt(diag(oo_.var));
SDM=100*SD./MEAN;
SDY=100*SD/SD(1);
COVY=oo_.var(:,1);
CORRY=(COVY./SD)/SD(1);
tablevalue=[SDY CORRY];
dyntable('STATISTICS SUMMARIZE (10000 PERIODS)',strvcat('VARIABLE (LOG LEVEL)','S.D. RELATIVE TO THAT OF Y in %', 'CORRELATION WITH Y'),M_.endo_names,tablevalue,10,8,4);

