Sorry, I was thinking about IRFS, not simulations. I guess you are looking for something like
[code]var c k y B h A y_h i w r g_A g_B;
varexo eps_a eps_b;
parameters delta psi alpha beta ;
% Parameter Values
beta = 0.95;
psi = 0.33;
alpha = 0.68;
delta = 0.05;
Model;
y = k(-1)^(1-alpha)*exp(eps_a+eps_b)^(alpha-1)h^alpha;
g_A=eps_a;
g_B=eps_b;
log(A)=log(A(-1)) + eps_a;
log(B)=log(B(-1)) + eps_b;
r = (1-alpha)(y/k(-1))exp(eps_a+eps_b);
w= alphay/h;
c^(-1)= exp(eps_a(+1)eps_b(+1))^(-1)betac(+1)^(-1)(1+r(+1)-delta);
h^(1/psi)*B^(-1) = c^(-1)*w;
y_h = y/h;
c + k = y + (1-delta)*k(-1)*exp(eps_a+eps_b)^(-1);
i = y-c;
end;
steady_state_model;
A=1;
B=1;
r=(1/beta)-(1-delta);
y_k=r/(1-alpha);
k_y=(1-alpha)/r;
i_y=(1/y_k)delta;
c_y=1-i_y;
h= (alpha1/c_y)^(psi/(1+psi));
k=((h^alpha)/y_k)^(1/alpha);
y=k^(1-alpha)h^alpha;
c=c_yy;
i=i_yy;
w=alpha(y/h);
y_h =y/h;
end;
shocks;
var eps_a; stderr 1;
var eps_b; stderr 1;
end;
resid(1);
steady;
check;
stoch_simul(order=1,periods=1000,drop=960) y h c w y_h i g_A g_B;
// Rebuild non-stationary time series by remultiplying with A_{t} and B_{t}
log_A_0=0; //Initialize Level of Technology at t=0;
log_B_0=0; //Initialize Level of Technology at t=0;
log_A(1,1)=log_A_0 ; //Level of Tech. after shock in period 1
log_B(1,1)=log_B_0; //Level of Tech. after shock in period 1
// reaccumulate the non-stationary level series
for ii=1:options_.periods
log_A(ii+1,1)=log_A(ii,1)+g_A(ii,1);
log_B(ii+1,1)=log_B(ii,1)+g_B(ii,1);;
y_nonstationary(ii+1,1)=y(ii,1)+log_A(ii,1)+log_B(ii,1);
h_nonstationary(ii+1,1)=h(ii,1)+log_B(ii,1);
c_nonstationary(ii+1,1)=c(ii,1)+log_A(ii,1)+log_B(ii,1);
w_nonstationary(ii+1,1)=w(ii,1)+log_A(ii,1)+log_B(ii,1);
y_h_nonstationary(ii+1,1)=y_h(ii,1)+log_A(ii,1)+log_B(ii,1);
end
[/code]