//Comment1: TFP has permanent shock Z(t) and transi shock A_T, and z(t) is defined as the log difference of Z(t) //
//Comment2: Almost all varibles are in log term, except z, a, mu_p, i,A_T //
//Comment3: output, wage, wage_star are de-trended by TFP Z //

var 
y i u z A_T mu_p 
f w_p pi pi_star n_d n_s g h 
m re_y re_d p_y p_d;

varexo e_u e_z e_A_T;

parameters 
beta psi omega 
rho_i i_bar phi_y phi_pi rho_u 
rho_z z_bar rho_A_T  
theta_p alpha_p 
y_SS i_SS u_SS z_SS A_T_SS mu_p_SS 
f_SS w_p_SS pi_SS pi_star_SS n_d_SS n_s_SS g_SS h_SS 
m_SS re_y_SS re_d_SS p_y_SS p_d_SS;

psi = 2;                // Risk aversion parameter
omega=0.35;             // Inverse of Frisch labor elasticity 0.35


rho_z = 0.391;           // Autocorrelation of Permanent Tech shock
z_bar = 0.0047;           // long run growth rate of productivity
rho_A_T = 0.95;

rho_i = 0.63;           // Interest rate smoothing coefficient in policy rule
i_bar = 0.029;          // average interest rate equal steady state -log (pricing kernal)
phi_y = 0.125;          // Response to output gap
phi_pi = 1.5;           // Response to inflation
rho_u = 0.564;          // Autocorrelation of policy shock

beta = exp(psi*(z_bar) - i_bar);             // Discount factor: i_bar = -log(beta)+psi*(z_bar)

theta_p = 6;            // Elasticity of substitution of differentiated goods
alpha_p = 0.63;         // Price rigidity


// The following are steady state value for endogenuous variables//

i_SS = i_bar;
u_SS = 0;
z_SS = z_bar;
A_T_SS = 1;
mu_p_SS = theta_p/(theta_p-1);

w_p_SS = log(A_T_SS*(theta_p-1)/theta_p);
f_SS = 0;
pi_SS = 0;
pi_star_SS = 0;
g_SS = log(1/(1-alpha_p*beta*exp((1-psi)*z_SS)));
h_SS = log(1/(1-alpha_p*beta*exp((1-psi)*z_SS)));

y_SS = (1/(omega+psi))*log(exp(w_p_SS)*exp(omega*(log(A_T_SS))));
n_d_SS = y_SS - log(A_T_SS);
n_s_SS = n_d_SS;

m_SS = -i_bar;
re_y_SS = -m_SS;
re_d_SS = -m_SS;
p_y_SS = log(exp(z_SS)/(exp(re_y_SS)-exp(z_SS)));
p_d_SS = log(exp(z_SS)/(exp(re_d_SS)-exp(z_SS)));



model;

// Wage Setting //
exp(w_p) = exp(omega*n_s)*exp(psi*y);

// Price Dispersion //
exp(f) = (1-alpha_p)*exp(-theta_p*pi_star) + alpha_p*exp(theta_p*pi)*exp(f(-1));

// Price Setting //
exp(pi_star)*exp(h)*A_T = (theta_p/(theta_p-1))*exp(w_p)*exp(g);
exp(h) = 1 + alpha_p*exp(m(+1))*exp((-theta_p-1)*pi(+1))*exp(y(+1)-y)*exp(z(+1))*exp(h(+1));
exp(g) = 1 + alpha_p*exp(m(+1))*exp(theta_p*pi(+1))*exp(y(+1)-y)*exp(z(+1))*exp(w_p(+1)-w_p)*exp(g(+1));

// Price Aggregators //
1 = (1-alpha_p)*exp((1-theta_p)*pi_star) + alpha_p*exp((theta_p-1)*pi);

// Aggregate Labor Supply and Demand //
exp(n_s) = exp(n_d);
exp(n_d) = exp(y)*exp(f)/A_T;

// Markups //
mu_p = A_T*exp(-f)*exp(-w_p);                //NOTICE: Markups are NOT defined in log term!!!!!

// Interest rate rule //
i = rho_i*i(-1) + (1-rho_i)*(i_bar + phi_pi*pi + phi_y*(y-y_SS)) + u;
u = rho_u*u(-1) + e_u;

i = -log(exp(m(+1)-pi(+1)));

// TFP //
z = (1-rho_z)*z_bar + rho_z*z(-1) + e_z;
log(A_T) = rho_A_T*log(A_T(-1)) + e_A_T;

// Pricing Kernal, Return, and Price-Cash ratio //
//exp(m) = beta^((1-gamma)/(1-psi))*exp(-(y-y(-1))*psi*(1-gamma)/(1-psi))*exp(-z*psi*(1-gamma)/(1-psi))*exp(re_q*(psi-gamma)/(1-psi));
//exp(re_q) = ((1+exp(p_q))/exp(p_q(-1)))*exp(z)*(exp(y) - ((1-psi)/(1+omega))*exp(w_p)*exp(n_s))/(exp(y(-1)) - ((1-psi)/(1+omega))*exp(w_p(-1))*exp(n_s(-1)));
//1 = exp(m(+1)+re_q(+1));
exp(m) = beta*exp(-psi*(y-y(-1)))*exp(-z*psi);

exp(re_y) = exp(z)*(1+exp(p_y))*exp(y)/(exp(p_y(-1)+y(-1)));
1 = exp(m(+1)+re_y(+1));

exp(re_d) = exp(z)*(1-1/mu_p)*(1+exp(p_d))*exp(y)/((1-1/mu_p(-1))*exp(p_d(-1)+y(-1)));
1 = exp(m(+1)+re_d(+1));

end;

initval;
i = i_SS;
u = u_SS;
z = z_SS;
A_T = A_T_SS;
y = y_SS;
mu_p = mu_p_SS;

f = f_SS;
w_p = w_p_SS;
pi = pi_SS;
pi_star = pi_star_SS;
n_d = n_d_SS;
n_s = n_s_SS;
g = g_SS;
h = h_SS;

m = m_SS;
re_y = re_y_SS;
re_d = re_d_SS;
p_y = p_y_SS;
p_d = p_d_SS;
end;

resid(1);
steady;
check;

shocks ;
var e_A_T;
stderr 0.001;
var e_z;
stderr 0.002;
var e_u;
stderr 0.0015;
end ;

stoch_simul(irf =15, order=2);