close all;

var c c_h im s y y_star ii i_star inflation inflation_h n w b a nx 
ex p p_star p_h c_f ss1 ptm cmt e r u v z;

varexo Taylor_shock output_shock philip_shock world_output_shock 
inflation_star_shock is_shock;

parameters theta phiy phip phiy_star phip_star mu
disc beta epsilon rhoi rhoy alpha phi psi gamma
etha delta rhoy_star rho_inflation lambda_ex etha_star rho_inflation_1
rho_is phii_star phii omega delta_star;

lambda_ex = 0.1;       
etha_star = 2.14;      
etha = 0.9;            
delta = 0.22;          
disc = 0.01;           
alpha = 0.58;          
gamma = 0.2572;        
phi = 1.6;             
psi = 0.01;            
rho_inflation_1 = 0.8; 
rhoi = 0.66;           
rhoy = 0.86;           
rhoy_star = 0.45;      
rho_inflation = 0.635; 
rho_is = 0.8;          
beta = exp(-disc);                            
phip = 1.50;           
phiy = 0.125;          
phip_star = 1.5;       
phiy_star = 0.5;       
theta = 2;             
epsilon = 6;           
omega = 0.75;                                 
phii_star = 0.50;      
phii = 0.43;           
mu = log(epsilon/(epsilon-1));
delta_star = delta*((1-alpha)^((1-alpha)/(phi+(1-alpha)*(theta-1)*(1-gamma))));  %Grado de Apertura de la Economica Externa

model;
#mut = omega/(omega-1);                                                          %mark-up firmas en Monopolio
#theta_gamma = gamma*(theta-1);
#lambda = ((1-omega)*(1-beta*omega)/omega)*((1-alpha)/(1-alpha+alpha*epsilon));  % Grado Stickness Precios
#kappa_y = lambda*((phi+alpha)/(1-alpha));                                       % Ecuacion 59 Parametros Curva Phillips NK con formacion de habitos
#kappa_a = lambda*((phi+1)/(1-alpha));                                           % Ecuacion 59 Parametros Curva Phillips NK con formacion de habitos
#ii_ss = (2/beta)-1;
#inflation_ss = 1;
#p_star_ss = 1;
#i_star_ss = ii_ss;
#y_star_ss=1;
#y_ss = (1-alpha)^((1-alpha)/(phi+(1-alpha)*(theta-1)*(1-gamma)));
#a_ss = 1;
#v_ss = 1;
#u_ss = 1;
#z_ss = 1;
#kst = (1-alpha)^(1/(phi+(1-alpha)*(theta-1)*(1-gamma)));

