var  
c,i,                  
R,R_s,
b,b_s, 
w,l,  
q,rk,
g1,g2,
n,d,    
k,PI,
PIH,ph,
PIs,re,
mc,y,
yH,g,
de,rer,
ca,tb,
x,x_s,
tax,ca_y;

varexo
phit y_ss b_ss
g_ss xs_ss;

parameters
v,sigmal,
sigmac,h,
bet,delta,
lambda,alp,
phi,deltah,
theta,thetah,
rhors,R_ss,
PI_ss,Rs_ss,
rhor,rhoy,
rhopi,
h1,h2,eta,
rhob,by,g_y,
gammah,gamma,
sx,sm;

sigmal=1.3;
sigmac=2;
h=0.6;
bet=0.998;
delta=0.025;
lambda=5;
alp=0.45;
phi=7;
deltah=0.5;
theta=0.67;
thetah=2;
eta=0.99;
rhob=0.80;
by=2.5;
rhors=0.98;
sx=0.15;
sm=0.12;
gammah=0.80;
gamma=2;
v=30;

rhor=0.90;
rhoy=0.10;
rhopi=0.7;

h1=-0.1;
h2=-0.4;

PI_ss=1.0075;
g_y=0.17;
Rs_ss=1.01;
R_ss=PI_ss/bet;


model;

v*l^sigmal=(c-h*c(-1))^(-sigmac)*w;

(c-h*c(-1))^(-sigmac)=bet*R*(c(+1)-h*c)^(-sigmac)/PI(+1);

k=(1-delta)*k(-1)+(1-(lambda/2)*(i/i(-1)-1)^2)*i;

1=q*((1-(lambda/2)*(i/i(-1)-1)^2)-lambda*(i/i(-1)-1)*(i/i(-1)))
+lambda*q(+1)*(PI(+1)/R)*(i(+1)/i-1)*(i(+1)/i)^2;

q*k=n+d;

re=PI*(rk+(1-delta)*q)/q(-1);

R(-1)*d(-1)/PI=phit*(rk+(1-delta)*q)*k(-1);

n=eta*(re*k(-1)*q(-1)/PI-R(-1)*d(-1)/PI);

yH=k(-1)^alp*l^(1-alp);

k(-1)/l=(alp*w/((1-alp)*rk));

mc=alp^(-alp)*(1-alp)^(-(1-alp))*rk^alp*w^(1-alp);

g1=(c-h*c(-1))^(-sigmac)*yH*PIs+bet*theta*(PIs/PIs(+1))
*(PIH^deltah/PI(+1))^(1-phi)*g1(+1);

g2=(phi/(phi-1))*(c-h*c(-1))^(-sigmac)*yH*mc+bet*theta
*(PIH^deltah/PI(+1))^(-phi)*g2(+1);

g1=g2;

ph^(1-phi)=theta*(ph(-1)*PIH(-1)^deltah/PI)^(1-phi)+(1-theta)*PIs^(1-phi);

1=gammah*ph^(1-thetah)+(1-gammah)*rer^(1-thetah);

ph/ph(-1)=PIH/PI;

b=d+rer*b_s;

yH=gammah*ph^(-thetah)*(c+i+g)+x;

y=c+i+g+tb;

g=tax+(de*R_s(-1)-R(-1))*rer(-1)*b_s(-1)/PI;

rer*b_s=de*rer(-1)*R_s(-1)*b_s(-1)/PI+tb;

tb=ph*x-(1-gammah)*rer^(1-thetah)*(c+i);

x=rer^gamma*x_s;

ca=rer*b_s-de*rer(-1)*b_s(-1)/PI;

ca_y=ca/y;

log(R_s)=rhors*log(R_s(-1))+(1-rhors)*log(Rs_ss);

log(x_s)=rhors*log(x_s(-1))+(1-rhors)*log(xs_ss);

R/R_ss=(R(-1)/R_ss)^rhor*((y/y_ss)^rhoy
*(PI/PI_ss)^rhopi)^(1-rhor);

g/g_ss=(y/y_ss)^h1*(b(-1)/b_ss)^h2;

b=rhob*b(-1)+(1-rhob)*by*y;

g=rhob*g(-1)+(1-rhob)*g_y*y;

end;

