var c k y b h g y_h i w r c_y log_y log_h log_y_h log_w log_c log_c_y ;

varexo eps_g eps_b;

parameters delta psi alpha beta mu_g mu_b rho_g rho_b ;

% Parameter Values 

beta = 0.95; 
psi = 0.33; 
alpha = 0.68; 
delta = 0.07; 
mu_g = log(1.0066);
mu_b = log(1.0066);
rho_g = 0.1;
rho_b = 0.1;


Model; 

y = k(-1)^(1-alpha)*exp(g+b)^(alpha-1)*h^alpha; 
log(g) = (1-rho_g)*log(mu_g)+rho_g*log(g(-1))+eps_g;  
log(b) = (1-rho_b)*log(mu_b)+rho_b*log(b(-1))+eps_b;  
r = (1-alpha)*(y/k(-1))*exp(g+b); 
w= alpha*(y/h);
c^(-1)= exp(g(+1)+b(+1))^(-1)*beta*c(+1)^(-1)*(1+r(+1)-delta); 
h^(1/psi) = c^(-1)*w;
y_h = y/h; 
c + k = y + (1-delta)*k(-1)*exp(g+b)^(-1); 
i = y-c;
c_y=c/y;

// use logarithm to get variables in percentage deviations
log_y=log(y);
log_h=log(h);
log_y_h=log(y_h);
log_w=log(w);
log_c=log(c);
log_c_y=log(c_y);
end; 


steady_state_model;
g = mu_g;
b= mu_b;
r=(1/(beta*exp(mu_g)^(-1)*exp(mu_b)^(-1))-(1-delta));
y_k=r/((1-alpha)*exp(mu_g)*exp(mu_b));
k_y=(1-alpha)/r;
i_y=(1/y_k)*(1-(1-delta)*exp(mu_g)^(-1)*exp(mu_b)^(-1));
c_y=1-i_y;
h= (alpha*1/c_y)^(psi/(1+psi));
k=(((exp(mu_g)^(alpha-1)*exp(mu_g)^(alpha-1)*(h^alpha))/y_k))^(1/alpha);
y=k^(1-alpha)*h^alpha*exp(mu_g)^(alpha-1)*exp(mu_g)^(alpha-1);
y_h =y/h;
c=c_y*y;
i=i_y*y;
w=alpha*(y/h);
log_y=log(y);
log_h=log(h);
log_y_h=log(y_h);
log_w=log(w);
log_c=log(c);
log_c_y=log(c_y);
end;


shocks;
var eps_g; stderr 0.01;
var eps_b; stderr 0.01;
end;
resid(1);
steady;
check;


stoch_simul(order=1, periods=1000, nofunctions) y h y_h c_y;

// Rebuild non-stationary time series by remultiplying with A_{t} and B_{t} 

log_Gamma_0=0; //Initialize Level of Technology at t=0;
log_B_0=0; //Initialize Level of Technology at t=0;
log_Gamma(1,1)=g(1,1)+log_Gamma_0; //Level of Tech. after shock in period 1
log_B(1,1)=b(1,1)+log_B_0; //Level of Tech. after shock in period 1
// reaccumulate the non-stationary level series  (non-stationary log-level variables)

for ii=1:options_.periods

    log_Gamma(ii+1,1)=g(ii,1)+log_Gamma(ii,1);
    log_B(ii+1,1)=b(ii,1)+log_B(ii,1);
    log_y_nonstationary(ii+1,1)=log_y(ii,1)+log_Gamma(ii,1)+log_B(ii,1);                
    log_h_nonstationary(ii+1,1)=log_h(ii,1)+log_B(ii,1);  
    log_y_h_nonstationary(ii+1,1)=log_y_h(ii,1)+log_Gamma(ii,1);
  
end