(c^(-theta_gamma-theta)) = beta*((c(-1)^(-theta_gamma))/(c(+1)^(theta)))*r(+1)*u; %FOC Households1
w = (n^phi)*(c(-1)^(-theta_gamma))*(c^theta);                                     %FOC Households2
y = a*(n^(1-alpha));                     %Firms Prod. Function
n = (1-alpha)*(y*kst)*(p_h/w);           %FOC Firms
c_h = (1-delta)*((p_h/p)^etha)*c;        %Domestic Consumption (Ph=P=1)
c_f = delta*((p/(e*p_star))^etha)*c;     %Foreign Consump. (Pf=P=1)
1 = z*(((1-delta)^(1/etha))*((c_h/c)^((etha-1)/etha)) + (delta^(1/etha))*((c_f/c)^((etha-1)/etha)))^(etha/(etha-1)); %Consumo Interno
ss1 = cmt + (beta*omega)*ss1(+1);
ptm = mut*(1-beta*omega)*ss1;                                                        %Precio de las (1-w) firmas que optimizan su precio
1 = (omega*((1/inflation)^(1-theta)) + (1-omega)*((ptm/p)^(1-theta)))^(1/(1-theta)); %Indice Precios Firmas (Min/Prodc)
1 = ((1-delta)*(p_h/p)^(1-etha) + delta*((e*p_star)/p)^etha)^(1/(1-etha));         %Indice Precios Total
(ex/y_star) = ((delta_star*(s^etha_star))^mu)*((ex(-1)/y_star(-1))^(1-mu));        %Eq12
nx = ex - im;
%im = c_f;
%p_f = e*p_star;
y = c_h + ex;
%q = e*p_star/p;
r = (1+ii(-1))/(1+inflation);
b = r*b(-1) + p_h(-1)*y(-1) - p(-1)*c(-1); 
(1+ii) = (1+i_star)*exp(-psi*(b/(y_ss*p_h)))*(e(+1)/e);
%s = p_f/p_h;
s = (e*p_star)/p_h;
(ii/ii_ss) = ((ii(-1)/ii_ss)^phii)*((inflation/inflation_ss)^phip)*((y/y_ss)^phiy)*v; %(23) Política Monetaria (Regla Taylor)
(i_star/i_star_ss) = ((i_star(-1)/i_star_ss)^phii_star)*((p_star/p_star_ss)^phip_star)*((y_star/y_star_ss)^phiy_star); %(3.4) Regla Taylor Política Externa
inflation = p/p(-1);
inflation_h = p_h/p(-1);
log(p_star) = (1-rho_inflation)*log(p_star_ss)+rho_inflation*log(p_star(-1)) + inflation_star_shock;    % Choque Precios Externa
log(v) = (1-rhoi)*log(v_ss) + rhoi*log(v(-1)) + Taylor_shock;                % Choque Politica Monetaria
log(a) = (1-rhoy)*log(a_ss) + rhoy*log(a(-1)) + output_shock;                % Choque Tecnologico Doméstico
log(z) = (1-rho_is)*log(z_ss) + rho_is*log(z(-1)) + is_shock;                % Choque por Demanda Doméstica
log(y_star) = (1-rhoy_star)*log(y_star_ss)+rhoy_star*log(y_star(-1)) + world_output_shock; % Choque de Producto Externo
log(u) = (1-rho_inflation_1)*log(u_ss) + rho_inflation_1*log(u(-1)) + philip_shock;        % Choque por la Curva Phillips
end;

write_latex_dynamic_model;

initval;
a=1;
p_h=1;
p=p_h;
p_star=p_h;
n=((1-alpha)*p_h)^(1/(phi+(1-alpha)*(theta-1)*(1-gamma))); 
%n=(1-alpha)^alpha;
%c=(1-alpha)^((-alpha*phi)/(gamma*(1+theta)+theta));
c=n^(1-alpha);
w=(1-alpha)*c;
%lambdat=(n^phi)/w;
y=n^(1-alpha);
%y=(1-alpha)^((1-alpha)/alpha);
%w=(1-alpha)*(y/n)*p_h; 
c_h=(1-delta)*c;
c_f=delta*c;
im=c_f;
ex=y-c_h;
%p_f=e*p_star;
nx=0;
e=1;
%p_f=e*p_star;
s=(e*p_star)/p_h;
%q=e*p_star/p;
%y_star=ex/delta_star;
y_star=1;
%mpt=(1-alpha)*a*(y/n);
%ptm=(omega/(omega-1))*cmt;
r=1/beta; 
ii=(2/beta)-1;
i_star=ii;
inflation=1;
inflation_h=1;
%inflation_star=1;
b=0;
cmt=p*((omega-1)/omega);
%cmt=(w/p_h)-(1-alpha)*a*(y/n);
ss1=cmt/(1-beta*omega);
%ss2=(y*p^omega)/(1-beta*omega);
ptm=(omega/(omega-1))*cmt;
p_star=1;
v=1;
u=1;
z=1;
end;

%shocks; %Volt. Historicas
% var is_shock = (0.029046)^2;
% var Taylor_shock = (0.04482)^2;
% var output_shock = (0.021281)^2;
% var philip_shock = (0.014155)^2;
% var world_output_shock = (0.021281)^2;
% var inflation_star_shock = (0.010215)^2;
%end;

shocks;
 var is_shock = (0.50)^2;
 var Taylor_shock = (0.50)^2;
 var output_shock = (0.50)^2;
 var philip_shock = (0.50)^2;
 var world_output_shock = (0.50)^2;
 var inflation_star_shock = (0.50)^2;
end;

check;
steady;
stoch_simul(order=1,periods=500,irf=20) w n c ex y ii inflation i_star;

%model_diagnostics(M_, options_, oo_)     
%M_.param_names(find(isnan(M_.params)),:) 
%resid(1);