var  c_h, c_f, c, lambda, p, p_h, p_f, w, l, ni, y_h, b,bstar,s, pi_h,c_hstar, p_hstar,m_c, m_d,m,tax,mc,profit,mm,
nistar,pstar,ystar,zeta,g;

varexo epsil_zeta, epsil_pstar, epsil_nistar,  epsil_ystar, epsil_m, epsil_g,epsil_mm;

parameters rho, beta, hb, sigma_l, alpha, gamma,mu_h,phim, thetah, chi,omega,alphay,alphai,delta,nibar,nistarbar,
rhonistar,rhopstar,rhoystar,rhozeta,rhog;


%------------------------------------------------------------------------------------------%
%                                 PARAMETERS                                               % 
%------------------------------------------------------------------------------------------%

rho=2;
beta=0.98;
hb=0.6;
sigma_l=1;
alpha=0.35;
gamma=0.5089;
thetah=0.5;
phim=0.5;
chi=1.03;
omega=0.5;
nibar=1/beta;
nistarbar=1.02;
alphay=0.6;
alphai=-0.5;
delta=0.27;
mu_h=1.2;

rhonistar=0.68;
rhozeta=0.85;
rhopstar=0.84;
rhoystar=0.88;
rhog=0.78;



model(linear);
%1
c_h=-rho*(p_h-p)+c;
%2
c_f=-rho*(p_f-p)+c;
%3
p+c=0.5089*(p_h+c_h)+(1-0.5089)*(p_f+c_f);
%4
mc=w+alpha/(1-alpha)*y_h-1/(1-alpha)*zeta;
%5
p=p_h/(1+(1-gamma)/gamma*2^(1-rho))+p_f/(1+gamma/(1-gamma)*(1/2)^(1-rho));
%6
% p+c=(l+w)*0.918-(b(+1)-nibar*b-nibar*ni(-1))*0.001/(nibar-1)+(profit*0.05)-(p+tax)*0.05-(m_d(+1)-nibar*m_d-nibar*ni(-1))*0.0008/(nibar-1)-(bstar(+1)+s-nistarbar*bstar-nistarbar*nistar(-1)-nistarbar*s)*0.0834/(nistarbar-1);
c+p=2*(m_c-m_c(+1))+0.6*l+0.6*w-0.3/(nibar-1)*(m_d(+1)+b(+1))+(profit*0.05)-(p+tax)*0.05+0.3/(nibar-1)*nibar*((m_d+b)+ni(-1))-0.1/(nistarbar-1)*(bstar(+1)+s)+0.1/(nistarbar-1)*nistarbar*(bstar+nistar(-1)+s);
%7
lambda=hb/(1-beta*hb)/(1-hb)*c(-1)+1/(1-hb)*c+beta*hb/(1-hb)/(1-beta*hb)*c(+1);
%8
lambda+w-p=sigma_l*l;
%9
lambda-p+p(+1)=0;
%10
% lambda+s-p=lambda(+1)+s(+1)+nistar-p(+1);
lambda+s+p(+1)-p-nistar-lambda(+1)-s(+1)=(0.05*bstar(+1)+0.05*s-0.05*p-0.05*y_h)*(2+1/(1/2+0.05));
%11
lambda-p=lambda(+1)+ni-p(+1);
%12
p_f=s+pstar;
%13
p_hstar=p_h-s;
%14
c_hstar=-chi*omega*(p_hstar-pstar)+omega*ystar+(1-omega)*c_hstar(-1);
%15
pi_h=p_h-p_h(-1);
%16
pi_h=1/(1+beta)*pi_h(-1)+beta/(1+beta)*pi_h(+1)+(1-thetah)*(1-beta*thetah)/(1+beta)/thetah*mc;
%17
m=0.0987*m_c+(1-0.0987)*m_d;
%18
m=m_c+mm; 
%19
mm=alphay*y_h+alphai*ni+epsil_mm; 
%20
m_c=delta*y_h(-1)+epsil_m;
%21
tax+p=36*(m-m(+1))+p+g+nibar*32*(m_d(+1)+ni)-32*m_d;
%22
y_h=0.1887*c_h+(1-0.1887-0.171)*c_hstar+0.171*g;
%23
y_h=zeta+(1-alpha)*l;
%24
profit=y_h+p_h/(1-1/mu_h)-mc/(mu_h-1);

% Shock Equations
nistar=rhonistar*(nistar(-1))+epsil_nistar;

zeta=rhozeta*(zeta(-1))+epsil_zeta;

ystar=rhoystar*(ystar(-1))+epsil_ystar;

pstar=rhopstar*(pstar(-1))+epsil_pstar;

g=rhog*(g(-1))+epsil_g;


end;


resid(1);
steady;
check;

shocks; 

var epsil_zeta; stderr 0.1; 
var epsil_nistar; stderr 0.1;
var epsil_ystar; stderr 0.1; 
var epsil_m; stderr 0.1; 
var epsil_mm; stderr 0.1;
var epsil_pstar; stderr 0.1;
var epsil_g; stderr 0.1;

end;


stoch_simul(linear,irf=100) y_h c;