% Calcula o estado estacionario do modelo
% Parametros
clear;
global beta a1 phi delta tau T eta omega alpha fi nu psi1 A1 A2
% Parametros principais
beta    = 0.985;
a1      = 0.5;
phi     = 6.7;
delta   = 0.025;
% Parametros para funcao
tau     = 0.02;
T       = 0;
eta     = 2; 
omega   = 1.5;
alpha   = 0.39;
fi      = 6.7;
nu      = 0.4;
psi1    = 0.040228;
A1      = 1;
A2      = 1;


xss0 = [1.6134  % Y_1
        1.6134  % Y_2
        9.4115  % K_1
        9.4115  % K_2
        0.6242  % L_1
        0.6242  % L_2
        0.8608  % N
        1       % P_1
        1       % P_2
        1.3791  % W
        1.00097 % R_fd
        ];
options = optimset('Display','off');

[xss,Fval,exitflag] = fsolve(@SS_k,xss0,options);
exitflag;
Y_1 = xss(1);
Y_2 = xss(2);
K_1 = xss(3);
K_2 = xss(4);
L_1 = xss(5);
L_2 = xss(6);
N = xss(7);
P_1 = xss(8);
P_2 = xss(9);
W = xss(10);
R_fd = xss(11);

% Equacoes estaticas - SS
    % Eq. 3.6
    R_n = 1/beta;
    % Eq. 3.7
    Y = ((1-a1)^(1/phi)*Y_1^((phi-1)/phi) + a1^(1/phi)*Y_2^((phi-1)/phi))^(phi/(phi-1));
    % Eq. 3.9 
    P = ((1-a1)*P_1^(1-phi) + a1*P_2^(1-phi))^(1/(1-phi));
    % Eq. 3.31
    Q = - P;
    % Eq. 3.30
    R = Q*(beta*(1-delta) -1)*(1/beta);
    % Eq. 3.32
    R_f = R_n;
    % Eq. 3.37
    L = L_1 + L_2;
    % Eq. 3.38
    K = K_1 + K_2;
    % Eq. 3.26
    I = delta*K;
    % Eq. 3.36
    C = Y - I;
    
% Guardando o resultado final no mesmo vetor
xss(12) = R_n;
xss(13) = Y;
xss(14) = P;
xss(15) = Q;
xss(16) = R;
xss(17) = R_f;
xss(18) = L;
xss(19) = K;
xss(20) = I;
xss(21) = C;

% Guarda valores em arquivo
csvwrite('steady_state.csv',xss)