initval;
q=1;
PI=PI_ss;
R=PI/bet;
phit=0.3;
rer=(sm/((1-gammah)*(1-g_y-sx+sm)))^
(1/(1-thetah));
ph=((1-(1-gammah)*rer^(1-thetah))/
gammah)^(1/(1-thetah));
PIH=PI;
PIs=(ph^(1-phi)*(1-theta*PI^((deltah-1)*
(1-phi)))/(1-theta))^(1/(1-phi));
mc=((phi-1)/phi)*((1-bet*theta*PI^
((deltah-1)*(-phi)))/(1-bet*theta*PI^
((deltah-1)*(1-phi))))*PIs;
rk=(eta*(1-phit)+bet*phit)^(-1)-(1-delta);
R_s=Rs_ss;
w=(mc*alp^alp*(1-alp)^(1-alp)*rk^(-alp))^(1/(1-alp));
re=PI*(rk+1-delta);
y=((1-h)^(-sigmac)*(sm/((1-gammah)*rer^(1-thetah))-
(delta*(gammah*ph^(-thetah)*(sm/((1-gammah)*
rer^(1-thetah))+g_y)+sx/ph)*(alp*w/((1-alp)*rk))^
(1-alp)))^(-sigmac)*w/(v*((gammah*ph^(-thetah)*
(sm/((1-gammah)*rer^(1-thetah))+g_y)+sx/ph)*
(alp*w/((1-alp)*rk))^(-alp))^sigmal))^(1/(sigmal+sigmac));
yH=(gammah*ph^(-thetah)*(sm/((1-gammah)*
rer^(1-thetah))+g_y)+sx/ph)*y;
l=yH*(alp*w/((1-alp)*rk))^(-alp);
k=(alp*w/((1-alp)*rk))*l;
i=delta*k;
c=(sm/((1-gammah)*rer^(1-thetah))-i/y)*y;
g=g_y*y;
x=yH-gammah*ph^(-thetah)*(c+i+g);
x_s=rer^(-gamma)*x;
tb=(sx-sm)*y;
d=PI*phit*k*(rk+1-delta)/R;
n=k-d;
b=by*y;
g1=(c-h*c)^(-sigmac)*yH*PIs/
(1-bet*theta*PI^((deltah-1)*(1-phi)));
g2=g1;
b_s=(b-d)/rer;
de=PI*(1-tb/(rer*b_s))/R_s;
ca=(1-de/PI)*b_s*rer;
tax=g-rer*b_s*(de*R_s-R)/PI;
ca_y=ca/y;
b_ss=b;
y_ss=y;
g_ss=g;
xs_ss=x_s;
end;

steady;

endval;
q=1;
PI=PI_ss;
R=PI/bet;
phit=0.2;
rer=(sm/((1-gammah)*(1-g_y-sx+sm)))^
(1/(1-thetah));
ph=((1-(1-gammah)*rer^(1-thetah))/
gammah)^(1/(1-thetah));
PIH=PI;
PIs=(ph^(1-phi)*(1-theta*PI^((deltah-1)*
(1-phi)))/(1-theta))^(1/(1-phi));
mc=((phi-1)/phi)*((1-bet*theta*PI^
((deltah-1)*(-phi)))/(1-bet*theta*PI^
((deltah-1)*(1-phi))))*PIs;
rk=(eta*(1-phit)+bet*phit)^(-1)-(1-delta);
R_s=Rs_ss;
w=(mc*alp^alp*(1-alp)^(1-alp)*rk^(-alp))^(1/(1-alp));
re=PI*(rk+1-delta);
y=((1-h)^(-sigmac)*(sm/((1-gammah)*rer^(1-thetah))-
(delta*(gammah*ph^(-thetah)*(sm/((1-gammah)*
rer^(1-thetah))+g_y)+sx/ph)*(alp*w/((1-alp)*rk))^
(1-alp)))^(-sigmac)*w/(v*((gammah*ph^(-thetah)*
(sm/((1-gammah)*rer^(1-thetah))+g_y)+sx/ph)*
(alp*w/((1-alp)*rk))^(-alp))^sigmal))^(1/(sigmal+sigmac));
yH=(gammah*ph^(-thetah)*(sm/((1-gammah)*
rer^(1-thetah))+g_y)+sx/ph)*y;
l=yH*(alp*w/((1-alp)*rk))^(-alp);
k=(alp*w/((1-alp)*rk))*l;
i=delta*k;
c=(sm/((1-gammah)*rer^(1-thetah))-i/y)*y;
g=g_y*y;
x=yH-gammah*ph^(-thetah)*(c+i+g);
x_s=rer^(-gamma)*x;
tb=(sx-sm)*y;
d=PI*phit*k*(rk+1-delta)/R;
n=k-d;
b=by*y;
g1=(c-h*c)^(-sigmac)*yH*PIs/
(1-bet*theta*PI^((deltah-1)*(1-phi)));
g2=g1;
b_s=(b-d)/rer;
de=PI*(1-tb/(rer*b_s))/R_s;
ca=(1-de/PI)*b_s*rer;
tax=g-rer*b_s*(de*R_s-R)/PI;
ca_y=ca/y;
b_ss=b;
y_ss=y;
g_ss=g;
xs_ss=x_s;
end;

steady;


shocks;
var phit;
periods 1:5;
values 0.3;
end;

simul(periods=200);

var_str={'c','i','R','b','b_s','n','d', ...
'k','PI','PIs','re','y','tb','ca','ca_y'};

for ii=1:size(var_str,2)
    if ii>=1 && ii<10
	figure(1)
	subplot(3,3,ii)
	plot(1:200,oo_.endo_simul(strmatch(var_str{ii}, ...
	M_.endo_names,'exact'),2:201),'b-','linewidth',1.5);
	title(var_str(:,ii));
	
	elseif ii>=10 && ii<19
	figure(2)
	subplot(3,3,ii-9)
	plot(1:200,oo_.endo_simul(strmatch(var_str{ii}, ...
	M_.endo_names,'exact'),2:201),'b-','linewidth',1.5);
	title(var_str(:,ii));
	
	end
end