%Solves steady state system of equations 

% October 2019

function F=Rafaelss(w,pa,la,ka,ca,lb,kb,cb,pb,PED,Rb)

%global M_
 
% Here we load the values of the deep parameters in a loop.
 
%NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
%for i = 1:NumberOfParameters                                  % Loop...
  %paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
 % eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
%end                                                           % End of the loop.  
%check = 0;


%% Variables

%Não estava funcionando fazer F=Rafaelss(x)

%w=x(1);
%pa=x(2);
%la=x(3);
%ka=x(4);
%ca=x(5);
%lb=x(6);
%kb=x(7);
%cb=x(8);
%pb=x(9);
%PED=x(10);
%Rb = x(11);


%% Parameter definitions [Calibration from Li (2019)]

phi = 1.45;
beta = 0.998;          
alpha = 0.33;       
betaE = 0.97;         
delta = 0.025;      
a=0.5;            
b=0.5;
rho = -10;                  % quero aqui uma CES com bens complementares                     
eps = 6;                     
ma=0.8;                                          
mb=0.8;                      
N=3;                                       
Rs=0.02;                           

%% Steady state system of equations

F = [w-((la+lb)*ma*ka*((1-delta)/Rs)*((la+lb)*(beta-Rs)-beta+Rs - (phi*(a*(((ka^alpha)*((la)^1-alpha))^rho)+b*((kb^alpha)*((lb)^1-alpha))^rho)^1/rho - ca - cb - delta*(ka+kb))/((la + lb)*ma*ka* (1-delta)/Rs)))/(((la + lb)^3) - (la + lb)^2);
  pa - (((betaE*alpha)/(1-betaE*(1-delta)*(1+ (ma/betaE*Rs) - ma) ))^-alpha)*((w/1-alpha)^(1-alpha));
  la - (((1-alpha)*ka^alpha/w)^1/alpha)*(((betaE*alpha)/(1-betaE*(1-delta)*(1+ (ma/betaE*Rs) - ma) ))^-1)*((w/1-alpha)^((1-alpha)/alpha));
  ka - ((betaE * alpha *(ka^alpha)*((la)^1-alpha))/(1-betaE*(1-delta)*(1+ (ma/betaE*Rs) - ma)))*pa;
  ca - (1-betaE)*(alpha * (ka^alpha)*((la)^1-alpha))*pa + (1-betaE)*((1-delta)*ka - Rs*ma*ka*((1-delta)/Rs));
  lb - (((1-alpha)*kb^alpha/w)^1/alpha)*(((betaE*alpha)/(1-betaE*(1-delta)*(1+ (mb/betaE*Rb) - mb) ))^-1)*((w/1-alpha)^((1-alpha)/alpha));
  kb - (   (betaE * alpha *(kb^alpha)*((lb)^1-alpha))/(1-betaE*(1-delta)*(1+ (mb/betaE*Rb) - mb)))*(((betaE*alpha)/(1-betaE*(1-delta)*(1+ (mb/betaE*Rb) - mb) ))^-1)*((w/1-alpha)^((1-alpha)/alpha));
  cb - (1-betaE)*(alpha * (kb^alpha)*((lb)^1-alpha))*(((betaE*alpha)/(1-betaE*(1-delta)*(1+ (mb/betaE*Rb) - mb) ))^-alpha)*((w/1-alpha)^(1-alpha)) + (1-betaE)*((1-delta)*kb - Rb*mb*kb*((1-delta)/Rb));
  (1-eps)/eps -  (pa/(a + b * (pa * b/pb * a)^(-rho/rho-1)) + pb/(a * ( a * pa/b * pa)^(rho/1-rho) + b))^(-eps);
  PED - (1 + (1/(1-alpha))*((mB * (1 - delta))/(Rb * beta * alpha * pb * ((kb^(alpha-1)) * (lb^(1-alpha))  )    )));
  Rb - beta/(1 - (1/N) * (1/PED))];
  
%% Now put this on command window:

%Initial guess  -->>  x0 = [1.1;1;0.3;4;0.5;1;4;0.5;1;30;1]
%options = optimoptions('fsolve','Display','iter');
%Call fsolve -->> [x,fval] = fsolve(Rafaelss,x0,options)






