function [yss,check] = cc_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
% beta					贴现因子
% 1/yita                 跨期替代弹性
% xi_l		xi_h			折旧率
% // alpha_l	alpha_h		资本在总产出中所占的比重
% // sigma_l	sigma_h		资本与劳动的替代弹性
% // gamma_l	gamma_h      数字经济对资本扩展的影响程度
% // delta_l	delta_h      数字经济对劳动扩展技术的影响程度
% // rho_l		rho_h        技术前面的系数
% mu_h mu_l 计算的中间量
%M_h M_l 计算的中间量
%A_lss A_l的稳态
%A_hss A_h的稳态
check = 0;
beta=0.99;	
yita=1/0.1526;
xi_l=0.05;
xi_h=0.05;
sigma_l =0.75;
sigma_h =0.5;
alpha_l=0.549;
alpha_h =0.535;
%epsilon = 0.284;
epsilon = 0.5; 
w1   = 0.329;
w2   = 0.671;

gamma_l =0.1;
gamma_h =0.2;
rho_l =0.95;
rho_h =0.95;
A_l=1;
A_h=1;
H=1;
chi=4.4;
delta_l=-0.01;
delta_h= 0.01;


r_h  =1/beta+xi_h-1;
r_l  =1/beta+xi_l-1;

x0=[1;1;0.0857138;0.115864];
[xx,fval]=fsolve(@qumodelNfunction,x0,optimoptions('fsolve','TolFun',1e-19, 'MaxIter', 1e12),alpha_h,alpha_l,sigma_h,sigma_l,gamma_h,gamma_l,A_h,A_l,H,w1,w2,xi_h,xi_l,delta_h,delta_l,yita,chi,r_h,r_l,epsilon);

W_h=xx(1);
W_l=xx(2);
L_h=xx(3);
L_l=xx(4);

C_l=(chi/W_l)^(-1/yita);
C_h=(chi/W_h)^(-1/yita);

lambda_h = C_h^(-yita);
lambda_l = C_l^(-yita);

K_h=(alpha_h*W_h/(1-alpha_h)/r_h)^sigma_h*H^((gamma_h-delta_h)*(sigma_h-1))*L_h;
K_l=(alpha_l*W_l/(1-alpha_l)/r_l)^sigma_l*H^((gamma_l-delta_l)*(sigma_l-1))*L_l;

Y_h = A_h * ( alpha_h * (H^gamma_h * K_h)^((sigma_h - 1)/sigma_h) + (1-alpha_h)*(H^delta_h * L_h)^((sigma_h-1)/sigma_h) )^(sigma_h/ (sigma_h - 1));
Y_l = A_l * ( alpha_l * (H^gamma_l * K_l)^((sigma_l - 1)/sigma_l) + (1-alpha_l)*(H^delta_l * L_l)^((sigma_l-1)/sigma_l) )^(sigma_l / (sigma_l - 1));


Y=(w1*Y_h^((epsilon-1)/epsilon)+w2*Y_l^((epsilon-1)/epsilon))^(epsilon/(epsilon-1));

P_h = w1*(Y/Y_h)^(1/epsilon);
P_l = w2*(Y/Y_l)^(1/epsilon);

%K= K_h + K_l;
%L= L_h + L_l;



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
