% SOE-RBC model for emerging -country business cycles
% This is the same model applied in chapter 5 - Open Economy
% Macroeconomics, 2014 - Uribe & Schmitt-Grohe

% University of Brasilia - UnB
% Department of Economics 
% Macroeconomics III
% Professor Joaquim Andrade
% Student Raphael Barcelos

%=========================
% 0. Preparation
%=========================

close all;

%=========================
% 1. Variables
%=========================

var
y               % Gross Domestic Product
c               % Consumption
h               % Worked hours
lambda          % Lagrange multiplier
k               % stock of physical capital
d               % stock of debt
i               % investiment
r               % domestic interest rate
vu              % preference shocks
a               % TFP
g               % gross growth of X
S               % aggregate shifts in domestic absorption
s               % Proportional S
X               % Stochastic trend (Public spending shocks)
mu              % Stochastic term on interest-rate
;

varexo
epsm            % Stochastic term on interest-rate
epsa            % Stochastic term on TFP
epsv            % Stochastic term on preferences
epsg            % Stochastic term on X
epss            % Stochastic term on public spending            
;

parameters
    beta        % subjective discount factor
    omega       % labor supply elasticity
    gama        % inverse of intertemporal elasticity of substitution
    alpha       % capital elasticity of output
    eta         % fraction of wage bill subject to working-capital constraint 
    g
    delta       % depreciation rate
    phi         % capital adjustment cost parameter
    rstar       % international interest-rate
    psi         % debt adjustment cost parameter
    y_bar       
    d_bar       
    s_bar
    rhoa
    rhom
    rhov
    rhog
;

%=========================
% 2. Calibration
%=========================

    beta  = 0.9286;     % subjective discount factor
    omega = 1.6000;     % labor supply elasticity
    gama  = 2.0000;     % inverse of intertemporal elasticity of substitution
    alpha = 0.3200;     % capital elasticity of output
    eta   = ;     % fraction of wage bill subject to working-capital constraint 
    g_bar = 1.0107;
    delta = 0.1255;     % depreciation rate
    phi   = ;     % capital adjustment cost parameter
    rstar = 0.1000;     % international interest-rate
    psi   = ;     % debt adjustment cost parameter
    y_bar = ;     
    d_bar = ;     
    s_bar = ;
    rhoa  = ;
    rhom  = ;
    rhov  = ;
    rhog  = ;
    rhos  = ;
    s_bar = ;

% Steady-state equations (calculated)

r_ss       = rstar;
g_ss       = g_bar;
h_ss       = g_bar*((g^gama/beta-1+delta)*(1/(alpha)))^(1/(alpha))* ((1+eta*r_ss/(1+r_ss))*(1/((1-alpha)*g_bar^(1-alpha)))^(-1/alpha);
k_ss       = (h^(alpha+omega-1)/((1-alpha)*g_bar^(1-alpha))*(1+eta*r_ss/(1+r_ss)))^(1/alpha);
i_ss       = delta*k_ss;
y_ss       = (k_ss^alpha)*(h_ss^(1-alpha));
d_ss       = d_bar;
c_ss       = y_ss-i_ss-s_ss-d_ss*(r_ss/(1+r_ss));

%=========================
% 3. Model
%=========================

model;
% Eq. 1
exp(vu)*(exp(c)-omega^(-1)*exp(h)^omega)^(-gama)=exp(lambda);
% Eq. 2
exp(h)^(omega-1)=(1-alpha)*exp(a)*exp(g)^(1-alpha)*(exp(k)/exp(h))^alpha*(1+eta*exp(r)/(1+exp(r)))^(-1);
% Eq. 3
exp(lambda)=(beta/exp(g)^gama)*(1+exp(r))*exp(lambda(+1));
% Eq. 4
(1+phi*((exp(k(+1))/exp(k))*exp(g)-g_bar))*exp(lambda)=(beta/exp(g)^gama)*exp(lambda(+1))*(1-delta+alpha*exp(a(+1))*(exp(g(+1))*exp(h(+1))/exp(k(+1)))^(1-alpha) 
    +phi*(((1-delta)*exp(k(+1))+exp(i(+1)))/exp(k(+1)))*exp(g(+1))*((((1-delta)*exp(k(+1))+exp(i(+1)))/exp(k(+1)))*exp(g(+1))-g_bar)-(phi/2)*((((1-delta)*exp(k(+1))+exp(i(+1)))/exp(k(+1)))*exp(g(+1))-g_bar)^2);
% Eq. 5
(d(+1)/(1+exp(r)))*exp(g)=d-exp(y)+exp(c)+exp(s)+exp(i)+(phi/2)*((exp(k(+1))/exp(k))*exp(g)-g_bar)^2*exp(k); 
% Eq. 6
exp(r)=rstar+psi*(exp((d(+1)-d_bar)/y_bar)-1)+exp(mu-1)-1;
% Eq. 7
exp(k(+1))*exp(g)=(1-delta)*exp(k)+exp(i);
% Eq. 8
exp(y)=exp(a)*exp(k)^alpha(exp(g)*exp(h))^(1-alpha);
% Eq. 9
a(+1) = rhoa*a + epsa;
% Eq. 10
mu(+1) = rhom*mu + epsm;
% Eq. 11
vu(+1) = rhov*vu + epsv;
% Eq. 12
g(+1)/g_bar = rhog*g/g_bar + epsg;
% Eq. 13
exp(s) = exp(S)/exp(X(-1));
% Eq. 14
exp(g) = exp(X)/exp(X(-1));
% Eq. 15
exp(s(+1))/s_bar = rhos*exp(s) + epss;
end;

%=========================
% 4. Initial values
%=========================

initval;
r = log(g_bar^(gama)/beta - 1);
d = d_ss;
h = log(h_ss);
k = log(k_ss);
y = log(y_ss);
c = log(c_ss);
i = log(i_ss);


lambda = ;
end;
%=========================
% 5. Calculations
%=========================