
function [ ys,params, check ] = MODELFNG_steadystate(ys,exo,M_,options_)
%main_steadystate()
% this function solves for the steady state of the model in "main.mod"

% declare the global variable M:


% make the fsolve functions available to the steady state file
%addpath('\\srv9dnbfil002\userhome\Kons-FNG\Desktop\SS\solvers');

% OTD_steadystate() Summary of this function goes here
%   STEADY STATE VERSION 1: Solves for the steady state (and 3 parameter
%   calibrations), given all other parameters and 3 steady state targets. 

%% Set the Global Variable
global M_

% file must find check = 0 to verify steady state
check = 1;

%% load parameters:
%for jj = 1:M_.param_nbr 
    %paramname = deblank(M_.param_names(jj,:));
  
   %eval([paramname ' = M_.params(' int2str(jj) ');']);

%end

%%
for ii = 1:M_.param_nbr
  eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);
end
%% set a tolerance level
tol = 1e-10;

%% 1st: KNOWN STEADY STATE VALUES
% ANALYTICALLY 
LAM     = BETA;         % FIRMS DISCOUNT FACTORS EQUAL BETA IN SS
R       = 1/BETA;       % FROM EULER
r       = R-1;          % NET INTEREST RATE
PI      = 1;            % PRICES ARE EQUAL IN SS
pistar  = 1;            % PRICES ARE EQUAL IN SS
X       = 1;            % LAND IS FIXED
rk      = R-1+DK;       % Rental rate on capital is net interest + depreciation premium
rkh     = R - 1+DKH ;   % Rental rate on capital is net interest + depreciation premium
Mp_SS   = 1.15;          % Markup of 10 pct (implicitly epsilon = 11)

Ac      = 0;
Ah      = 0;
eR      = 0;

%% AUXILIARY PARAMETERS
% FROM RENTAL RATES
ZETA0   = MUC*BETA/(1-BETA*(1-DK))/Mp_SS; 
ZETA1   = MUH*BETA/(1-BETA*(1-DKH));

% FROM HOUSING FOC
ZETA2   = JEI / (1-BETA*(1-DH));
ZETA3   = JEI / (1-BETA1*(1-DH)-(BETA-BETA1)*LTV);

% FROM BORROWING CONSTRAINT
ZETA4   = (R-1)*LTV/R;

%% SOLVER
x0 = [7; 1; 1];    %// Initial guess for fsolve
[x,fval] = fsolve(@ssfun,x0,[],OMEGAH,EPS,ZETA3,BETA,DH) ;

H1 = x(2);
HZ = x(3);

%%
% DEFINITION OF HTHILDE
HTHILDE = (OMEGAH^(1/EPS)*H1^((EPS-1)/EPS)+(1-OMEGAH)^(1/EPS)*HZ^((EPS-1)/EPS))^(EPS/(EPS-1));
%%
CHI1    = 1+DH*ZETA2*(1-r*ZETA1-RHOH-ALPHAH*(1-MUH-RHOH)); % OK

CHI2    = DH*ZETA3*(1/HTHILDE)*(HTHILDE*OMEGAH/H1)^(1/EPS)*H1*(r*ZETA1+RHOH+ALPHAH*(1-MUH-RHOH)) +ZETA4*ZETA3*(1/HTHILDE)*H1*(HTHILDE*OMEGAH/H1)^(1/EPS)+HZ*(HTHILDE*(1-OMEGAH)/HZ)^(1/EPS); %FSOLVE

CHI3    = (Mp_SS-1+r*ZETA0*Mp_SS+ALPHAC*(1-MUC))/Mp_SS; % OK

CHI4    = 1 + DH*ZETA3*(1/HTHILDE)*(HTHILDE*OMEGAH/H1)^(1/EPS)*H1*(1-(1-ALPHAH)*(1-MUH-RHOH))+ZETA4*ZETA3*(1/HTHILDE)*(HTHILDE*OMEGAH/H1)^(1/EPS)*H1+(HTHILDE*(1-OMEGAH)/HZ)^(1/EPS)*HZ; %FSOLVE

