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

close all;

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

%% Endogenous variables - 25
var c c_H c_F h y_H lambda GDP; // 7
var w pH pF mc_H pH_tilde f1_H f2_H Delta; //8
var PI PI_H PI_S PI_I; // 4
var rer c_H_star tb d_star R_star R; // 6
var g g_H g_F pG bG tax;

%% Exogenous processes - 6
var R_W y_star PI_star z em g_shock; // 6

%% Disturbances variables - 6
varexo u_R_W u_y_star u_PI_star u_z u_em u_g;

%% Parameters 
parameters beta psi sigma varphi omega eta eps_H theta_H phi_D eta_star alpha_pi alpha_y mu chi;
parameters rho_R_W rho_y_star rho_PI_star rho_z rho_em;
parameters c_ss c_H_ss c_F_ss h_ss y_H_ss lambda_ss w_ss pH_ss pF_ss mc_H_ss pH_tilde_ss f1_H_ss f2_H_ss Delta_ss;
parameters PI_ss PI_H_ss PI_S_ss PI_I_ss rer_ss c_H_star_ss tb_ss d_star_ss  R_star_ss R_ss R_W_ss y_star_ss PI_star_ss z_ss em_ss GDP_ss d_bar;
parameters sigma_R_W sigma_y_star sigma_PI_star sigma_z sigma_em;
parameters s_g s_bG omega_G eta_G gamma_g gamma_y phi_bG phi_y rho_g;
parameters g_ss g_H_ss g_F_ss pG_ss bG_ss tax_ss g_shock_ss;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------
% Load the parameters file

params = load('parameters.mat');

% Get all parameter names (field names of the struct)
param_names = fieldnames(params);

% Loop through each parameter
for i = 1:numel(param_names)
    current_variable = param_names{i};
    current_value = params.(current_variable);

    % Set the parameter in Dynare
    set_param_value(current_variable, current_value);
end


sigma_R_W     = 0.01;
sigma_y_star  = 0.01;
sigma_PI_star = 0.01;
sigma_z       = 0.01;
sigma_em      = 0.01;
sigma_g       = 0.01;

%----------------------------------------------------------------
% 3. Model - non linear version 
%----------------------------------------------------------------
model; 
% Households' part %%%%%%%%%%%%%%%%%% 

% (1) FOC consumption
exp(lambda) = exp(c)^(-sigma);  

% (2) FOC labor                                           
psi*exp(h)^varphi = exp(lambda)*exp(w);                                    

%(4) Euler equation foreign bonds (UIP)
exp(lambda) = beta*exp(R_star)*exp(lambda(+1))*exp(PI_S(+1))/exp(PI(+1));      
 
%(3) Euler equation domestic bonds
exp(lambda) = beta*exp(R)*exp(lambda(+1))/exp(PI(+1)); 

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Consumption aggregator %%%%%%%%%%%%%
% (5) Aggregate consumption
exp(c) = (omega^(1/eta)*(exp(c_H)^(1-1/eta)) + (1-omega)^(1/eta)*(exp(c_F)^(1-1/eta)))^(eta/(eta-1)); 

% (6) NoN-tradable consumption
exp(c_H) = omega*exp(pH)^(-eta)*exp(c);                                

% (7) Tradable consumption 
exp(c_F) = (1-omega)*(exp(pF)^-eta)*exp(c);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sticky prices in the home-good sector %%%%%%%%%%%%%%%%%%%%%%%

% (8) Marginal cost 
exp(pH)*exp(mc_H)*exp(z) = exp(w);

% (9) Equation 1 for optimal price 
exp(f1_H) = (exp(pH_tilde)^(1-eps_H))*exp(y_H)*((eps_H-1)/eps_H) 
+ theta_H*beta*(exp(lambda(+1))/(exp(lambda)*exp(PI(+1))))*(exp(pH_tilde)*(exp(PI_I(+1))^mu)/exp(pH_tilde(+1)))^(1-eps_H)*(exp(PI_H(+1))^eps_H)*exp(f1_H(+1));
                     
% (10) Equation 2 for optimal price 
exp(f2_H) = (exp(pH_tilde)^(-eps_H))*exp(y_H)*exp(mc_H) 
+ theta_H*beta*(exp(lambda(+1))/(exp(lambda)*exp(PI(+1))))*(exp(pH_tilde)*(exp(PI_I(+1))^mu)/exp(pH_tilde(+1)))^(-eps_H)*(exp(PI_H(+1))^(1+eps_H))*exp(f2_H(+1));

% (11) Equations 3 for optimal price 
exp(f1_H) = exp(f2_H);

% (12) Equations for price evolution
1 = theta_H*(exp(PI_H)/(exp(PI_I)^mu))^(eps_H-1) + (1-theta_H)*(exp(pH_tilde)^(1-eps_H));

% (13) Price dispersion 
exp(Delta) = (1-theta_H)*exp(pH_tilde)^(-eps_H) + theta_H*(exp(PI_H)/exp(PI_I)^mu)^(eps_H)*exp(Delta(-1)); //E17

% (14) Indexation rule  
exp(PI_I) = (exp(PI(-1))^chi)*(PI_ss^(1-chi));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% National identities and definitions 
% (15) Foreign interest rate 
%exp(R_star) = exp(R_W) + phi_D*(exp(d_star - d_bar)-1);
exp(R_star) = exp(R_W) + phi_D*(exp(d_star - d_bar) - 1);


% (16) Real exchange rate 
exp(rer) = exp(pF);

