var x i r mpy pie z n dm;
varexo eps_mu eps_m;
parameters beta rho epsilon fi theta rhor taupi kappa rhom;

// Real
beta        = 0.99;
rho         = 1/beta-1;
epsilon     = 11;
fi          = 1;

// Nominal
theta       = 0.75;
rhor        = 0.7;
taupi       = 1.5;
lambda      = (1-beta*theta)*(1-theta)/theta;
kappa            = lambda*(1+fi);
rhom        = 0.5;


model;
x        = x(+1) - (r(0) - rho);
pie(0)  = beta*pie(+1)+ kappa*x(0);
n        = x;
r        = i - pie(+1);
i(0)    = rho + (1/4)*x(0)-(1/4)*mpy(0);
mpy(0)  = mpy(-1) - pie + dm(0) -z;
z        = eps_mu(0);
dm(0)   = rhom*dm(-1)+eps_m(0);
end;

shocks;
var        eps_mu;  stderr 1;
var        eps_m;  stderr 1;
end;
                
 
steady;
check;

stoch_simul(order=1,irf=30) pie x n i r z;
y_gap = x_eps_mu;
y     = y_gap + ones (30,1); 
n     = n_eps_mu;
 
data_m  = [ 4*pie_eps_m  x_eps_m  4*r_eps_m 4*i_eps_m ];
data_mu = [ 4*pie_eps_mu y_gap y n ];

Tbig=26;
tid = linspace(0,Tbig-1,Tbig)';          
zeroline = zeros(Tbig);             

figure;
subplot(4,1,1);
plot( tid,data_m(1:Tbig,1),'k-',...
      tid,zeroline,'k-');
title('Inflation');
subplot(4,1,2);
plot( tid,data_m(1:Tbig,2),'k-',...
      tid,zeroline,'k-');   
title('Output');
subplot(4,1,3);
plot( tid,data_m(1:Tbig,3),'k-',...
      tid,zeroline,'k-');   
title('Real rate');
subplot(4,1,4);
plot( tid,data_m(1:Tbig,4),'k-',...
      tid,zeroline,'k-');
title('Nominal rate');

figure;
subplot(4,1,1);
plot( tid,data_mu(1:Tbig,1),'k-',...
      tid,zeroline,'k-');
title('Inflation');
subplot(4,1,2);
plot( tid,data_mu(1:Tbig,2),'k-',...
      tid,zeroline,'k-');   
title('Output gap');
subplot(4,1,3);
plot( tid,data_mu(1:Tbig,3),'k-',...
      tid,zeroline,'k-');   
title('Output');
subplot(4,1,4);
plot( tid,data_mu(1:Tbig,4),'k-',...
      tid,zeroline,'k-');
title('Employment');