%Declare endogenous variables 
var ak az cb ce ch cs d epse epss he hh hs ke kh lambdab lambdae lambdas le ls nh ns q re rh rk rm rs rv wh ws x y zke zkh;
%Declare exogenous variables
varexo ve vk vs vz;

%Declare parameters and their values
parameters alpha betab betae betah betas bk delta eta j gammae gammas mh mk mn ms mu rhoak rhoaz rhobe rhobh rhoe rhod rhos sigma sigmave sigmavk sigmavs sigmavz tauh taus v zetae zetah;
alpha=0.35;
betab=0.945;
betae=0.94;
betah=0.9925;
betas=0.94;
bk=0.0122104818154666;
delta=0.035;
gammae=0.9;
gammas=0.9;
j=0.075;
mh=0.9;
mk=0.9;
mn=1;
ms=0.9;
mu=0.464053954028799;
rhoak=0.907564408672319;
rhoaz=0.990599723824274;
rhobe=0.930835177208793;
rhobh=0.969199134850878;
rhoe=0.05;
rhod=0.650291811303701;
rhos=0.700910340904156;
sigma=0.326564481162563;
sigmave=0.00105366310939923; 
sigmavk=0.00675523019012514; 
sigmavs=0.00132513108400193; 
sigmavz=0.00699646242063357;
tauh=0.2;
taus=0.2;
v=0.0378989733672577;
zetae=0.417264168830014;
zetah=0.376780917346019;

model;
    %Equations for household saver (1-5)
exp(ch)=exp(rm+zkh+kh(-1))+(exp(rh(-1))*exp(d(-1)))+(exp(wh)*exp(nh))-exp(d)-(exp(q)*(exp(hh)-exp(hh(-1))))-exp(kh);
1/exp(ch)=betah*exp(rh-ch(+1));
exp(wh-ch)=tauh/(1-exp(nh));
exp(q-ch)=(j/exp(hh))+betah*exp(q(+1)-ch(+1));
(1/exp(ak))*(1/exp(ch))=betah*((exp(rm(+1)+zkh(+1)))*(1/exp(ch)));

    %Equations for household borrower (6-10)
exp(cs)=exp(ls)+(exp(ws)*exp(ns))-exp(rs(-1)+ls(-1))-(exp(q)*(exp(hs)-exp(hs(-1))))-exp(epss);
exp(ls)=(rhos*exp(ls(-1)))+(1-rhos)*((ms)/exp(rs))*exp(q(+1))*exp(hs);
(1/exp(cs))=betas*exp(rs-cs(+1))+exp(lambdas+rs);
exp(ws-cs)=taus/1-exp(ns);
exp(q-cs)=(j/exp(hs))+betas*exp(q(+1)-cs(+1))+exp(lambdas+q(+1));

    %Equations for bank (11-15)
exp(cb)=exp(d)+(exp(re)*exp(le(-1)))+(exp(rs)*exp(ls(-1)))-(exp(rh(-1))*exp(d(-1)))-exp(le)-exp(ls)-exp(epse)-exp(epss);
exp(le)+exp(ls)-exp(d)-exp(epse(-1))-exp(epss(-1))=rhod*(exp(le(-1))+exp(ls(-1))-exp(d(-1))-exp(epse)-exp(epss))+(1-gammae)*(1-rhod)*(exp(le)+exp(ls)-exp(epse(+1))-exp(epss(+1)));
(1-exp(lambdab))*(1/exp(cb))=betab*((exp(rh)-rhod*exp(lambdab(+1)))*(1/exp(cb(+1))));
(1-(gammae*(1-rhod)+rhod)*exp(lambdab))*(1/exp(cb))=betab*((exp(re(+1))-rhod*exp(lambdab(+1)))*(1/exp(cb(+1))));
(1-(gammas*(1-rhod)+rhod)*exp(lambdab))*(1/exp(cb))=betab*((exp(rs)-rhod*exp(lambdab(+1)))*(1/exp(cb)));

    %Equations for entrepreneurs (16-21)
