function [yss,check] = cn_aging_debt20230109_steadystate(yss, exo)
% function [ys,check] = alpha_steadystate(ys,exo)
% computes the steady state for the alpha.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;

alpha=0.4;
%% 可直接求


%% 给定k求出其余
x0=[1.05];
options = optimoptions('fsolve','TolFun',1e-15, 'MaxIter', 1e5);
[xxx]=fsolve(@kss,x0,options,omega_W, omega_r, omega_d, omega_y, betaH, sigmaH, kai, alpha,  delta, tau, ss, tau0,  gcss);
R=xxx(1);


g_n=omega_y +1-omega_W;
ns_W=1/(1+(g_n -(1 - omega_r))/omega_W+omega_r/(g_n - (1 - omega_d)));
ns_y=(g_n -(1 - omega_r))/omega_W*ns_W;
ns_r=omega_r*ns_W/(g_n - (1 - omega_d));
% R=g_n/betaH;
r_k=  delta + R-1;
xi_r=1 -g_n^(kai*sigmaH/(1 - sigmaH))*(betaH*R)^(1/(1 - sigmaH))*((1 - omega_d)/R);
xi_w0=[0.1];
options = optimoptions('fsolve','TolFun',1e-15, 'MaxIter', 1e6);
[xx]=fsolve(@xi_wss,xi_w0,options,omega_r,xi_r,sigmaH,betaH,R,g_n,kai);
xi_w=xx(1); %%%%%%
omega= 1 - omega_r + omega_r*(xi_r/xi_w)^((sigmaH - 1)/ sigmaH);


yk=r_k/( alpha*(1-tau)*g_n);
k=(yk/(ns_W^(1-alpha))*g_n^alpha)^(1/(alpha-1));
y=yk*k;
W = (1-alpha)*(1-tau)*y/ns_W;
I=delta*k;
tau_W = tau0*y/ns_W;
H_W = W*ns_W/( 1- (1 - omega_r)/(R*omega));
T_W = tau_W*ns_W/(1- (1 - omega_r)/(R*omega));
sss=ss*y;
s_r=sss/ns_r;
S_r = s_r*ns_r/(1- (1 - omega_d)/R);
S_W = ((omega-(1-omega_r))/(R*omega)*S_r)/(1-(1 - omega_r)/(R*omega) );
% kg=(y/((k/g_n)^alpha*ns_W^phi))^(1/(1-alpha-phi))*g_n;
% gi=delta*kg;
PiF=0;
d_W = 0;
d_r = 0;
D_r = d_r*ns_r /(1- (1 - omega_d)/R);
D_W =( d_W*ns_W + (omega-(1-omega_r))/(R*omega)*D_r)/(1-(1 - omega_r)/(R*omega));
A_W=((1 - omega_r)*(d_W*ns_W + W*ns_W- tau_W*ns_W- xi_w*( D_W + H_W+S_W - T_W)))/(1-(1 - omega_r)*(R/g_n-xi_w*R/g_n));
C_W = xi_w*(R*A_W/g_n + D_W + H_W+S_W - T_W);
A_r= (d_r*ns_r+ s_r*ns_r - xi_r*( D_r+S_r) + omega_r*A_W/(1-omega_r))/(1-R/g_n+xi_r*R/g_n);
C_r = xi_r*(R*A_r/g_n + D_r+S_r);
bb=A_r+A_W - k;
gc=bb+tau*y-R*bb/g_n-sss;


%% 

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(['yss(' int2str(ii) ') = ' varname ';']);
end