CHI5    = (1-ALPHAH)*(1-MUH-RHOH)*DH*ZETA2; % OK

CHI6    = (1-MUC)*(1-ALPHAC)/Mp_SS; % OK


%% RATIOS
CtoY    = (CHI3*CHI4+CHI2*CHI6)/(CHI1*CHI4-CHI2*CHI5);

C1toY   = (CHI1*CHI6+CHI3*CHI5)/(CHI1*CHI4-CHI2*CHI5);

QIHtoY   = DH*ZETA2*CtoY + DH*ZETA3*(HTHILDE*OMEGAH/H1)^(1/EPS)*H1;


%% LABOR
% AUX RATIOS
RATIOL = (1-MUH-RHOH)*ALPHAH/((1-MUC)*ALPHAC)*Mp_SS*QIHtoY;
RATIOL1 = (1-MUH-RHOH)*(1-ALPHAH)/((1-MUC)*(1-ALPHAC))*Mp_SS*QIHtoY;

% HOURS WORKED
LC  = ( ((1-TAUW)*(1-MUC)*ALPHAC/CtoY/Mp_SS)/((1+RATIOL)^(ETA-1)) )^(1/ETA);
LH  = RATIOL*LC;

LC1 = (((1-TAUW)*(1-MUC)*(1-ALPHAC)/C1toY/Mp_SS)/((1+RATIOL1)^(ETA-1)))^(1/ETA) ;
LH1 = RATIOL1*LC1;

%% PRODUCTION
Y   = (LC^ALPHAC)*(LC1^(1-ALPHAC))*ZETA0^(MUC/(1-MUC)) ;

IH  = LH^(ALPHAH*(1-MUH-RHOH)) * LH1^((1-ALPHAH)*(1-MUH-RHOH)) * ZETA1^MUH*(QIHtoY*Y)^MUH ;

qh  = QIHtoY*Y/IH;
QIH = QIHtoY / Y;
qrent = (1-BETA*(1-DH))*qh;

%% CAPITAL
kc  = ZETA0*Y;
kh  = ZETA1*Y;
IK  = DK*kc + DKH*kh;

%% CONSUMPTION
C   = CtoY * Y;
C1  = C1toY * Y;

%% HOUSING 
H   = ZETA2*C/qh;
%H1  = (ZETA3*C1/qh)^EPS;

%% OTHER
% BORROWING
b   = LTV*qh*H1/R; 

% MULTIPLIER
PSI = (BETA-BETA1)/C1;

% MARKUP
Mp = VAREPS/(VAREPS-1); % VAREPS SHOULD BE 11

% GDP
GDP = qh*IH + C + C1 + IK;

% 
rx = qh*RHOH*IH;

% FROM PATIENT FOC
qx = rx/(1-BETA);

% WAGES FROM FIRM FOC
wc = (1-MUC)*(ALPHAC)*Y/Mp/LC;
wc1 = (1-MUC)*(1-ALPHAC)*Y/Mp/LC1;

wh = (1-MUH-RHOH)*ALPHAH*QIH/LH;
wh1 = (1-MUH-RHOH)*(1-ALPHAH)*QIH/LH1;

% PRICE AUX
% 26) - Auxiliary price setting recursion 1
z1 = THETA*LAM/Mp*Y;

% 27) - Auxiliary price setting recursion 2
z2 = THETA*LAM*Y;

% 28) - Reset inflation
pistar = (VAREPS/(VAREPS-1))*z1/z2;

%% STEADY STATE MATRIX
 ys = [LH1 IK IH kc kh Y GDP X b  qh qrent qx ...
R rk rx rkh wc wc1 wh wh1 PI Mp Ac Ah eR  z1 z2 pistar LAM PSI]';


%% Step 4: Update parameters and variables
params=NaN(M_.param_nbr,1);

%%
for iter = 1:M_.param_nbr %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end
%%
for ii = 1:M_.orig_endo_nbr %auxiliary variables are set automatically
  eval(['ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
end
 %%
save('Steady_state_external.mat','ys');

%save('calculatedSS_DSGE1.mat', 'ys');
end