exp(ce)=exp(y)+exp(le)-exp(ke)-(exp(q)*(exp(he)-exp(he(-1))))-(exp(re)*exp(le(-1)))-(exp(wh)*exp(nh))-(exp(ws)*exp(ns))-exp(rm+zkh+kh(-1))+exp(epse);
y=exp(az)*((exp(zkh)*exp(kh(-1)))^(alpha*(1-mu)))*((exp(zke)*exp(ke(-1)))^(alpha*mu))*((exp(he(-1)))^v)*((exp(nh))^((1-alpha-v)*(1-sigma)))*((exp(ns))^((1-alpha-v)*sigma));
exp(le)=(rhoe*exp(le(-1)))+(1-rhoe)*(mh/exp(re(+1)))*exp(q(+1)+he)+mk*(exp(ke))-mn*((exp(wh)*exp(nh))+(exp(ws)*exp(ns)));
(exp(q)-exp(lambdae)*(1-rhoe)*mh*exp(q(+1))/exp(re(+1)))*(1/(exp(ce)))=betae*(exp(q(+1))*(1+exp(rv(+1)))*(1/(exp(ce(+1)))));
(1-exp(lambdae))*(1/exp(ce))=betae*((exp(re(+1))-(rhoe*exp(lambdae)))*(1/(exp(ce(+1)))));
((1-exp(lambdae))*(1-rhoe)*mk)*(1/exp(ce))=betae*(((1+exp(rk(+1))+exp(zke(+1))))*(1/exp(ce(+1))));

    %Other equations (22-29)
alpha*mu*exp(y)=exp(rk+zke+ke(-1));
alpha*(1-mu)*exp(y)=exp(rm+zkh+kh(-1)+x);
v*exp(y)=exp(rv+q+he(-1));
(1-alpha-v)*(1-sigma)*exp(y)=exp(wh+nh)*(1+mn*exp(lambdae));
(1-alpha-v)*sigma*exp(y)=exp(ws+ns)*(1+mn*exp(lambdae));
bk=(1/betae)+1-delta;
exp(rk)=bk*(((zetae/(1-zetae))*exp(zke))+1-(zetae/(1-zetae)));
exp(rm)=((1/(betah-1+delta))*((zetah/(1-zetah))*zkh)) + 1-((zetah)/(1-zetah));

    %Shock Equations (30-33)
exp(epse)=rhobe*(exp(epse(-1)))+ve; //% Shock to repayment of subprime guy
exp(epss)=rhobh*(exp(epss(-1)))+vs; //% Shock to repayment of subprime guy
exp(ak)=rhoak*exp(ak(-1))+vk;
exp(az)=rhoaz*exp(az(-1))+vz;

    %Equilibrium (34)
exp(hh)+exp(hs)+exp(he)=1;
end;

initval;
ak=0.01;
az=0.01;
cb=0.01;
ce=0.01;
ch=0.01;
cs=0.01;
d=0.01;
epse=0.01;
epss=1;
he=0.01;
hh=0.01;
hs=0.01;
ke=0.01;
kh=0.01;
lambdab=0.01;
lambdae=0.01;
lambdas=0.01;
le=0.03;
ls=0.03;
nh=0.0378989733672577;
ns=0.01;
q=0.01;
re=0.05;
rh=0.03;
rk=0.01;
rm=0.01;
rs=0.03;
rv=0.01;
wh=0.01;
ws=0.01;
x=0;
y=0.01;
zke=0.01;
zkh=0.01;
end;


    %calculations of Steady State
steady;
check;


    %stochastic properties of shocks
shocks;
var ve; stderr sigmave;
var vs; stderr sigmavs;
var vk; stderr sigmavk;
var vz; stderr sigmavz;
end;

stoch_simul(order=1,irf=20,nomoments,noprint) y rh rs re ch ce cb cs q hh he hs le ls d;