% February 18, 2016

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

var   c_h c_f c_ht c_ft c_hh c_hf b_h c_fh c_ff b_f p_h p_f p_ht p_ft et w r_h r_f y_h y_f y_hn y_fn l_h f g pi_f pi_g ;

parameters betaa gammaa theta omega lambda alphaa gammaF g_ss sigma_g;
varexo eps_g;

betaa  =  0.99; % discount factor
gammaa =  0.4; % weight on tradable consumption
theta  =  3  ; % CES between tradable and non-tradable
omega  =  0.5; % home-bias in tradable consumption
lambda =  1.5; % CES between H and F tradable consumption
alphaa =  0.7; % exponent on production
gammaF =  0.1; % FIs credit constraint tightness
g_ss   =  0.005; % exogenous growth rate
sigma_g = 0.05; % standard deviation of g shock

%----------------------------------------------------------------
% 2. Model
%----------------------------------------------------------------

% equations to characterize equilibrium

model;

c_ht = ( omega^(1/lambda)*c_hh^((lambda-1)/lambda) + (1-omega)^(1/lambda)*c_hf^((lambda-1)/lambda) )^(lambda/(lambda-1));
c_ft = ( (1-omega)^(1/lambda)*c_fh^((lambda-1)/lambda) + omega^(1/lambda)*c_ff^((lambda-1)/lambda) )^(lambda/(lambda-1));

c_h  = ( gammaa^(1/theta)*c_ht^((theta-1)/theta) + (1-gammaa)^(1/theta)*y_hn^((theta-1)/theta) )^(theta/(theta-1));
c_f  = ( gammaa^(1/theta)*c_ft^((theta-1)/theta) + (1-gammaa)^(1/theta)*y_fn^((theta-1)/theta) )^(theta/(theta-1));

p_ht = ( omega*p_h^(1-lambda) + (1-omega)*p_f^(1-lambda) )^(1/(1-lambda));
p_ft = ( (1-omega)*p_h^(1-lambda) + omega*p_f^(1-lambda) )^(1/(1-lambda));

c_ht/y_hn = gammaa/(1-gammaa) * (p_ht)^(-theta);
c_hh/c_hf = omega/(1-omega)   * (p_h/p_f)^(-lambda);
(1/c_h)*(c_h/y_hn)^(1/theta) = betaa*(1+r_h) * (1/(c_h(+1)*(1+g_ss)))*((c_h(+1)*(1+g_ss))/y_hn(+1))^(1/theta);
p_h*c_hh + p_f*c_hf + y_hn + b_h = p_h*y_h + y_hn + (1+r_h(-1))/(1+g_ss)*b_h(-1) + et*pi_g;

% (10,11)

c_ft/y_fn = gammaa/(1-gammaa) * (p_ft/et)^(-theta);
c_fh/c_ff = (1-omega)/omega   * (p_h/p_f)^(-lambda);
(1/c_f)*(c_f/y_fn)^(1/theta) = betaa*(1+r_f) *(1/(c_f(+1)*(1+g_ss)))*((c_f(+1)*(1+g_ss))/y_fn)^(1/theta);
p_h*c_fh + p_f*c_ff + et*y_fn + et*b_f = p_f*y_f + et*y_fn + et*(1+r_f(-1))/(1+g_ss)*b_f(-1) + et*pi_f;

y_h  = l_h^alphaa;
alphaa * l_h^(alphaa-1) = w;
y_hn = (1-l_h)^alphaa;
alphaa * (1-l_h)^(alphaa-1) = p_h*w;

y_f  = 1;
y_fn = 1;

% (20,21)

f = 1/gammaF * (1 - (1+r_h)/(1+r_f)*(et/et(+1)));
g = 0.8 * g(-1) + eps_g;

pi_f = ((1+r_f(-1)) - (1+r_h(-1))*et(-1)/et) * f(-1)/(1+g_ss);
pi_g = ((1+r_f(-1)) - (1+r_h(-1))*et(-1)/et) * g(-1)/(1+g_ss);

y_h = c_hh+c_fh;
y_f = c_hf+c_ff;
b_f + f + g = 0;

% b_h + et*b_f = 0;

end;

%----------------------------------------------------------------
% 3. Computation
%----------------------------------------------------------------

% Instead of initval, steady state variables are calculated by m-file

shocks;
var eps_g = sigma_g^2;
end;

steady;
check;

%stoch_simul(order=2,nograph,nomoments,nocorr);
%stoch_simul(order=2,irf=50,irf_shocks=(eps_r));
stoch_simul(order=2);