% PC3 - Modelo RBC estándar de SOE con costos de ajustes de capital y portafolio

%Declaración de variables endógenas
var  c h y i k a lambda ${\lambda}$ util;  
%Declaration of exogenous variables
varexo e e_r;   

%Calibración de parámetros                                     
parameters  gamma ${\gamma}$
            omega ${\omega}$
            rho ${\gamma}$
            sigma_tfp ${\sigma_{a}}$
            sigma_r ${\sigma_{r}}$
            delta ${\delta}$
            alpha ${\alpha}$
            psi_3 ${\psi_3}$
            phi ${\phi}$
            r_bar ${\bar r}$
            d_bar ${\bar d}$;

%Calibración
gamma  = 2; %aversión al riesgo 
omega  = 1.455; %elasticidad Frisch
alpha  = 0.32; %labor share
r_bar    = 0.04; %world interest rate		
delta  = 0.1; %depreciation rate
rho    = 0.42; %autocorrelation TFP 
sigma_tfp = 0.0129; %standard deviation TFP
sigma_r = 0.01; %standard deviation Interest Rate
phi = 0.28; %costo de ajuste de capital
d_bar  = 0.7442;
psi_3  = 0.00074;   

%Declaración de variables endógenas
var d tb_y, ca_y, r;

%Parametrización del factor de descuento
parameters beta ${\beta}$;

%Declaración del modelo
model;
[name='Eq. (1), Evolución de la deuda']
    d = (1+exp(r(-1)))*d(-1)- exp(y)+exp(c)+exp(i)+phi*(exp(k(+1))-exp(k))^2+psi_3*(d-d_bar)^2;
[name='Eq. (4), Función de produciión']
    exp(y) = exp(a)*(exp(k(-1))^alpha)*(exp(h)^(1-alpha));
[name='Eq. (2), Ley del movimiento del capital']
    exp(k) = exp(i)+(1-delta)*exp(k(-1)); 
[name='Eq. (9),  Ecuación de euler']
    exp(lambda)*(1-psi_3*(d-d_bar))= beta*(1+exp(r))*exp(lambda(+1)); 
[name='Eq. (6), Utilidad marginal']
    (exp(c)-((exp(h)^omega)/omega))^(-gamma)   = exp(lambda);  
[name='Eq. (7), CPO Trabajo']
    ((exp(c)-((exp(h)^omega)/omega))^(-gamma))*(exp(h)^(omega-1))  = exp(lambda)*(1-alpha)*exp(y)/exp(h);
[name='Eq. (8), CPO Inversión']
    exp(lambda) = beta*exp(lambda(+1))*(exp(a(+1))*alpha*exp(y(+1))/exp(k)+(1-delta)+phi*(exp(k(+1))-exp(k)));
[name='Eq. (*), Ley de movimiento de la PTF']
    a = rho*a(-1)+sigma_tfp*e; 
[name='Eq. (3), Tasa de interés']
    exp(r) = r_bar + sigma_r*e_r;

[name='Definición de balanza comercial/producto']
    tb_y = 1-(exp(c)+exp(i))/exp(y);
    ca_y = (1/exp(y))*(d(-1)-d); 
    [name='Definición de función de utilidad']    
    util=(((exp(c)-omega^(-1)*exp(h)^omega)^(1-gamma))-1)/(1-gamma);
end;

%Declaración del estado estacionario
steady_state_model;
    beta  = 1/(1+r_bar);
    r     = log((1-beta)/beta);
    d     = d_bar;
    h     = log(((1-alpha)*(alpha/(r_bar+delta))^(alpha/(1-alpha)))^(1/(omega-1)));
    k     = log(exp(h)/(((r_bar+delta)/alpha)^(1/(1-alpha))));
    y     = log((exp(k)^alpha)*(exp(h)^(1-alpha)));
    i     = log(delta*exp(k));
    c     = log(exp(y)-exp(i)-r_bar*d); 
    tb_y    = 1-((exp(c)+exp(i))/exp(y));
    util=(((exp(c)-omega^(-1)*exp(h)^omega)^(1-gamma))-1)/(1-gamma);
    lambda= log((exp(c)-((exp(h)^omega)/omega))^(-gamma));
    a     = 0;
    ca_y    = 0;
