close all
clear all
clc 

%% Parameters
r     = 1/100;
sigma = 2;
psi   = 0.07;
omega = 3;
chi   = 0.44;
S_Co  = 7.98/100;
S_D   = 17.89/100;
mu    = 0.5;
alpha_T = 0.79;
alpha_N = 0.25;

R  = r;

% Parameters for exogenous process
rho_pCo = 0.6868;
sigma_pCo = 0.0933;


%% Steady state 
AN = 1;
AT = 1;
pCo = 1;
p = 1;
pN = 1;

%S.1
beta = 1/(1+r);

%S.2
N1N_N2N = (((1-alpha_T)/(1-alpha_N))*(alpha_N/alpha_T)^(alpha_T/(alpha_T-1)))^((alpha_T-1)/(alpha_T-alpha_N));

%S.3
N1T_N2T = ((alpha_N/alpha_T)*(N1N_N2N)^(alpha_N-1))^(1/(alpha_T-1));

%S.5
N1 = (alpha_T*(N1T_N2T)^(alpha_T-1))^(1/(omega-1));

%S.6
N2 = ((1-alpha_T)*(N1T_N2T)^(alpha_T))^(1/(omega-1));

%S.7
N2N = (N1 - N1T_N2T*N2)/(N1N_N2N - N1T_N2T);

%S.8
N1N = N1N_N2N*N2N;

%S.9
N1T = N1 - N1N;

%S.10
N2T = N2 - N2N;

%S.11
YT = (N1T^alpha_T)*(N2T^(1-alpha_T));

%S.12
YN = (N1N^alpha_N)*(N2N^(1-alpha_N));

%S.13
Y = (YT + YN)/(1-S_Co);

%S.14
Co = S_Co*Y;

%S.17
D = S_D*Y;

%S.18
w1 = alpha_T*YT/N1T;

%S.19
w2 = (1-alpha_T)*YT/N2T;

%S.34
d2 = D/2;
d1 = D/2;

%S.20
c1 = w1*N1 - r*d1 + (1-mu)*pCo*Co;

%S.21
c2 = w2*N2 - r*d2 + mu*pCo*Co;

%S.22
C = c1 + c2;

%S.23
phi = YN/C;

%S.24
c1_N = phi*c1;

%S.25
c1_T = (1-phi)*c1;

%S.26
c2_N = phi*c2;

%S.27
c2_T = (1-phi)*c2;

%S.30
C_N = c1_N + c2_N;

%S.31
C_T = c1_T + c2_T;

%S.32 
wT = (N1T/(N1T + N2T))*w1 + (N2T/(N1T + N2T))*w2;

%S.33 
wN = (N1N/(N1N + N2N))*w1 + (N2N/(N1N + N2N))*w2;


% Total debt 
D = d1 + d2;

%% check that system of equation holds in steady state
% %Eq 1:
% eq1 = (c1 - (1/omega)*(N1^omega))^(-sigma) - beta*((c1 - (1/omega)*(N1^omega))^(-sigma)*(p/p)*(1+R))
% 
% %Eq 2:
% eq2 = (c2 - (1/omega)*(N2^omega))^(-sigma) - beta*((c2 - (1/omega)*(N2^omega))^(-sigma)*(p/p)*(1+R))
% 
% %Eq 3:
% eq3 = c1_T + pN*c1_N + (1+R)*d1 - (alpha_T*YT*N1/N1T + d1 + (1-mu)*pCo*Co)
% 
% %Eq 4:
% eq4 = c2_T + pN*c2_N + (1+R)*d2 - ((1-alpha_T)*YT*N2/N2T + d2 + mu*pCo*Co)
% 
% %Eq 5:
% eq5 = c1_N - phi*(pN/p)^(-chi)*c1
% 
% %Eq 6:
% eq6 = c2_N - phi*(pN/p)^(-chi)*c2
% 
% %Eq 7:
% eq7 = c1_T - (1-phi)*(1/p)^(-chi)*c1
% 
% %Eq 8:
% eq8 = c2_T - (1-phi)*(1/p)^(-chi)*c2
% 
% %Eq 9:
% eq9 = C - (c1 + c2)
% 
% %Eq 10:
% eq10 = C_N - (c1_N + c2_N)
% 
% %Eq 11:
% eq11 = C_T - (c1_T + c2_T)
% 
% %Eq  12: 
% eq12 = p - (phi*(pN^(1-chi)) + (1-phi))^(1/(1-chi))
% 
% %Eq 13:
% eq13 = YT - AT*(N1T^alpha_T)*(N2T^(1-alpha_T))
% 
% %Eq 14:
% eq14 = YN - AN*(N1N^alpha_N)*(N2N^(1-alpha_N))
% 
% %Eq 15:
% eq15 = N1^(omega-1) - alpha_T*YT/(N1T*p)
% 
% %Eq 16:
% eq16= N1^(omega-1) - pN*alpha_N*YN/(N1N*p) 
% 
% %Eq 17:
% eq17 = N2^(omega-1) - (1-alpha_T)*YT/(N2T*p)
% 
% %Eq 18:
% eq18 = N2^(omega-1) - pN*(1-alpha_N)*YN/(N2N*p) 
% 
% %Eq 19:
% eq19 = R - (r+psi*(exp(d1+d2-d1-d2)-1))
% 
% %Eq 20:
% eq20 = YN - C_N 
% 
% %Eq 21:
% eq21 = N1 - (N1T + N1N)
% 
% %Eq 22:
% eq22 = N2 - (N2T + N2N)
% 
% %Eq 23:
% eq23 = log(pCo) - ((1-rho_pCo)*log(pCo) + rho_pCo*log(pCo))


%% Export 
% Rename variables (add "_ss" at the end)
variables = {'YT', 'YN', 'C', 'C_N', 'C_T', 'c1', 'c2', 'c1_N', 'c1_T', 'c2_N', 'c2_T',...
             'N1N', 'N1T', 'N2N', 'N2T', 'N1', 'N2', 'd1', 'd2',...
             'pN', 'p', 'pCo', 'w1', 'w2', 'D','R'};

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


save parameters r beta sigma psi omega chi phi alpha_T alpha_N Co...
                S_Co S_D rho_pCo sigma_pCo AT AN mu...
                YT_ss YN_ss C_ss C_N_ss C_T_ss c1_ss c2_ss c1_N_ss c1_T_ss c2_N_ss c2_T_ss...
                N1N_ss N1T_ss N2N_ss N2T_ss N1_ss N2_ss d1_ss d2_ss...
                pN_ss p_ss pCo_ss D_ss R_ss ;
       
       

%% Dynare code 
dynare code 
