//RBC model

/// VARIABLES 

var c k a h d y i t_b lambda t_b_y g_y g_c g_i g r ;

varexo epsilon_a epsilon_g  epsilon_nu epsilon_s epsilon_mu; //epsilon_y_ME 

/// PARAMETERS

parameters alpha phi delta gamma theta omega beta psi d_bar g_bar rho_a rho_g rho_nu rho_s rho_mu; 

/// CALIBRATION 

alpha = 0.32;
phi = 0.588;
delta = 0.076;
gamma = 2;
theta = 2.24;
omega = 1.5;
beta = 0.99; // Households' discount factor 
psi = 0.001;
d_bar = 0.007;
g_bar = 1.0050;

// SHOCK PARAMETERS

rho_a = 0.75; // Persistence of stationary productivity shock
rho_g = 0.75;  // Persistence of nonstationary productivity shock 
//rho_nu = 0.75; // Persistence of preference shock 
//rho_s = 0.75; // Persistance of domestic spending shock
//rho_mu = 0.75; // Persistance of country premium shock 

/// MODEL      
                      
model(linear);

#r_star=1/beta+gamma*g_bar;

// Output 
y=a+alpha*k+(1-alpha)*(g+h);
 
// Resource constraint 
y=t_b+c+i+phi*(k(+1)-g_bar); 

// Resource constraint in the case of RBCFF
//// y=d-d(+1)+c+i+r++s+phi*(k(+1)-g_bar); 

// Capital stock 
//i=k(+1)*g - (1-delta)*k;
k*g=i+(1-delta)*k(-1);

// Marginal utility of concumption  
gamma*(g(-1)-c)-theta*h=lambda; 

// First order condition for labour
(1+gamma)*g(-1)-gamma*c+((theta*omega)-alpha)*h=(1-alpha)*a+alpha*k+(1-alpha)*g+lambda;

// Euler equation 
lambda(-1)=beta*r-gamma*g+lambda;   
 
// First order condition for capital 
lambda+(beta*gamma)*g-phi*g_bar=lambda(+1)+(1-delta)+alpha*a(+1)+(1-alpha)*h(+1)-(1-alpha)*k(+1)+phi*k(+2); 

// World interest rate  
r+r_star+psi*(d(+1)-g-d_bar); 

// DEFINITIONS 

// Trade balance
t_b = d - d(-1) + g - r; 

// Trade balance to output ratio
t_b_y = t_b - y; 

// Output growth
g_y = y - y(-1) +g(-1); 

// Consumption growth
g_c = c - c(-1) +g(-1); 

// Investment growth
g_i = i - i(-1) +g(-1); 

// EXOGENOUS AR(1) PROCESS  
 
// Stationary productivity shock
a = rho_a*a(-1) + epsilon_a; //epsilon_y_ME; 

// Nonstationary productivity shock
g-g_bar=rho_g*(g(-1)-g_bar)+ epsilon_g; 

// Preference shock in the case of RBCFF    
////nu=rho_nu*nu(-1) + epsilon_nu; 

// Domestic spending shock in the case of RBCFF  
////s-sbar=rho_s*(s(-1)-sbar)+ epsilon_s; 

// Exogenous stochastic country premium shock in the case of RBCFF  
////mu=rho_mu*mu(-1) + epsilon_mu; 

end;   

shocks;

var epsilon_a; stderr 0.75;
//var epsilon_y_ME; stderr 0.75; 
var epsilon_g; stderr 0.75;
var epsilon_nu; stderr 0.75;
var epsilon_s; stderr 0.75;
var epsilon_mu; stderr 0.75;

end;

//resid(1);
steady;
check; 
  
//stoch_simul(order=1,irf=0)g_y g_c g_i t_b_y; //

//ESTIMATION

varobs g_y g_c g_i t_b_y;

estimated_params;

// Average growth rate
g_bar, , , , uniform_pdf,(1+1.03)/2,sqrt(12)^(-1)*(1.03-1), 1, 1.03;
stderr epsilon_g, , , , uniform_pdf,(0+0.2)/2,sqrt(12)^(-1)*(0.2-1), 0, 0.2;
// Inverse of elasticity of labour
gamma, , , , uniform_pdf,(1.478+2.447)/2,sqrt(12)^(-1)*(2.447-1.478), 1.478, 2.447; //***

// Capital adjustment cost
phi, , , , uniform_pdf,(0.588+0.174)/2,sqrt(12)^(-1)*(0.174-0.588), 0.588, 0.174; //**

// Country risk premium
psi, , , , uniform_pdf,(0+5)/2,sqrt(12)^(-1)*(5-0), 0, 5;

//Persistance parameters
rho_a, , , , uniform_pdf,(-0.99+0.99)/2,sqrt(12)^(-1)*(-0.99-0.99), -0.99, 0.99;
rho_g, , , , uniform_pdf,(-0.99+0.99)/2,sqrt(12)^(-1)*(-0.99-0.99), -0.99, 0.99;
stderr epsilon_a, , , , uniform_pdf,(0+0.2)/2,sqrt(12)^(-1)*(0.2-1), 0, 0.2;
//rho_nu, , , , uniform_pdf,(0.75+0.1)/2,sqrt(12)^(-1)*(0.1-0.75), 0.75, 0.1;
//rho_s, , , , uniform_pdf,(0.75+0.1)/2,sqrt(12)^(-1)*(0.1-0.75), 0.75, 0.1;
//rho_mu, , , , uniform_pdf,(0.75+0.1)/2,sqrt(12)^(-1)*(0.1-0.75), 0.75, 0.1;

///Structural shocks
//Measurement Errors Using Dynare Capabilities 
stderr g_y, , , , uniform_pdf,(0.0001+0.013)/2,sqrt(12)^(-1)*(0.013-0.0001),0.0001, 0.013; 
stderr g_c, , , , uniform_pdf,(0.0001+0.019)/2,sqrt(12)^(-1)*(0.019-0.0001), 0.0001, 0.019; 
stderr g_i, , , , uniform_pdf,(0.0001+0.051)/2,sqrt(12)^(-1)*(0.051-0.0001), 0.0001, 0.051; 
stderr t_b_y, , , , uniform_pdf,(0.0001+0.013)/2,sqrt(12)^(-1)*(0.013-0.0001), 0.0001, 0.013; 

//stderr epsilon_y_ME, , , , uniform_pdf,(1+inf)/2,sqrt(12)^(-1)*(inf-1), 1, inf; 

//stderr epsilon_nu, , , , uniform_pdf,(1+inf)/2,sqrt(12)^(-1)*(inf-1), 1, inf;
//stderr epsilon_s, , , , uniform_pdf,(1+inf)/2,sqrt(12)^(-1)*(inf-1), 1, inf;
//stderr epsilon_mu, , , , uniform_pdf,(1+inf)/2,sqrt(12)^(-1)*(inf-1), 1, inf;

end; 



estimated_params_init(use_calibration); 
    stderr g_y, 0.000114861607534;
    stderr g_c, 0.000130460798135;
    stderr g_i, 0.001436553959841;
    stderr t_b_y, 0.000109652068259;

end;

estimation (datafile=Book1, mode_check, mode_compute=8, moments_varendo, mh_nblocks=1, mh_replic=2000000) g_y g_c g_i t_b_y;