var dc dr dm dpi din r h x a e z v y k in c d m w r q n tau lambd ksi pi c_obs m_obs in_obs pi_obs r_obs;
varexo epsilona epsilone epsilonx epsilonz epsilonv;
parameters rhoa rhoe rhox rhoz rhov alph phip phik lne lnz sigmaa sigmae sigmax sigmaz sigmav omegapi omegay omegatau bet et delt thet gamm ctrend constepi conster intrend mtrend;

rhoa=0.75;
rhoe=0.75;
rhox=0.75;
rhoz=0.75;
rhov=0.50;
alph=0.3;
phip=50;
phik=32.14;
%lne=1.39;
%lnz=8.13;
sigmaa=0.01;
sigmae=0.01;
sigmax=0.01;
sigmaz=0.01;
sigmav=0.01;
omegapi=1.3;
omegay=0.125;
omegatau=0.38;
bet=0.99;
et=1.5;
delt=0.025;
thet=6;
gamm=0.042;

ctrend=0;
constepi=0;
conster=0;
ctrend=0;
intrend=0;
mtrend=0;

model;
k=(1-delt)*k(-1)+x*in;
ln(x)=rhox*ln(x(-1))+epsilonx;
ln(a)=rhoa*ln(a(-1))+epsilona;
ln(e)=(1-rhoe)*ln(steady_state(e))+rhoe*ln(e(-1))+epsilone;
ln(z)=(1-rhoz)*ln(steady_state(z))+rhoz*ln(z(-1))+epsilonz;
ln(v)=rhov*ln(v(-1))+epsilonv;
a=lambd*c^(1/gamm)*(c^((gamm-1)/gamm)+(e^(1/gamm))*(m^((gamm-1)/gamm)));
lambd*w*(1-h)=et;
c*e=m*(1-1/r)^gamm;
lambd=bet*r*((lambd(+1))/(pi(+1)));
lambd*(1/x+phik*(k/k(-1)-1))=bet*(lambd(+1))*(q(+1)+(1-delt)/x(+1))-bet*phik/2*(lambd(+1)*(k(+1)/k-1)^2)+bet*phik*lambd(+1)*(k(+1)/k-1)*(k(+1)/k);
y=c+in+(phik/2)*((k/k(-1)-1)^2)*k(-1)+(phip/2)*(((pi/steady_state(pi))-1)^2)*y;
d=y-w*h-q*k(-1)-(phip/2)*((pi/steady_state(pi)-1)^2)*y;
lambd*w*h=(1-alph)*ksi*y;
lambd*q*k(-1)=alph*ksi*y;
phip*lambd*(pi/steady_state(pi)-1)*(pi/steady_state(pi))=(1-thet)*lambd+thet*ksi+bet*phip*(lambd(+1))*((pi/steady_state(pi))-1)*(pi/steady_state(pi))*((y(+1))/y);
y=(k(-1)^alph)*(z*h)^(1-alph);
ln(r/steady_state(r))=omegatau*ln(tau/steady_state(tau))+omegapi*ln(pi/steady_state(pi))+omegay*ln(y/steady_state(y))+ln(v);
n=y/h;
tau=(m/m(-1))*pi;
c_obs=log(c)-log(steady_state(c));
m_obs=log(m)-log(steady_state(m));
in_obs=log(in)-log(steady_state(in));
pi_obs=log(pi);
r_obs=log(r)*400;

din=in-in(-1)+intrend;
dc=c-c(-1)+ctrend;
dpi=1*pi+constepi;
dr=1*r+conster;
dm=m-m(-1)+mtrend;

end;

steady_state_model;

a=1;
e=1;
x=1;
z=1;
v=1;
pi=1;
r=pi/bet;
q=1/bet-1+delt;
lambd=(et+(1-alph)*((1+e*(r/(r-1))^(gamm-1))^(-1))*(thet/(thet-1)-delt*(alph/q))^(-1))/((1-alph)*z*((thet-1)/thet)^((1)/(1-alph))*(alph/q)^(alph/(1-alph)));
c=(1+e*((r/(r-1))^(gamm-1)))^(-1)*(1/lambd);
ksi=((thet-1)/thet)*lambd;
m=c*e*((r/(r-1))^(gamm));
y=c*(1-delt*(alph/q)*((thet-1))/thet)^(-1);
k=(alph/q)*((thet-1)/thet)*y;
in=delt*k;
h=(1/z)*(y/k^alph)^(1/(1-alph));
w=(1-alph)*((thet-1)/thet)*(y/h);
d=y-w*h-q*k;
n=y/h;
tau=pi;
c_obs=0;
m_obs=0;
in_obs=0;
pi_obs=log(pi);
r_obs=log(r)*400;
din=intrend;
dc=ctrend;
dpi=1*pi+constepi;
dr=1*r+conster;
dm=mtrend;
end;

steady;

resid;

estimated_params;
bet, 0.99;
gamm, 0.042;
phip, 50;
phik, 32.14;
rhoa, 0.75;
rhoe, 0.75;
rhox, 0.75;
rhoz, 0.75;
rhov, 0.5;
omegapi, 1.3;
omegay, 0.125;
omegatau, 0.38;
stderr epsilone, 1.39;
stderr epsilona, 0.01;
stderr epsilonx, 0.01;
stderr epsilonz, 0.01;
stderr epsilonv, 0.01;
end;

shocks;
var epsilona=sigmaa^2;
var epsilone=sigmae^2;
var epsilonx=sigmax^2;
var epsilonz=sigmaz^2;
var epsilonv=sigmav^2;
end;

% stoch_simul(graph_format = fig);

varobs in_obs c_obs m_obs r_obs pi_obs;

% identification;

estimation(optim=('MaxIter',200),datafile=Data,mode_compute=6,
first_obs=12,presample=4,lik_init=2,mode_check,prefilter=0,
mh_replic=0,mh_nblocks=2,mh_jscale=1.08,mh_drop=0.2,xls_range=C1:G118,graph_format=fig,prior_trunc=0,order=1);