% este es el model del paper de 2003 incluyendo su funci—n de costos, 
%incluyendo theta del paper del 2009 y los par‡metros

%AKI YA DIFERENCIASTE ENTRE R (QUE ES LA TASA DE RETORNO DEL CAPITAL) E I (QUE
% ES EL INTERES QUE PAGAS POR LA DEUDA MAS EL INTERES MUNDIAL)

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

close all;

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


var  y, c, h, d, x, k, a, mc,  tby,  i, tb, cay, i_debt, r, g,w;   

varexo e_a,e_g;                                    
                                             
parameters  gamma, omega, psi,alpha, phi,i_star,delta, rho_a, sigma_a, beta,  d_bar, rho_g, tao_c,sigma_g;
 

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

gamma  = 2;
omega  = 1.45; %exponent of labor in utility function
psi    = 0.000742;
alpha  = 0.32;
rho_a= 0.42;
rho_g= 0.42;
tao_c= .14;
phi    = 0.028;
i_star = 0.04;
delta  = .1;	
sigma_a = .0129;
sigma_g = .0129;
d_bar  = 0.7442; %foregin debt
			
        	
%Steady states of the variables

d_barra   = d_bar;   %From first order condition wrt to consumption , substituing b(1+r)=1


h_barra   = (((1-alpha)/(1+tao_c))*(alpha / (i_star +delta)) ^ (alpha/(1-alpha)))^ (1/(omega-1)); %From consumption/leisure equation

k_barra   = h_barra *( ( alpha/(i_star + delta) )^(1 / (1-alpha) )  ); %From euler equation 2

x_barra   = delta*k_barra;      %From capital accumulation equation  

y_barra   = (k_barra^alpha)*(h_barra^(1-alpha));      %From Production function    
                                             
c_barra   = (y_barra- x_barra)/(1+tao_c);  
%From resource constraint, sustituyendo inversion, trade balance
% ahi deber’a de haber una g_barra pero esta toma el valor de 1 por eso desaparece sino seria 
% y-x-g-id=c

tb_barra  = (d_barra*i_star); % Foreign net debt position

tby_barra= tb_barra/y_barra; 

cay_barra= (-i_star*d_barra + tb_barra)/y_barra ; %Current account relative to output que es cero

%set the subjective discount factor equal to the world interest rate
beta   = 1/(1+i_star);	

%g_barra=tao_c*c_barra+d_barra-i_star*d_barra; pero yo creo que mas bien es como una variable ex—gena como 
% la aÉ que en el steady state es uno



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

model;



%%%%1.- Marginal utility of consumption 
	mc=  ((c-1/omega*(h^omega))^(-gamma))/(1+tao_c);


%%%%2.- Consumption/leisure trade off

(h^(omega-1))  = w*((1-alpha)/(1+tao_c));

%2.1 Marginal product of labor

w=a*(k(-1)^alpha)*(h^(-alpha));


%3.- marginal product of capital
r=a*alpha*(k(-1)^(alpha-1))*(h^(1-alpha));

%4.-Euler equation with capital cost adjustment 
mc*(1+phi*(k -k(-1))) = beta*(mc(+1))*((r(+1))+1-delta+phi*(k(+1)-k)); 


%5.-Production Function
   
   	 y = a*(k(-1)^alpha)*(h^(1-alpha));

%6.-Capital evolution

  	 k = x+(1-delta)*k(-1); 

%7.-Rate that depends on foreign debt level
 	i_debt= psi*(exp(d-d_bar)-1);

%8.-interest rate 
   	i = i_star +  i_debt;

%9.-technology process
	log(a) = rho_a*log(a(-1))+e_a; 

%%%%10 Y 11.- aggregate resource constraint show trade balance: y-c-x-g (incluyendo costos capital adj)= nx
	tb= y -c - x -((phi/2)*(k-k(-1))^2)-g;

	tb= ((1+i)*d(-1))-d;

%%%%%12.- Government constraint

	g+(1+i)*d(-1)= d + tao_c*c;

%%%%%13.-Government shock process

	log(g) = rho_g*log(g(-1))+e_g;

% 14.- net exports divided by the output
    	tby = (tb)/y;

%15.- current account
  	cay = (-i*d(-1))/y+ (tby);


end;

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

initval;

%For each of your variables you need an initial value (take the one from the steady state)
% y, c, h, d, x, k, a, mc,  tby,  i, tb, cay, i_debt;   

y     = y_barra;
c     = c_barra;
h     = h_barra;
d     = d_barra;
x     = x_barra;
k     = k_barra;


cay= cay_barra; %equal to zero in steady state
tby= tby_barra;


mc=  ((c-1/omega*(h^omega))^(-gamma))/(1+tao_c);


i_debt=0;
i = (1-beta)/beta;

a=1;
e_a=0;
e_g=0;
g= .2; %porque incrementa como en shock sino seria cero


end;


shocks;
var e_a = sigma_a^2;
var e_g = sigma_g^2;
var e_a,e_g = 0.01*sigma_a*sigma_g;              % to calculate impulse responses set correlation = 0
end;



steady; 




stoch_simul(hp_filter = 1600, order = 1);