function [ys,check] = NDF_steadystate(ys,exo)
% function [ys,check] = NDF_steadystate(ys,exo)
% computes the steady state for the NDF.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;


%% Enter model equations here

options = optimset('Display','iter','TolFun',1e-16); % set options for numerical solver

% the steady state computation follows 
pi =1.15;
A =1;
Y_oil =1;
G_c =1;
mu =eta*pi;
i =((pi*eta)/betta)-1;
r_kp =((1+i)/pi)+ delta_p -1;
mc =(1/theta)*((1-betta*eta^2)* phi_p *(pi-1)* pi + theta-1);

[W,fval]=fzero(@(W) (1-(eta+delta_p-1)*(alpha_kp*mc/r_kp)-(phi_p/2)*(pi-1)^2)...
+((W/alpha_lg)^(alpha_kg/alpha_lg)-(eta+delta_p-1)*(1+phi*(delta_g/2)))*(W/(alpha_lp*mc))^(alpha_lp/alpha_kgp)*(r_kp/(alpha_kp*mc))^(alpha_kp/alpha_kgp)-(W/psi_l) ,0.2,options);

C = W / psi_l;
m =((1+i)/i) * C* psi_m;
K_g = (W/(alpha_lp*mc))^(alpha_lp/alpha_kgp)*(r_kp/(alpha_kp*mc))^(alpha_kp/alpha_kgp);
Y_g =  (W/alpha_lg)^(alpha_kg/alpha_lg)* K_g;
Y_p = (W/(alpha_lp*mc))^(alpha_kgp/alpha_lp)*(r_kp/(alpha_kp*mc))^(alpha_kgp/alpha_kp)* K_g;
K_p = ( alpha_kp * mc / r_kp )* Y_p;
I_ge = (eta + delta_g - 1)* K_g;
I_g = I_ge * ( 1 + phi * ( delta_g / 2 ));

L_g = ( alpha_lg / W )* Y_g;
L_p = ( alpha_lp * mc / W )* Y_p;

I_pa = (eta + delta_p - 1) * K_p;
L =  L_p + L_g;
I = I_g + I_pa;
Y_noil = Y_p + Y_g;
Y = Y_noil + 1;
 
NDF = phi_f / ((1-(1/eta)) - r -(((1+rd)*alpha_nd /(1+alpha_nd-1/eta))-1) * alpha_f); 
F = alpha_f * NDF ;
Z = r * NDF;   
ND = ((1+rd)/(1+alpha_nd-1/eta)) * F;
I_p = I_pa - F ;
mubar = mu;
pibar = pi;
Y_pbar = Y_p;

%% end own model equations

for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
