%=========================================================================%
% 
% 
%=========================================================================%


var  d, c, c_s, c_w, h, h_w, h_s, y, i, k, a, tb, ca, riskpremium, r, v, s, q, w ;  

varexo e;               


parameters  gamma, omega, rho, sigmae, delta, psi, alpha, phi, beta, r_w, ksi, chi_w, chi_s, mu, d_bar ;


ksi=0.5;  // (ratio of participant)
chi_w=0.6; // (weight for labor elasticity)
chi_s=0.5; // (weight for labor elasticity)
alpha  = 0.3;  // (production funtion power) 
rho    = 0.42;   //  (ar1 tech shock) 
phi    = 0.028;  //(adjustment cost) 
r_w    = 0.04;   //(world interest rate)
gamma  = 2;  // (degree of risk aversion)
omega  = 2;  //  (labor elasticity)
psi    = 0.000742;  //  (weight for risk premium)
delta  = 0.025;  //  (depriciation rate)
sigmae = 0.0129; // 
beta   = 1/(1+r_w);
h_w_ss =((1-alpha)/(chi_w*omega))^(1/(omega-1))*((1/alpha)*(1/beta-1+delta))^(alpha/((alpha-1)*(omega-1)));
h_s_ss=((1-alpha)/(chi_s*omega))^(1/(omega-1))*((1/alpha)*(1/beta-1+delta))^(alpha/((alpha-1)*(omega-1)));
w_ss=(1-alpha)*((1/alpha)*(1/beta-1+delta))^(alpha/(alpha-1));
h_ss=ksi*h_s_ss+(1-ksi)*h_w_ss; // (steady state for labor)
k_ss=((1/alpha)*(1/beta-1+delta))^(1/(alpha-1))*h_ss;  //(steady state for capital)
i_ss   = delta*k_ss;  // (steady state for investment)
y_ss   = (k_ss^alpha)*(h_ss^(1-alpha));  // (steady state for GDP)
mu=0.3;
d_bar=mu*k_ss; //(steady state for foreign debt)
d_ss   = d_bar;   // (steady state for foreign debt)
v_ss=y_ss-w_ss*h_ss-i_ss-r_w*d_ss;
c_w_ss =(1-alpha)* ((1-alpha)/(chi_w*omega))^(1/(omega-1))*((1/alpha)*(1/beta-1+delta))^((alpha*omega)/((alpha-1)*(omega-1)));
c_s_ss = (1-alpha)* ((1-alpha)/(chi_w*omega))^(1/(omega-1))*((1/alpha)*(1/beta-1+delta))^((alpha*omega)/((alpha-1)*(omega-1)))+v_ss;
c_ss=ksi*c_s_ss+(1-ksi)*c_w_ss; // (steady state for consumption)
tb_ss  = y_ss-c_ss-i_ss;// (steady state for trade balance)
q_ss=1;
s_ss=(beta/(1-beta))*v_ss;


model;

//shareholders
q=1;
w*h_s+q*(exp(s)+exp(v))=c_s+q(+1)*exp(s(+1));

w=chi_s*omega*(exp(h_s)^(omega-1));

exp(s(+1))/(exp(s(+1))+exp(v(+1)))=beta*(exp(c_s(+1))-chi_s*(exp(h_s(+1))^omega))/((exp(c_s)-chi_s*(exp(h_s)^omega))^(-gamma));


//workers

w*exp(h_w)=c_w;
w=chi_w*omega*(exp(h_w)^(omega-1));


//aggregation
exp(c)=ksi*exp(c_s)+(1-ksi)*exp(c_w);
exp(h)=ksi*exp(h_s)+(1-ksi)*exp(h_w);


//firms
    exp(y) = exp(a)*(exp(k(-1))^alpha)*(exp(h)^(1-alpha)); // production function

    exp(k(+1)) = exp(i)+(1-delta)*exp(k)-(phi/2)*(exp(k(+1))-exp(k))^2; // capital law of motion

exp(v(+1))=exp(y)-w*exp(h)-exp(i)+d(+1)-(1+exp(r))*d; // dividend


(1/beta)=(1+exp(r(+1)));

1+phi*(exp(k(+1))-exp(k))=beta*(alpha*exp(a)*(exp(k(-1))^alpha-1)*(exp(h)^(1-alpha)))+1-delta+phi*(exp(k(+2))-exp(k(+1)));

w=(1-alpha)*exp(a)*(exp(k(-1))^alpha)*(exp(h)^(-alpha)); //



    a = rho*a(-1)+e; // technology shock

    tb = 1-((exp(c)+exp(i))/exp(y)); // trade balance

    ca = (1/exp(y))*(d-d(-1));  // current account      


// risk premium

 riskpremium = psi*(exp(d-d_bar)-1);
    exp(r) = r_w+riskpremium;




end;


initval;
    r     = log((1-beta)/beta);
    q     = 1;
    d     = d_ss;
    w     = log(w_ss);
    h     = log(h_ss);
    h_w     = log(h_w_ss);
    h_s     = log(h_s_ss);
    k     = log(k_ss);
    y     = log(y_ss);
    c_s     = log(c_s_ss);
    c_w     = log(c_w_ss);
    i     = log(i_ss);
    v     = log(v_ss);
    s     = log(s_ss);
    tb    = 1-((exp(c)+exp(i))/exp(y));

end;


check;
steady; 


shocks;
    var e; stderr sigmae;
end;


stoch_simul(order=1, irf=24, nograph);