var ci cs ki ks bi bs q r mu w rb a kappa_t BETAi BETAs ;

varexo nui nuk;

parameters EPSILON  KAPPA_bar OMEGA N SIGMA RHO_a RHO_kappa SIGMANUi SIGMANUk ETA ZETAi ri ZETAs Z1 ; 
%ZETAs 
parameters ci_s cs_s ki_s ks_s bi_s bs_s q_s r_s mu_s w_s rb_s a_s kappa_t_s BETAi_s BETAs_s Z1_s ZETAs_s;


ETA       = 0.0125;

EPSILON   = 0.39;

ri = 0.01;     % the world interest rate

OMEGA     = 0.105;

KAPPA_bar = .75; 

N         = 0.5;

SIGMA     = 2; 

RHO_a     = 0.9;

RHO_kappa = 0.9; 

SIGMANUi  = 0.006;

SIGMANUk  = 0.011;

corr_nui_nuk = 0.00;

model;

%1
BETAi =  ZETAi*(1+ N*exp(ci))^(-ETA);
%2  
BETAs = ZETAs*(1+ (1-N)*exp(cs))^(-ETA);


%3
exp(ci)+(exp(q)*exp(ki)) = exp(w)+(exp(q)+exp(r))*exp(ki(-1))+ bi-exp(rb(-1))*bi(-1); 



%4
bi = exp(kappa_t)*exp(q)*exp(ki); 


%5 
exp(cs) + exp(q)*exp(ks) = exp(w)+ exp(q)*exp(ks(-1))+ bs - exp(rb(-1))*bs(-1) + Z1*(exp(ks(-1)))^OMEGA; 


%6
exp(ci)^(-SIGMA)= BETAi *exp(ci(+1))^(-SIGMA)*((exp(q(+1))+exp(r(+1)))/exp(q)) + exp(kappa_t)*mu; 


%7
exp(ci)^(-SIGMA) = BETAi * exp(ci(+1))^(-SIGMA)*exp(rb) + mu; 


%8
exp(cs)^(-SIGMA)= BETAs *exp(cs(+1))^(-SIGMA)*((exp(q(+1))+ OMEGA*Z1*(exp(ks))^(OMEGA-1))/exp(q));  


%exp(cs)^(-SIGMA) = BETAs * exp(cs(+1))^(-SIGMA)*exp(rb); 


%9
exp(r)= exp(a)*(1-EPSILON)*N^(1-EPSILON)*(N*exp(ki(-1)))^(-EPSILON); 


%10
exp(w)= (exp(a)*(N*exp(ki(-1)))^(1-EPSILON))- exp(r)*N*exp(ki(-1));  


%11
N*exp(ki)+(1-N)*exp(ks) = 1; 

%12
N*bi + (1-N)*bs = 0; 



%13
 N*exp(ci) + (1-N)*exp(cs) = exp(a)*(N*exp(ki(-1)))^(1-EPSILON)+(1-N)*Z1*exp(ks(-1))^OMEGA; 

%14
a = RHO_a*a(-1) + nui;

%15
exp(kappa_t) = (1-RHO_kappa)*KAPPA_bar + RHO_kappa*exp(kappa_t(-1)) + nuk;
end;

initval;
BETAi = BETAi_s;
BETAs = BETAs_s;
ci= log(ci_s);
cs = log(cs_s);
ki = log(ki_s);
ks = log(ks_s);
r = log(r_s);
rb = log(rb_s);
w = log(w_s);
%a= log(a_s);
a = 0;
%kappa_t = log(kappa_t_s);
kappa_t = log(KAPPA_bar);
bi = bi_s;
bs = bs_s;
%the value of bs_s is negative, therefore no log()
%mu = log(mu_s);
mu = mu_s;
q = log(q_s);
nui = 0;
nuk = 0;
end;

steady;
resid (1);
check(qz_zero_threshold= 1e-16);
shocks;
var nui =SIGMANUi^2;
var nuk = SIGMANUk^2;
end;
stoch_simul(order =1,irf=40, hp_filter=1600);