% Replication of Closing small open economy models
% Stephanie Schmitt-Grohe , Martin Uribe 
% Journal of International Economics 61 (2003) 163–185


% Endogenous Discount Factor Model without Internalization


%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

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

var y c k h d invest beta tb ca lambda a eta;
varexo ea;

parameters alpha delta gamma phi psi1 omega r rhoa sigma_a;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

gamma = 2; 	      %intertemporal elasticity of substitution
omega = 1.455; 	  %exponent of labor in utility function
delta = 0.1; 	  %Depreciation rate
alpha = 0.32; 	  %Capital elasticity of the production function
psi1 = 0.11135;	  %0.11; %Minus elasticity of discount factor with respect to its argument
rhoa = 0.42; 	  %Serial correlation of productivity shock
phi = 0.028; 	  %Parameter of adjustment cost function
r = 0.04; 	      %world interest rate
sigma_a = 0.0129; %STD of innovation in transitory technology shock

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Steady-Sate Values of Variables  %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Steady-state value of hours(h):
hss = ((1 - alpha)*(alpha/(r+delta))^(alpha/(1-alpha)) )^(1/(omega-1) );

% Steady-state value of consumption:
css = (1+r)^(1/psi1) + (hss^omega)/omega - 1;

% Steady-state value of capital:
kss = hss/( ((r + delta)/alpha)^(1/(1-alpha)) ); 

% Steady-state value of investment:
iss = delta*kss;

% Steady-state value of output:
yss = (kss^alpha)*(hss^(1-alpha));

// steady state trade balance
tbss        = yss-css-iss;

// steady state financial assets 		
dss         = tbss/r;
 
%U = (( exp(c) - exp(h)^omega/omega-gamma)-1)  / (1-gamma);         %Utility Function
%Uc = ((exp(c) - exp(h)^omega/omega)^(-gamma));                     %marginal utility of consumption
%Uh = ((exp(c) - exp(h)^omega/omega)^(-gamma))*(exp(h)^(omega-1));  %marginal utility of hours.

%betac = -psi1*( (1+ exp(c) - ((exp(h)^omega) /omega ))^ ( -psi1-1));
%betah = -psi1*( (1+ exp(c) - ((exp(h)^omega) /omega ))^ ( -psi1-1))*(exp(h)^(omega-1));


%----------------------------------------------------------------
% 3. Model  %%%  
%----------------------------------------------------------------

model; 



%%FOCs:
%U = (( exp(c) - exp(h)^omega/omega-gamma)-1)  / (1-gamma); 
%Uh = ((exp(c) - exp(h)^omega/omega)^(-gamma))*(exp(h)^(omega-1));
%betah = -psi1*( (1+ exp(c) - ((exp(h)^omega) /omega ))^ ( -psi1-1))*(exp(h)^(omega-1));

% 1. Equation: Subjective discount factor
beta =  ( (1+ exp(c) - ((exp(h)^omega) /omega ))^ ( -psi1));

% 2. Equation: market clearing condition
d  = (1+r)*d(-1)- exp(y)+exp(c)+exp(invest)+((phi/2)*(exp(k)-exp(k(-1)))^2);

% 3. Equation: Production Function
exp(y) = exp(a)*(exp(k(-1))^alpha)*(exp(h)^(1-alpha));

% 4. Equation: Evolution of capital stock                                          
exp(invest) = exp(k) - (1-delta)*exp(k(-1)); 

% 5. Equation: Euler Equation
exp(lambda) = beta*exp(lambda(+1))*(1+r);  

% 6. Equation: Marginal Utility of Consumption
((exp(c) - exp(h)^omega/omega)^(-gamma)) - exp(eta)*(-psi1*( (1+ exp(c) - ((exp(h)^omega) /omega ))^ ( -psi1-1))) = exp(lambda);

% 7. Equation
exp(eta) = -( (( exp(c(+1)) - exp(h(+1))^omega/omega-gamma)-1)  / (1-gamma) ) + exp(eta(+1))*( (1+ exp(c(+1)) - ((exp(h(+1))^omega) /omega ))^ ( -psi1));

% 8. Equation: Marginal Utility of Hours worked                                    
((exp(c) - exp(h)^omega/omega)^(-gamma))*(exp(h)^(omega-1))+ exp(eta)*(-psi1*( (1+ exp(c) - ((exp(h)^omega) /omega ))^ ( -psi1-1))*(exp(h)^(omega-1))) = exp(lambda)*(1-alpha)*( exp(a)*(exp(k(-1))^alpha)*(exp(h)^(-alpha)) ); 

% 9. Equation: Marginal returns to capital
exp(lambda)*(1 + phi*(exp(k) - exp(k(-1)))) = beta*exp(lambda(+1))*(alpha*exp(y(+1))/exp(k(+1))+1-delta+phi*(exp(invest(+1))-delta*exp(k)));

% 10. Equation: Trade balance
tb = 1-((exp(c)+exp(invest))/exp(y));

% 11. Equation:
ca     = (1/exp(y))*(d-d(-1));

% 12. Equation:
a = rhoa*a(-1) + ea;
 
end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval;

beta  = 1/(1+r);
a     = 0;
ea     = 0;

h     = log(hss);
k     = log(kss);
y     = log(yss);
c     = log(css);
invest     = log(iss);
d     = (exp(y)-exp(c)-exp(invest))/r;
tb    = 1-((exp(c)+exp(invest))/exp(y));
ca    = 0;

%U = (( css - hss^omega/omega-gamma)-1)  / (1-gamma);      %Utility Function
%Uc = ((css - hss^omega/omega)^(-gamma));                  %marginal utility of consumption
%Uh = ( (css - hss^omega/omega)^(-gamma))*(hss^(omega-1));  %marginal utility of hours.

%betac = -psi1*( (1+ css - ((hss^omega) /omega ))^ ( -psi1-1));
%betah = -psi1*( (1+ css - ((hss^omega) /omega ))^ ( -psi1-1))*(hss^(omega-1));

lambda= log((exp(c)-((exp(h)^omega)/omega))^(-gamma) / (-( ((( css - hss^omega/omega-gamma)-1)  / (1-gamma)) / (1-beta) ))* (-psi1*( (1+ css - ((hss^omega) /omega ))^ ( -psi1-1))));
eta = -( ((( css - hss^omega/omega-gamma)-1)  / (1-gamma)) / (1-beta) );

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%specify the innovations%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shocks;
var ea; stderr sigma_a; %comment Philipp: you can just specify stderr instead of variance. 
end;

%check;
resid(1); 
steady;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%simulation%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

stoch_simul(order=1, periods=1000, irf=40) y c h k d invest tb ca a;