close all

var
    u us      	
    v vs      	
    ev evs     	
    vpg vpgs        
    m ms      	
    y ys      	
    c cs      	
    i is      	
    k ks      	
    l ls      	
    w ws      	
	z zs      	
    q
;

varexo
	epsilon
    epsilons
;

parameters
    upsilon
    beta
    gamma
    psi
	alpha
    delta
    big_omega   // Auxiliary variable to compute the steady state
    big_phi     // Auxiliary variable to compute the steady state      
    core_e
    rho
;



% calibrate parameters
upsilon     = 0.34;
beta        = 0.99;
delta       = 0.025;
gamma       = 10;
psi         = 0.5;
alpha		= 0.36;
big_omega   = ((1/alpha) * ((1/beta)-1+delta)) ^ (1/(alpha-1));
big_phi     = upsilon/(1-upsilon)*(1-alpha) * (big_omega ^ alpha);
rho			= 0.906;
core_e		= 0.258;
zeta        = 0;


% Some values for the steady state
z_ss        = 0;
l_ss 		= big_phi / (big_omega^alpha - delta*big_omega + big_phi);
k_ss 		= l_ss * big_omega;
y_ss        = exp(z_ss)*(k_ss^alpha)*(l_ss^(1-alpha));
c_ss 		= l_ss * (big_omega^alpha - delta*big_omega);
i_ss        = delta * k_ss;
u_ss        = c_ss^upsilon * (1-l_ss)^(1-upsilon);
v_ss        = u_ss;
ev_ss       = v_ss^(1-gamma);
w_ss        = (1-alpha) * (k_ss^alpha) * (l_ss^(-alpha));
r_ss        = alpha * (k_ss^(alpha-1)) * (l_ss^(1-alpha));
lambda_ss   = (v_ss^psi) * (1-beta) * u_ss^(1-psi) * upsilon * c_ss^(-1);
q_ss        = 1;
m_ss        = beta;


model;
% recursive preferences
    v = ( (1-beta)*(u)^(1-psi) + beta*(ev)^((1-psi)/(1-gamma))) ^ (1/(1-psi));
    vs = ( (1-beta)*(us)^(1-psi) + beta*(evs)^((1-psi)/(1-gamma))) ^ (1/(1-psi));

    ev = vpg(+1);
    evs = vpgs(+1);

    vpg = (v)^(1-gamma);
    vpgs = (vs)^(1-gamma);

% period utility
	u = c^upsilon * (1-l)^(1-upsilon);
	us = cs^upsilon * (1-ls)^(1-upsilon);

% stochastic discount factor
    m = beta * (u(+1)/u)^(1-psi) * (c/c(+1)) * (v(+1)/ev^(1/(1-gamma)))^(psi-gamma);
	ms = beta * (us(+1)/us)^(1-psi) * (cs/cs(+1)) * (vs(+1)/evs^(1/(1-gamma)))^(psi-gamma);

    y = exp(z) * ((k(-1))^alpha) * (l^(1-alpha));
    ys = exp(zs) * ((ks(-1))^alpha) * (ls^(1-alpha));

% technology
    z = rho*z(-1) + epsilon;
    zs = rho*zs(-1) + epsilons;

% LOM capital
    k = (1-delta)*(k(-1)) + i;
    ks = (1-delta)*(ks(-1)) + is;

% wage
    w = (1-alpha) * exp(z) * ((k(-1))^alpha) * (l ^(-alpha));
    ws = (1-alpha) * exp(zs) * ((ks(-1))^alpha) * (ls ^(-alpha));

	upsilon * c^(-1) * w = (1-upsilon) * (1-l)^(-1);
	upsilon * cs^(-1) * ws = (1-upsilon) * (1-ls)^(-1);

% euler eq. capital
    m * (alpha * exp(z(+1)) * (k^(alpha-1)) * (l(+1))^(1-alpha) + 1 - delta) = 1;
    ms * (alpha * exp(zs(+1)) * (ks^(alpha-1)) * (ls(+1))^(1-alpha) + 1 - delta) = 1;
 
% resource constraint
    c+i = y;
	cs + is = ys;

% exchange rate
    
    q(+1)/q = ms/m;
    
    // COMMENT OUT THIS LINE FOR AUTARKY
    //q=1;
end;


initval;
    u 		= u_ss;
	us 		= u_ss;
    v 		= v_ss;
    vs 		= v_ss;
    vpg     = v_ss^(1-gamma);
    vpgs    = v_ss^(1-gamma);
    ev 		= ev_ss;
    evs 	= ev_ss;
    k 		= k_ss;
    ks 		= k_ss;
    l 		= l_ss;
    ls 		= l_ss;
    y 		= y_ss;
	ys 		= y_ss;
    c 		= c_ss;
    cs 		= c_ss;
    i 		= i_ss;
    is 		= i_ss;
    m 		= m_ss;
    ms 		= m_ss;
    w 		= w_ss;
    ws 		= w_ss;
    z 		= z_ss;
    zs 		= z_ss;
    q       = 1;
end;


shocks;
  	var epsilon = (0.00852)^2;
  	var epsilons = (0.00852)^2;
end;


check;
steady;
model_diagnostics;

stoch_simul(order = 1, hp_filter=1600, irf=00);