%%DYNARE CODES FOR MASTER


var k, z, c, Y, r, rtilde, R, q, w, wtilde, b, taok, z_idf, phi_A, phi_B, PHI_A, PHI_B, P_A, P_B, v_D,psi,theta,miu,xsi;

varexo epis;

parameters alpha, delta, beta, rho, sigma, m_g, taon, eta, lambda;

beta=0.98;

alpha=0.32;

delta=0.1;

rho=0.95;

sigma=0.007;

m_g=0.2;               %    Ratio of government expenditure to GDP

taon = 0.3;

eta=0;  

lambda=2;

model;

r = alpha*z*k(-1)^(alpha-1);

rtilde(-1) = (1-taok(-1))*alpha*k(-1)^(alpha-1);

R = 1+z*rtilde(-1)-delta;

Y =z*k(-1)^alpha;

w = (1-alpha)*Y;

wtilde = (1-taon)*w;

z_idf=(1/q-1+delta)/rtilde;
phi_A = normpdf((log(z_idf)-rho*log(z))/sigma);

PHI_A = 0.5*(1+erf((log(z_idf)-rho*log(z))/sigma));

P_A = 1+(lambda-1)* PHI_A;

phi_B = normpdf((log(z_idf)-rho*log(z)-sigma^2)/sigma);

PHI_B = 0.5*(1+erf((log(z_idf)-rho*log(z)-sigma^2)/sigma));

P_B = 1+(lambda-1)* PHI_B;

v_D = (1-delta-1/q)*P_A+z^rho*exp(sigma^2/2)*rtilde*P_B;

c+k+m_g*Y=Y+(1-delta)*k(-1);

m_g*Y+b(-1)+wtilde+z*rtilde(-1)*k(-1)=Y+q*b;


R(+1)/c(+1)+eta*v_D=1/(beta*c);

q/c=beta/c(+1);

(miu*R(+1)+xsi*c-miu(+1))/(c(+1)^2)+1/c(+1)=theta(+1)+beta*xsi(+1)/c(+2);

theta/beta= theta(+1)*(r(+1)+1-delta)+psi(+1)*(r(+1)-z(+1)*rtilde)+eta*v_D;

((lambda-1)*((1-delta-1/q)*phi_A+z^rho*exp(sigma^2/2)*rtilde*phi_B)/(sigma*rtilde)-z^rho*exp(sigma^2/2)*P_B)*eta*(k-miu)=z(+1)*(-miu/c(+1)-psi(+1)*k);

psi*q=beta*psi(+1);

psi*b+xsi=((lambda-1)*((1-delta-1/q)*phi_A+z^rho*exp(sigma^2/2)*rtilde*phi_B)/(sigma*(1/q-1+delta))-P_A)*eta*beta*(k-miu)/q^2;

log(z)=rho*log(z(-1))+sigma*epis;

end;

initval;
taok = 0;
k = 4.20975;
c =  0.846245;
r = alpha*k^(alpha-1);
rtilde = (1-taok)*r;
R=1+rtilde-delta;
q = beta;
Y = k^alpha;
w = (1-alpha)*Y;
wtilde = (1-taon)*w;
b = (taon*w+taok*r*k-m_g*Y)/(1-q);
z=1;
z_idf=(1/q-1+delta)/rtilde;
phi_A = normpdf((log(z_idf)-rho*log(z))/sigma);

PHI_A = 0.5*(1+erf((log(z_idf)-rho*log(z))/sigma));

P_A = 1+(lambda-1)* PHI_A;

phi_B = normpdf((log(z_idf)-rho*log(z)-sigma^2)/sigma);

PHI_B = 0.5*(1+erf((log(z_idf)-rho*log(z)-sigma^2)/sigma));

P_B = 1+(lambda-1)* PHI_B;

v_D = (1-delta-1/q)*P_A+z^rho*exp(sigma^2/2)*rtilde*P_B;

psi = 4.33885;
theta = 0.708712;
miu = 3.11257;
xsi = 2;
end;
steady;

shocks;
var epis = 1;
end;


check;
stoch_simul (drop = 100, periods = 1000, irf=0, order=1);