clear all
close all
clc

%% Parameters

%Calibration
sigma    = 1;
varphi   = 1;
omega    = 0.65;
eta      = 1.3;
eps_H    = 11;
theta_H  = 0.7;
phi_D    = 0.3;
eta_star = 0.3;

mu = 1;
chi = 0.9;

alpha_pi = 1.5;
alpha_y = 0.5;

%Targets - normalization
PI   = 1.018;
h    = 1/3;
pH   = 1;
R    = 1.019;
PI_S = 1.00;
s_tb = 0;     %trade balance to GDP
s_g  = 0.3;   % government spending to GDP 
s_bG = 0.2;   % government debt to GDP

%Exogenous processes
z = 1;
em = 1;
g_shock = 1;

rho_z = 0.9;
rho_PI_star = 0.519; 
rho_R_W = 0.966;
rho_em = 0.9;
rho_y_star = 0.878;
rho_g = 0.9;

% Government parameters
omega_G = 0.65;
eta_G   = 1.3;

gamma_g = 0.1;
gamma_y = -0.1;
phi_bG = 0.1;
phi_y = 0.05;


%% Steady state 
beta = PI/R;
PI_star = PI/PI_S;
R_star = R/PI_S;
R_W = R_star;
PI_H = PI;
pH_tilde = 1;
Delta = pH_tilde;
mc_H = ((eps_H-1)/eps_H)*(pH_tilde^(1+eps_H));
y_H = z*h;
f1_H = ((eps_H-1)/eps_H)*pH_tilde*y_H/(1-beta*theta_H);
f2_H = (pH_tilde^-eps_H)*y_H*mc_H/(1-beta*theta_H);
pF = ((1-omega*(pH^(1-eta)))/(1-omega))^(1/(1-eta));

GDP = y_H;
pG = (omega_G*(pH^(1-eta_G)) + (1-omega_G)*(pF^(1-eta_G)))^(1/(1-eta_G));

tb = s_tb*GDP;
g  = s_g*GDP;
bG = s_bG*GDP;

c  = GDP - g - tb;

rer = pF;

c_H = omega*(pH^-eta)*c;
c_F = (1-omega)*(pF^-eta)*c;
g_H = omega_G*((pH/pG)^-eta_G)*g;
g_F = (1-omega_G)*((pF/pG)^-eta_G)*g;

c_H_star = y_H - c_H - g_H;
y_star = c_H_star*(pH/rer)^eta_star;

d_star = tb*(rer*(R_star/PI_star -1))^(-1);
d_bar = d_star;

w = pH*z*mc_H;
lambda = c^-sigma;

psi = w*lambda/(h^varphi);
PI_I = PI;

tax = g - bG*(1-R/PI);

%% Check that system of equation holds in steady state
%Eq 1:
eq1 = lambda - c^-sigma

%Eq 2:
eq2 = w*lambda - psi*(h^varphi)

%Eq 3:
eq3 = lambda - beta*R*lambda/PI

%Eq 4:
eq4 = lambda - beta*R_star*PI_S*lambda/PI

%Eq 5:
eq5 = c - (omega^(1/eta)*(c_H^(1-1/eta)) + (1-omega)^(1/eta)*(c_F^(1-1/eta)))^(eta/(eta-1)) 

%Eq 6:
eq6 = c_F - (1-omega)*(pF^-eta)*c

%Eq 7:
eq7 = c_H - omega*(pH^-eta)*c

%Eq 7:
eq8 = pH*mc_H*z - w

%Eq 9:
eq9 = f1_H - ((pH_tilde^(1-eps_H))*y_H*((eps_H-1)/eps_H) + theta_H*beta*(pH_tilde*(PI_I^mu)/pH_tilde)^(1-eps_H)*(PI_H^eps_H)*f1_H);

%Eq 10:
eq10 = f2_H - ((pH_tilde^(-eps_H))*y_H*mc_H + theta_H*beta*(pH_tilde*(PI_I^mu)/pH_tilde)^(-eps_H)*(PI_H^(1+eps_H))*f2_H);

%Eq 11:
eq11 = f1_H - f2_H

%Eq 12:
eq12 = 1 - (theta_H*(PI_H/(PI_I^mu))^(eps_H-1) + (1-theta_H)*(pH_tilde^(1-eps_H)))

%Eq 13:
eq13 = Delta - ((theta_H*(PI_H/(PI_I^mu))^eps_H)*Delta + (1-theta_H)*(pH_tilde^-eps_H))

%Eq 14:
eq14 = PI_I - (PI^chi)*(PI^(1-chi))

%Eq 15:
eq15 = R_star - (R_W*exp(phi_D*(d_star - d_bar)))

