% Small-Open Economy Model - Real Business Cycle, 2014
% Base model from Chapter 4 - 2014, Uribe, Martin & Stephanie Schmitt-Grohe
% Simulation to Peruvian Economy, from 1960-2013
% This archive intents inference Bayesian Estimation
% This is a Dynare archive.

%??????????????????????????????????
% 0. Preamble
%??????????????????????????????????

close all;

%??????????????????????????????????
% 1. Variables and parameters
%??????????????????????????????????

% Variables generated by Dynare
var
y				% Gross Domestic Product
c				% Household consumption
h 				% Worked hours
k				% Capital stock
d				% Household?s debts
i				% Investment
tb				% Trade balance
ca 				% Current account
a				% Total Productivity Factors
r				% Domestic interest rate
p				% Risk premium
lambda			% Lagrange multiplier

% Variables observable 
y_obs			% GDP observable
%c_obs			% Household?s consumption observable
%i_obs			% Investment observable
%tby_obs			% Trade balance observable
%r_obs			% Interest rate observable
;

% Stochastic variables
varexo
e				% Productivity shock
;

parameters
betta				% inter temporal discount factor
gamma			% risk aversion degree
ommega			% wage elasticity
alppha			% capital participation in GDP
deltta			% capital depreciation rate
db  			% d bar
rhho				% TFP persistence shock
phhi				% capital adjustment cost coefficient
pssi				% risk premium coefficient
r_star			% internacional interest rate
sigmae
;

%??????????????????????????????????
% 2. Calibration
%??????????????????????????????????

% calibrated parameters
r_star = .05;
betta = 1/(1+r_star);
gamma = 2;
alppha = .41;
deltta = .11;
sigmae=.0129;

% Posteriors average
 db = 1.3973;
 ommega = 2.2367;
 rhho = .99;
 phhi = 2.3907; 
 pssi = .3622;

h_ss = ((1-alppha)*(alppha/(r_star+deltta))^(alppha/(1-alppha)))^(1/(ommega-1)); 
k_ss   = h_ss/(((r_star+deltta)/alppha)^(1/(1-alppha)));
i_ss   = deltta*k_ss;                                                     
y_ss   = (k_ss^alppha)*(h_ss^(1-alppha));                                   
%d_bar  = 0.7442;
d_ss   = db;                                                        
c_ss   = y_ss-i_ss-r_star*d_ss;
tb_ss  = y_ss-c_ss-i_ss;

%??????????????????????????????????
% 3. Model
%??????????????????????????????????

model; 
% eq. 4.15
    d = (1+exp(r(-1)))*d(-1)- exp(y)+exp(c)+exp(i)+(phhi/2)*(exp(k)-exp(k(-1)))^2;
% eq. product function
    exp(y) = exp(a)*(exp(k(-1))^alppha)*(exp(h)^(1-alppha));
% eq. capital?s law of movement
    exp(k) = exp(i)+(1-deltta)*exp(k(-1));
% eq. 4.7
    exp(lambda)= betta*(1+exp(r))*exp(lambda(+1)); 
% eq. 4.8
    (exp(c)-((exp(h)^ommega)/ommega))^(-gamma) = exp(lambda);
% eq. 4.9
    ((exp(c)-((exp(h)^ommega)/ommega))^(-gamma))*(exp(h)^ommega) = exp(lambda)*(1-alppha)*exp(y); 
% eq. 4.10, replacing the capital's law of motion k(+2) - k(+1) 
    exp(lambda)*(1+phhi*(exp(k)-exp(k(-1)))) = betta*exp(lambda(+1))*(alppha*exp(y(+1))/exp(k)+1-deltta+phhi*(exp(i(+1))-deltta*exp(k))); 
% eq. 4.12 TFP
    a = rhho*a(-1)+e; 
% eq. 4.19
    tb = 1-((exp(c)+exp(i))/exp(y));
% eq. 4.20
    ca = (1/exp(y))*(d-d(-1));      
% eq. functional form of risk premium                           
    p = pssi*(exp(d-db)-1);
% eq. 4.13
    exp(r) = r_star+p;
% Equations from observable variables
    y_obs=log(y);
    %c_obs=log(c);
    %i_obs=log(i);
    %tby_obs=tb;
    %r_obs=r;

end;

initval;
    r     = log((1-betta)/betta);
    d     = d_ss;
    h     = log(h_ss);
    k     = log(k_ss);
    y     = log(y_ss);
    c     = log(c_ss);
    i     = log(i_ss);
    tb    = 1-((exp(c)+exp(i))/exp(y));
    lambda= log((exp(c)-((exp(h)^ommega)/ommega))^(-gamma));
end;

shocks;
    var e; stderr sigmae;
end;
steady;

estimated_params;
db,gamma_pdf,1.3973,.7; 
ommega,gamma_pdf,2.2367,.7;
rhho,gamma_pdf,.99,.7;
phhi,gamma_pdf,2.3907,.7; 
pssi,gamma_pdf,.3622,.7;
end;

varobs y_obs; %c_obs i_obs tby_obs r_obs;

estimated_params_init(use_calibration);
 db,1.3973;
 ommega,2.2367;
 rhho,.99;
 phhi,2.3907; 
 pssi,.3622;
end;

estimation(datafile=peru,nobs=54,mh_replic=2000,mh_jscale=0.2);