% (17) Foreign demand for domestic good 
exp(c_H_star) = ((exp(pH)/exp(rer))^(-eta_star))*exp(y_star);

% (18) Home good market clearing 
exp(y_H) = exp(c_H) + exp(c_H_star) + exp(g_H);

% (19) Aggregate output
exp(y_H) = exp(z)*exp(h)/exp(Delta);

% (20) Trade balance 
tb = exp(pH)*exp(c_H_star) - exp(rer)*(exp(c_F) + exp(g_F));

% (21) Aggregate market clearing 
exp(rer)*d_star + tb = (exp(rer)/exp(PI_star))*exp(R_star(-1))*d_star(-1);                                                           // E19

% (22) GDP
exp(GDP) = exp(c) + tb + exp(g);

% Relative prices %%%%%%%%%%%%%%%%%%%%%%

% (23) Home goods inflation 
exp(pH)/exp(pH(-1)) = exp(PI_H)/exp(PI(-1)); 

% (24) Foreign good inflation
exp(rer)/exp(rer(-1)) = exp(PI_S)*exp(PI_star)/exp(PI);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% (25) Monetary rule
%(exp(R)/R_ss) = (exp(PI)/PI_ss)^(alpha_pi)*(exp(GDP)/GDP_ss)^alpha_y*exp(em);
exp(PI_S)/PI_S_ss = 1; 

%Government block 

% (26)Government spending price index 
exp(pG) = (omega_G*(exp(pH)^(1-eta_G)) + (1-omega_G)*(exp(pF)^(1-eta_G)))^(1/(1-eta_G));

% (27) NoN-tradable gov consumption
exp(g_H) = omega*(exp(pH)/exp(pG))^(-eta)*exp(g);                                

% (28) Tradable gov consumption 
exp(g_F) = (1-omega)*(exp(pF)/exp(pG))^(-eta)*exp(g);

% (29) Government budget constrtaint
bG  = (exp(R(-1))/exp(PI))*bG(-1) + exp(g) - exp(tax); 

% (30) Government spending rule 
g - log(g_ss) = gamma_g*(g(-1)-log(g_ss)) + gamma_y*(GDP-log(GDP_ss)) + g_shock;

% (31) Government tax rule 
%tax - log(tax_ss) = phi_bG*(bG(-1) - bG_ss) + phi_y*(GDP-log(GDP_ss));
tax - log(tax_ss) = phi_bG*(bG(-1)/exp(GDP(-1)) - s_bG) + phi_y*(GDP-log(GDP_ss));


% Exogenous processes 
% (32) Monetary policy shock
em-log(em_ss) = rho_em*(em(-1)-log(em_ss))-u_em;            

% (33) Foreign output  
y_star-log(y_star_ss) = rho_y_star*(y_star(-1)-log(y_star_ss))-u_y_star;             

% (34) Foreign interest rate 
R_W-log(R_W_ss) = rho_R_W*(R_W(-1)-log(R_W_ss))+u_R_W;             

% (35) Foreign inflation 
PI_star-log(PI_star_ss) = rho_PI_star*(PI_star(-1)-log(PI_star_ss))+u_PI_star; 

% (36) Productivity shock
z-log(z_ss) = rho_z*(z(-1)-log(z_ss))+u_z; 

% (37) Government spending shock
g_shock-log(g_shock_ss) = rho_g*(g_shock(-1)-log(g_shock_ss))+u_g; 


end; 

%----------------------------------------------------------------
% 4. Steady State
%----------------------------------------------------------------

steady_state_model;
c = log(c_ss);
c_H = log(c_H_ss);
c_F = log(c_F_ss);
h = log(h_ss); 
y_H = log(y_H_ss);
lambda = log(lambda_ss);
GDP = log(GDP_ss); 
w = log(w_ss);
pH = log(pH_ss); 
pF = log(pF_ss); 
mc_H = log(mc_H_ss); 
pH_tilde = log(pH_tilde_ss);
f1_H = log(f1_H_ss); 
f2_H = log(f2_H_ss); 
Delta = log(Delta_ss); 
PI = log(PI_ss); 
PI_H = log(PI_H_ss); 
PI_S = log(PI_S_ss); 
PI_I = log(PI_I_ss); 
rer = log(rer_ss); 
c_H_star = log(c_H_star_ss); 
tb = tb_ss;
d_star = d_star_ss;
R_star = log(R_star_ss);
R = log(R_ss);
R_W = log(R_W_ss); 
y_star = log(y_star_ss);
PI_star = log(PI_star_ss);
z = log(z_ss);
em = log(em_ss);


g = log(g_ss);
g_H = log(g_H_ss);
g_F = log(g_F_ss);
pG  = log(pG_ss);
bG  = bG_ss;
tax = log(tax_ss);
g_shock = log(g_shock_ss);

end;

model_diagnostics;
%resid(1)
resid;
steady;
check;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------
steady;
check;

%----------------------------------------------------------------
% 7. Estimation
%----------------------------------------------------------------
shocks;
var u_R_W     = sigma_R_W^2;
var u_y_star  = sigma_y_star^2;
var u_PI_star = sigma_PI_star^2;
var u_z       = sigma_z^2;
var u_em      = sigma_em^2;
var u_g       = sigma_g^2;
end;


%--------------------------------------------------------------------------
% 8. Simulation
%--------------------------------------------------------------------------


// RUNNING
stoch_simul(periods=0, irf=100, order=1, ar=50, nograph);
save model_GM_gov_fixed.mat M_ oo_ options_;