%NK with sticky price(about subsidy)


var c i pi l k lamda w mc a y vp pisharp x1 x2 r tau;
varexo ea eg;
parameters eta sigma beta sita delta alfa phi epsilon rhoa rhog;
parameters cs is pis ls ks lamdas ws mcs as ys vps pisharps x1s x2s rs taus ;

beta =.9;
sigma = 1;
eta = 1;
epsilon =2;
rhoa = .95;
rhog = 0;
alfa=0.3;
delta=0.1;
%sita is subsity index
sita=0.2;
%the stickiness parameter
phi = .75;

%steady state calculation

as = 1;
%zero inflation steady state
pis = 1;
taus=1;
pisharps=((pis^(1-epsilon) - phi)/(1-phi))^(1/(1-epsilon));
vps =(1-phi)*(pis/pisharps)^epsilon/(1- pis^epsilon*phi);
mcs = (1-phi*beta*pis^epsilon)/(1-phi*beta*pis^(epsilon-1))
           *pis/pisharps*(epsilon-1)/epsilon;
rs=1/beta-1;
ws=(mcs*(alfa^alfa)*((1-alfa)^(1-alfa))*(rs^(-alfa)))^(1/(1-alfa));
ls = ((1-taus)*(as/vps*(ws*alfa/((1-alfa)*rs))^alfa-delta*ws*alfa/((1-alfa)*rs))*ws)^(1/(eta+sigma));
ks=ls*ws*alfa/(rs*(1-alfa));
ys = as*(ks^alfa)*(ls^(1-alfa))/vps;
is=delta*ks;
cs= ys-is;
x1s = cs^(-sigma)*mcs*ys/(1-beta*phi*pis^epsilon);
x2s =cs^(-sigma)*ys/(1-beta*phi*pis^(epsilon-1));
lamdas=cs^(-sigma)/(1-beta*sita);

model;
%(1) home Euler equation
c^(-sigma) = lamda-beta*sita * lamda(+1);

%(2) labor supply
l^eta =(1-tau)*c^(-sigma)*w;

%(3)capital suply
k(+1)=i+(1-delta)*k;

%(4)lamda 
lamda=beta*lamda(+1)*(1+r)*pi;

%(5) labor demand
w =mc*(1-alfa)*a*(k/l)^alfa;

%(6) capital-labor
r/w = (alfa/(1-alfa))*(l/k);



%(7) accounting identity
c+i=y;

%(8) the production technology
y = a *k^alfa*l^(1-alfa)/vp;

%(9) the price dispersion
vp = (1-phi)*pisharp^(-epsilon)*(pi^epsilon)
                + (pi^epsilon) *phi*vp(-1);

%(10) inflation evolution
pi^(1-epsilon) = (1-phi)*pisharp^(1-epsilon) + phi;

%(11)the sticky price equation
pisharp = epsilon/(epsilon -1 )*pi*x1/x2;

%(12) the auxiliary x1
x1 =c^(-sigma)*y*mc+phi*beta*(pi(+1))^epsilon*(x1(+1));

%(13) the auxiliary x2
x2 =c^(-sigma)*y +phi*beta*(pi(+1))^(epsilon-1)*(x2(+1));

%(14) technology shock
log(a) = rhoa*log(a(-1)) + ea;

%(15)government balance
tau*w*l=sita*(c(-1));

%(16)tax shock
log(tau)=(1-rhog)*log(taus)+rhog*log(tau(-1))+eg;

end;
initval;
c = cs;
i = is;
pi=pis;
l=ls;
k=ks;
w = ws;
tau=taus;
mc = mcs;
a = as;
y = ys;
vp = vps;
pisharp = pisharps;
x1 = x1s;
x2 = x2s;
lamda=lamdas;
r=rs;
end;

shocks;
var ea = .01^2;
var eg =.01^2;
end;

resid(1);
steady;
check;
%when used in loops, noprint is a good option.
stoch_simul(order=1) i pi l k r mc a y ;