var y c i k n
lb $\lambda$ 
a z w r
eta $\eta$ 
mu $\mu$
u
delta_u ${\delta(u)}$
x; 

varexo 
    eps_a_news ${\varepsilon_a^{news}}$
    eps_z_news ${\varepsilon_z^{news}}$

    eps_a_surprise ${\varepsilon_a^{surprise}}$
    eps_z_surprise ${\varepsilon_z^{surprise}}$;

parameters 
    beta $\beta$
    theta $\theta$
    sigma $\sigma$
    gamma $\gamma$
    alpha $\alpha$
    delta0 $\delta_0$
    delta1 $\delta_1$
    phi $\phi$
    psi $\psi$
    omega $\omega$ 
    chi $\chi$
    ra $\rho_a$
    rz $\rho_z$;

% set parameter values
beta   = 0.985;
theta  = 1.4;
sigma  = 1;
gamma  = 0.001;
alpha  = 0.64;
phi    = 1.5;
omega  = 1.3;
delta0 = 0.008;
delta1 = 0.08;
chi    = 1.009;
ra     = 0.95; // Persistence of Technology Shock
rz     = 0.95;

 
us  =  ((-delta1-delta1/phi)/(delta1/phi-(beta/chi)*delta1))^(1/phi);
delta_us = delta0+delta1/phi+(delta1*us)/phi;
rs  = delta1*us^phi;
ksy = (1-alpha)/rs;
isy = (chi-1+delta_us)*ksy;
csy = 1-isy;
ns  = 1/3;
ys  = ns*ksy^((1-alpha)/alpha)*us^(alpha/(1-alpha));
cs  = csy*ys;
is  = isy*ys;
ks  = ksy*ys;
ws  = (1-alpha)*ys/ns;
xs  = cs;
psis= ((beta*(1-gamma)-1)*ws)/(theta*ns^(theta-1)*cs*(beta*(1-gamma)-1)-ns^theta*gamma*ws);
lbs = ((cs-psis*ns^theta*cs)^(-sigma)*theta*psis*ns^(theta-1)*cs)/ws;
mus = ((cs-psis*ns^theta*cs)^(-sigma)*psis*ns^theta)/(beta*(1-gamma)-1);
etas= lbs;

model;
(c-(psi*(n^theta)*x))^(-sigma) + (mu*gamma*(x(-1)/(c))^(1-gamma)) = lb;

psi*(n^theta)*(c-(psi*(n^theta)*x))^(-sigma) + mu = beta*mu(+1)*(1-gamma)*((c(+1)/x)^gamma);

theta*psi*x*(n^(theta-1))*(c-(psi*(n^theta)*x))^(-sigma) = lb*w;

eta = (beta/chi)*(lb(+1)*r(+1) + eta(+1)*(1-delta_u(+1)));

lb*(1-alpha)*a*(u^(-alpha))*(k(-1))^(1-alpha)*(n^alpha) = eta*k*delta1*(u^(phi-1));

lb/z = eta*(1-(omega/2)*(i/i(-1)-1)^2-omega*(i/i(-1)-1)*(i/i(-1))) + beta*eta(+1)*omega*(i(+1)/i-1)*(i(+1)/i)^2;

x = (c^gamma)*x(-1)^(1-gamma);

y = c + i/z;

y = a*(u*k(-1))^(1-alpha)*(n^alpha);

chi*k = (1-delta_u)*k(-1) + i*(1 - (omega/2)*(i/i(-1)-1)^2);

delta_u = delta0+delta1/phi+(delta1*(u^phi))/phi;

w = alpha*a*y/n;

r = (1-alpha)*a*y/k(-1);

log(a)=ra*log(a(-1))+log(eps_a_surprise) + log(eps_a_news(-4));

log(z)=rz*log(z(-1))+log(eps_z_surprise) + log(eps_z_news(-4));
end;

initval;
a   = 1;
z   = 1;
y   = ys;
c   = cs;
i   = is;
k   = ks;
n   = ns;
lb  = lbs;
eta = etas;
mu  = mus;
w   = ws;
r   = rs;
x   = xs;
u   = us;
delta_u = delta_us;
end;


shocks;
    var eps_a_news=1;     //4 period anticipated TFP news shock
    var eps_a_surprise=1; //TFP surprise shock
    var eps_z_news=1;     //4 period anticipated TFP news shock
    var eps_z_surprise=1; //TFP surprise shock
end;

write_latex_static_model;
write_latex_dynamic_model;

steady;
check;

stoch_simul(order=1,irf=40);


//************ The following line generate the IRFs to a "pure" TFP news shock, i.e. the anticipated shock is 
//************ counteracted by an opposite surprise shock when it is supposed to realize

//initialize IRF generation

initial_condition_states = repmat(oo_.dr.ys,1,M_.maximum_lag);
shock_matrix = zeros(options_.irf,M_.exo_nbr); %create shock matrix with number of time periods in colums

// set shocks for pure news 
shock_matrix(1,strmatch('eps_z_news',M_.exo_names,'exact')) = 1; %set news shock to 1 (use any shock size you want)
shock_matrix(1+4,strmatch('eps_z_surprise',M_.exo_names,'exact')) = -1; %8 periods later use counteracting shock of -1

y2 = simult_(initial_condition_states,oo_.dr,shock_matrix,1);
y_IRF = y2(:,M_.maximum_lag+1:end)-repmat(oo_.dr.ys,1,options_.irf); %deviation from steady state


// manually select variables for figure
figure
subplot(2,3,1)
plot(y_IRF(strmatch('y',M_.endo_names,'exact'),:)); % use strmatch to select values
title('Output');
subplot(2,1,2)
plot(y_IRF(strmatch('n',M_.endo_names,'exact'),:));
title('Labor');
subplot(2,1,3)
plot(y_IRF(strmatch('x',M_.endo_names,'exact'),:));
title('X');
subplot(2,1,4)
plot(y_IRF(strmatch('i',M_.endo_names,'exact'),:));
title('Investment');
subplot(2,1,5)
plot(y_IRF(strmatch('k',M_.endo_names,'exact'),:));
title('Capital');
subplot(2,1,6)
plot(y_IRF(strmatch('c',M_.endo_names,'exact'),:));
title('Consumption');

// Automatically loop over variables for figure (may require different setting for subplots in larger models)
figure
for ii=1:M_.orig_endo_nbr
    subplot(3,3,ii)
    if max(abs(y_IRF(ii,:)))>1e-12 %get rid of numerical inaccuracies
        plot(y_IRF(ii,:));
    else
        plot(zeros(options_.irf,1));
    end
    title(deblank(M_.endo_names(ii,:)));
end