function yy=DJ_ss(xx,sigma_c,sigma_l,beta_t,beta,delta_k,delta_d,theta_p,alpha,theta,z,eta,epsilon_c, rhod,  rhog,rhoh,rhos)
%% 
k= xx(1);
N= xx(2);

beta_t = 0.99;
alpha =1/3;
beta=0.45;
theta = 0.11;
%% 
%% Enter model equations here
d=1;
h = 1; 
g = 1; 
s =1; 
pi=1;
i=1/beta_t-1;
r=i+1;
phi_c=epsilon_c;
%% 
I= delta_k*k;
pistar=((pi^(1-eta)-theta_p)/(1-theta_p))^1/(1-eta);
mc =(1-theta_p*beta_t*pi^eta)/(1-theta_p*beta_t*pi^(eta-1))*pi/pistar*(eta-1)/eta;
R =1/beta_t-(1- delta_k);
dp= (1-theta_p)*(pi/pistar)^eta/(1- pi^eta*theta_p);
p_d=(h+g+(theta-alpha-beta)*s)/(1-alpha-beta+theta);
D=((p_d-s)/((1-alpha-beta)*d*mc*k^alpha*N^beta))^(1/(theta-alpha-beta));
A=d*D^theta;
Y= A*k^alpha * N^beta * D^(1-alpha-beta)/dp;
c=Y-I;
w=N/(phi_c*c^(-sigma_c));
x1 = c^(-sigma_c)*phi_c*mc*Y/(1-beta_t*theta_p*pi^eta);
x2=c^(-sigma_c)*phi_c*Y/(1-beta_t*theta*pi^(eta-1));

%% 残差
f1=N-(c^-sigma_c*w*phi_c)^(1/sigma_l);
f2=k-w*alpha/(R*beta)*N;
yy=[f1,f2];



   
