Hi, I am trying to simulate a DSGE Model with ZLB using Dynare 5.0 and Occbin block. My question is how to do stochastic simulation and get for example the unconditional mean of production function (y) or Welfare function (U);

```
var c r k l w pi r_k y mc i z U;
varexo eps_z;
parameters varphi beta delta alpha phi epsilon rho_r phi_pi phi_y rho_z sigma;
varphi = 1.5;
beta = 0.99;
delta = 0.025;
alpha = 0.33;
phi = 106;
epsilon = 6;
rho_r = 0.6;
phi_pi = 1.5;
phi_y = 0.125;
rho_z = 0.95;
sigma = 0.01;
model;
#z_ss = 0;
#pi_ss = 0;
#r_ss = 1/beta;
#l_ss = 1/3;
#r_k_ss = 1/beta - (1 - delta);
#mc_ss = (epsilon-1)/epsilon;
#w_ss = (mc_ss*(alpha/r_k_ss)^alpha)^(1/(1-alpha))*(1-alpha);
#k_ss = w_ss/r_k_ss*alpha/(1-alpha)*l_ss;
#y_ss = exp(z_ss)*k_ss^alpha*l_ss^(1-alpha);
#i_ss = delta*k_ss;
#c_ss = y_ss - i_ss;
#chi = w_ss/c_ss*l_ss^(-varphi);
//Welfare
U = ln(c) - l^(1+varphi)/(1+varphi) + beta*U(+1);
//Euler equation
beta*c/c(+1)*r/(1+pi(+1)) = 1;
//Labor supply
chi*l^varphi = w/c;
//Capital supply
beta*c/c(+1)*(r_k(+1) + 1 - delta) = 1;
//Capital law motion
k = (1-delta)*k(-1) + i;
//Production function
y = exp(z)*k(-1)^alpha*l^(1-alpha);
//Marginal cost
mc = (w/(1-alpha))^(1-alpha)*(r_k/alpha)^alpha;
//Labor capital ratio
w/r_k = (1-alpha)/alpha*k(-1)/l;
//Phillips curve
phi*pi*(1+pi) = 1 +epsilon*(mc-1) + beta*(c/c(+1)*phi*pi(+1)*(1+pi(+1))*y(+1)/y);
//Resource constraint
(1-phi/2*pi^2)*y = c + i;
[name='Taylor rule',relax = 'ZLB']
r = r_ss + phi_pi*pi + phi_y*(y-y_ss);
[name='Taylor rule',bind = 'ZLB']
r = 1;
//TFP AR(1)
z = rho_z*z(-1) +eps_z;
end;
occbin_constraints;
name 'ZLB'; bind r<1;
end;
steady_state_model;
z = 0;
pi = 0;
r = 1/beta;
l = 1/3;
r_k = 1/beta - (1 - delta);
mc = (epsilon-1)/epsilon;
w = (mc*(alpha/r_k)^alpha)^(1/(1-alpha))*(1-alpha);
k = w/r_k*alpha/(1-alpha)*l;
y = exp(z)*k^alpha*l^(1-alpha);
i = delta*k;
c = y - i;
U = (ln(c) - l^(1+varphi)/(1+varphi))/(1-beta);
end;
//----------------------------------------
// Shocks
//----------------------------------------
shocks;
var eps_z; stderr 0.1;
end;
steady;
occbin_setup(simul_maxit = 100, simul_periodic_solution, simul_debug);
```