%Eq 16:
eq16 = rer - pF

%Eq 17:
eq17 = c_H_star - ((pH/rer)^(-eta_star))*y_star

%Eq 18:
eq18 = y_H - (c_H + c_H_star + g_H)

%Eq 19:
eq19 = y_H - z*h

%Eq 20:
eq20 = tb - (pH*c_H_star - rer*(c_F + g_F))

%Eq 21:
eq21 = (rer*d_star + tb) - (rer*d_star*R_star/PI_star)

%Eq 22:
eq22 = GDP - (c + g + tb)

%Eq 23:
eq23 = pH/pH - PI_H/PI

%Eq 24:
eq24 = rer/rer - PI_S*PI_star/PI

%Eq 25:
eq25 = R/R - (PI/PI)^(alpha_pi)*(GDP/GDP)^alpha_y

%Eq 26:
eq26 = pG - (omega_G*(pH^(1-eta_G)) + (1-omega_G))^(1/(1-omega_G));

%Eq 27:
eq27 = g_F - (1-omega_G)*((pF/pG)^-eta_G)*g

%Eq 28:
eq28 = g_H - omega_G*((pH/pG)^-eta_G)*g

%Eq 29
eq29 = bG - ((R/PI)*bG + g - tax) 

%Eq 30
eq30 = g - (g + gamma_g*(g-g) + gamma_y*(GDP-GDP)) 

%Eq 31
eq31 = tax - (tax + phi_bG*(bG - bG) + phi_y*(GDP-GDP))

%Eq 32: Monetary policy shock
eq32 = (em-em) -  rho_em*(em-em)            

%Eq 33: Foreign output  
eq33 = (y_star-y_star) -  rho_y_star*(y_star-y_star)  

%Eq 34: Foreign interest rate 
eq34 = (R_W-R_W) - rho_R_W*(R_W-R_W)           

%Eq 35: Foreign inflation 
eq35 = (PI_star-PI_star) - rho_PI_star*(PI_star-PI_star)

%Eq 36: Productivity shock
eq36 = (z-z) -  rho_z*(z-z) 

%Eq 37: Gov spending shock
eq37 = (g_shock-g_shock) -  rho_g*(g_shock-g_shock) 


%% Check steady state 
tol = 1e-6; 

% Check if any of the variables are above the threshold
if any([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13,...
        eq14, eq15, eq16, eq17, eq18, eq19, eq20, eq21, eq22, eq23, eq24, eq25,...
        eq26, eq27, eq28, eq29, eq30, eq31, eq32, eq33, eq34, eq35, eq36, eq37] > tol)
    disp('Could not find the steady state');
else
    % Check if all variables are below the threshold
    if all([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13,...
        eq14, eq15, eq16, eq17, eq18, eq19, eq20, eq21, eq22, eq23, eq24, eq25,...
        eq26, eq27, eq28, eq29, eq30, eq31, eq32, eq33, eq34, eq35, eq36, eq37] <= tol)
        disp('Done - we found the steady state');
%     else
%         disp('Some variables are above the threshold, but not all');
      end
end


 
%% Export 
% Rename variables (add "_ss" at the end)
variables = {'c', 'c_H', 'c_F', 'h', 'y_H','lambda',...
             'w', 'pH', 'pF', 'mc_H',...
             'pH_tilde', 'f1_H', 'f2_H', 'Delta',...
             'PI', 'PI_H', 'PI_S', 'PI_I',...
             'rer', 'c_H_star', 'tb','d_star', 'R_star','GDP',...
             'R',...
             'R_W', 'y_star', 'PI_star','z','em',...
             'g','g_H','g_F','pG','bG','tax','g_shock'};

% Iterate through each variable
for j = 1:numel(variables)
    % Construct the new variable name
    new_variable_name = [variables{j}, '_ss'];
    
    % Get the value of the current variable
    current_value = eval(variables{j});
    
    % Assign the value to the new variable name
    assignin('base', new_variable_name, current_value);
    
    % Clear the old variable
    eval(['clear ', variables{j}]);
end




save parameters beta psi sigma varphi omega eta eps_H theta_H phi_D eta_star alpha_pi alpha_y mu chi...
    rho_R_W rho_y_star rho_PI_star rho_z rho_em ...
    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...
    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 d_bar R_W_ss y_star_ss PI_star_ss z_ss em_ss GDP_ss...
    s_g s_bG omega_G eta_G gamma_g gamma_y phi_bG phi_y g_ss g_H_ss g_F_ss pG_ss bG_ss tax_ss g_shock_ss rho_g;
       


display('Doneeeee')   

disp(' ');


%% Run dynare file
dynare model_gov_fixed


%% Plots
%doing_plots