end;

%Print the residual of each equation
resid;

%Print los valores del estado estacionario
steady; 

%Declaración de los shocks
// Generate IRFs 
shocks;
     var e; stderr 1/sigma_tfp; %Normalize to unit shock
end;

stoch_simul(order=1, hp_filter=100, periods=0);

//Report results
y_pos=strmatch('y',M_.endo_names,'exact');
c_pos=strmatch('c',M_.endo_names,'exact');
i_pos=strmatch('i',M_.endo_names,'exact');
h_pos=strmatch('h',M_.endo_names,'exact');
tb_y_pos=strmatch('tb_y',M_.endo_names,'exact');
ca_y_pos=strmatch('ca_y',M_.endo_names,'exact');


fprintf('\nstd(y):              \t %2.1f \n',sqrt(oo_.var(y_pos,y_pos))*100)
fprintf('std(c):                \t %2.1f \n',sqrt(oo_.var(c_pos,c_pos))*100)
fprintf('std(i):                \t %2.1f \n',sqrt(oo_.var(i_pos,i_pos))*100)
fprintf('std(h):                \t %2.1f \n',sqrt(oo_.var(h_pos,h_pos))*100)
fprintf('std(tb/y):             \t %2.1f \n',sqrt(oo_.var(tb_y_pos,tb_y_pos))*100)
fprintf('std(ca/y):             \t %2.1f \n',sqrt(oo_.var(ca_y_pos,ca_y_pos))*100)
fprintf('corr(y_t,y_t-1):       \t %3.2f \n',oo_.autocorr{1}(y_pos,y_pos))
fprintf('corr(c_t,c_t-1):       \t %3.2f \n',oo_.autocorr{1}(c_pos,c_pos))
fprintf('corr(i_t,i_t-1):       \t %4.3f \n',oo_.autocorr{1}(i_pos,i_pos))
fprintf('corr(h_t,h_t-1):       \t %3.2f \n',oo_.autocorr{1}(h_pos,h_pos))
fprintf('corr(tb/y_t,tb/y_t-1): \t %3.2f \n',oo_.autocorr{1}(tb_y_pos,tb_y_pos))
fprintf('corr(ca/y_t,ca/y_t-1): \t %3.2f \n',oo_.autocorr{1}(ca_y_pos,ca_y_pos))
fprintf('corr(c_t,y_t):         \t %3.2f \n',oo_.gamma_y{1}(c_pos,y_pos)/sqrt(oo_.var(c_pos,c_pos)*oo_.var(y_pos,y_pos)))
fprintf('corr(i_t,y_t):         \t %3.2f \n',oo_.gamma_y{1}(i_pos,y_pos)/sqrt(oo_.var(i_pos,i_pos)*oo_.var(y_pos,y_pos)))
fprintf('corr(h_t,y_t):         \t %2.1f \n',oo_.gamma_y{1}(h_pos,y_pos)/sqrt(oo_.var(h_pos,h_pos)*oo_.var(y_pos,y_pos)))
fprintf('corr(tb/y_t,y_t):      \t %4.3f \n',oo_.gamma_y{1}(tb_y_pos,y_pos)/sqrt(oo_.var(tb_y_pos,tb_y_pos)*oo_.var(y_pos,y_pos)))
fprintf('corr(ca/y_t,y_t):      \t %4.3f \n',oo_.gamma_y{1}(ca_y_pos,y_pos)/sqrt(oo_.var(ca_y_pos,ca_y_pos)*oo_.var(y_pos,y_pos)))