var y c k l z r w i; 
varexo e a;


parameters beta psi sigma delta alpha rhoz gammax n x i_y k_y;

% set parameter values
sigma=1;                % risk aversion
alpha= 0.33;            % capital share
i_y=0.25;               % investment-output ration
k_y=10.4;               % capital-output ratio
x=0.0055;               % technology growth (per capita output growth)
n=0.0027;               % population growth
rhoz=0.97;              %technology autocorrelation base on linearly detrended Solow residual

model;
exp(c)^(-sigma)=beta/gammax*exp(c(+1))^(-sigma)* (alpha*exp(z(+1))*(exp(k)/exp(l(+1)))^(alpha-1)+(1-delta));
psi*exp(c)^sigma*1/(1-exp(l))=exp(w);
gammax*exp(k)=(1-delta)*exp(k(-1))+exp(i);
exp(y)=exp(i)+exp(c);
exp(y)=exp(z)*exp(k(-1))^alpha*exp(l)^(1-alpha);
exp(w)=(1-alpha)*exp(y)/exp(l);
r=4*alpha*exp(y)/exp(k(-1));
z=rhoz*z(-1)+e + a;
end;

steady_state_model; 
    gammax=(1+n)*(1+x); 
    delta=i_y/k_y-x-n-n*x;  
    beta=gammax/(alpha/k_y+(1-delta));
    l_ss=0.33;
    k_ss=((1/beta*gammax-(1-delta))/alpha)^(1/(alpha-1))*l_ss;
    i_ss=(x+n+delta+n*x)*k_ss;
    y_ss=k_ss^alpha*l_ss^(1-alpha);
    c_ss=y_ss-i_ss;
    psi=(1-alpha)*(k_ss/l_ss)^alpha*(1-l_ss)/c_ss^sigma; 
    i=log(i_ss);
    w=log((1-alpha)*y_ss/l_ss);
    r=4*alpha*y_ss/k_ss;
    y=log(y_ss);
    k=log(k_ss);
    c=log(c_ss);
    l=log(l_ss);
    z=0;
end;

shocks;
    var e=1;
    var a=1; 
end;

steady;
check;

stoch_simul(order=1,irf=20)y z;
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:2,strmatch('e',M_.exo_names,'exact'),strmatch('a',M_.exo_names,'exact')) = 2; %set news shock to 1 (use any shock size you want)
shock_matrix(1,strmatch('a',M_.exo_names,'exact')) = 1; %8 periods later use counteracting shock of -1
y2 = simult_(initial_condition_states,oo_.dr,shock_matrix,2);
y_IRF = y2(:,M_.maximum_lag+1:end)-repmat(oo_.dr.ys,2,options_.irf); %deviation from steady state

