var r, pi_var, a, c, z, theta, y, w, q, x, g, sigma_z;   

varexo eps_z, eps_a, eps_r, eps_theta, eps_sigma_z;

parameters beta_par,pi_ss, z_ss, psi,a_ss, r_ss, eta, phi, g_ss, c_ss, y_ss, theta_ss,

 q_ss, omega, rho_pi, rho_g, rho_x, rho_a, rho_e rho_theta, x_ss,

rho_sigma_z;
beta_par=0.99;
psi=0.1;
omega=0.0617;

pi_ss=1.0086;
z_ss=1.0048;
g_ss=z_ss;
r_ss=pi_ss*z_ss/beta_par;
a_ss=1;
theta_ss=1.05;
x_ss=((theta_ss-1)/theta_ss)^omega;
y_ss=x_ss;
c_ss=y_ss;
q_ss=a_ss;

eta=1/omega;
phi=eta*(theta_ss-1)/psi;


rho_pi=0.3597;
rho_g=0.2536;
rho_x=0.0347;
rho_a=0.9470;
rho_e=0.9625;
rho_sigma_z=0.9;
rho_theta=rho_e;


model;

0=log(beta_par*r*a(+1)*c*z_ss)-log(pi_var(+1)*a*c(+1));

theta-1=theta*y^(eta-1)*(c/a)-phi*(pi_var/pi_ss-1)*pi_var/pi_ss+beta_par*phi*(a(+1)/a)*(c/c(+1))*(pi_var(+1)/pi_ss-1)*y(+1)*pi_var(+1)/(y*pi_ss);

(y^(eta-1))*c/a=w;

q=a^omega;

y=c+(phi/2)*(pi_var/pi_ss-1)^2;

x=y*a^(-omega);

log(z)=log(z_ss)+sigma_z*eps_z;

g=y*z/y(+1);

log(r)=log(r(-1))+rho_pi*log(pi_var/pi_ss)+rho_g*log(g/g_ss)+rho_x*log(x/x_ss)+eps_r;

log(theta)=(1-rho_theta)*log(theta_ss)+rho_theta*log(theta(-1))+eps_theta;

log(a)=rho_a*log(a(-1))+eps_a;

sigma_z=rho_sigma_z*sigma_z(-1)+eps_sigma_z(-1);

end;

initval;


r=r_ss;

pi_var=pi_ss;

a=1;

c=c_ss; 

z=z_ss;

theta=theta_ss;

g=g_ss;
 
y=y_ss;

w=(y_ss^eta);

q=1; 

x=x_ss;

sigma_z=0;

end;



shocks;

var eps_r=0.0031^2;

var eps_z=1;

var eps_theta=(0.0012*phi) ^2;

var eps_a=0.0405^2;

var eps_sigma_z=0.01^2;

end;



%% so far I just wrote the model.

steady;


stoch_simul(order=3, periods=5000, drop=2500, nograph);




ex_ = zeros(60,5);
ex_(2,5) = 1;

ergo_mean = mean(oo_.endo_simul,2);

% with stoch_simul we simulated the model to compute the ergodic mean, which is then fed into simult_
%simult_ simulates the model using a perturbation method

%oo_.endo_simul is a variable that stores the result of a simulation. Note that the variables are stored ROW-by-ROW
%Pruning=Discard higher order terms when iteratively computing simulations of the solution, as in Kim, Kim, Schaumburg and Sims (2008)

irfergo3 = simult_(ergo_mean,oo_.dr,ex_,3);

%     simult_('initial value' , 'structure where the reduced form solution of the model is stored' , 'matrix of innovations', 'order of the Taylor Approximations')  


%  Definition of the IRF as deviations from steady state
%      "       "  "   "   "   "         "    ergodic mean
%-------------------------------------





rep_ergomean=repmat(ergo_mean,1,size(irfergo2,2));
rep_steady=repmat(oo_.steady_state,1,size(irfergo2,2));

IRF3_ERGO=irfergo3-rep_ergomean;    %1st row is r, 2nd row pi_var....




titles={'Interest Rate','inflation', 'preference', 'consumption', 'TFP', 'theta', 'gdp', 'wages', 'natural output', 'output gap', 'growth', 'TFP uncertainty'};

figure;

for i=1:12
subplot(3,4,i)
plot(IRF3_ERGO(i,:),'r')
title(titles(i))
end