clear all;

beta=0.998;
deltak=0.025;
deltah=0.005;
xi_ss=1;
psi_ss=1;
varphi=1;
alpha=0.5;
mu=0.2;
theta=0.7;
A_ss=1;
J_ss=1;
taus_ss=0.05;
tauy_ss=0.15;
tauw_ss=0.15;
G_ss=1;
Z_ss=1;
news_ss=1;

R_ss=1/beta-(1-deltak);

w_ss=(1-alpha)/alpha*(alpha*(1-tauy_ss)*G_ss^mu/R_ss)^(1/(1-alpha))*R_ss;

F= @(x) [xi_ss*x(1)-(1-beta*(1-deltah))*(1+taus_ss)*x(4)*x(2);
         (1-tauw_ss)*w_ss-psi_ss*x(3)^varphi*x(1);
         (1-alpha)*x(6)*R_ss-alpha*w_ss*x(3);
         x(7)-(theta*(1-taus_ss)*x(4)/R_ss)^(1/(1-theta))*Z_ss;
         (1-theta)*x(7)*R_ss-theta*x(8)*Z_ss;
         deltah*x(2)-x(7)^theta*Z_ss^(1-theta);
         x(6)^alpha*x(3)^(1-alpha)*G_ss^mu-x(1)-deltak*x(5)-G_ss;
         x(5)-x(6)-x(7);
        ];

x0 = [0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5];

options = optimoptions('fsolve','Display','iter','MaxFunctionEvaluations',10000,'MaxIterations',1000);

[x,fval] = fsolve(F,x0,options)

C_ss=x(1);
H_ss=x(2);
L_ss=x(3);
Ph_ss=x(4);
Kh_ss=x(5);
Kp_ss=x(6);
Ke_ss=x(7);
Pz_ss=x(8);

Y_ss=A_ss*Kp_ss^alpha*L_ss^(1-alpha)*G_ss^mu;

Yh_ss=J_ss*Ke_ss^theta*Z_ss^(1-theta);

q_ss=1/C_ss;

I_ss=deltak*Kh_ss;

eta_sy_ss=log(taus_ss)/log(Y_ss+Ph_ss*Yh_ss);

eta_wy_ss=log(tauw_ss)/log(Y_ss+Ph_ss*Yh_ss);

eta_yy_ss=log(tauy_ss)/log(Y_ss+Ph_ss*Yh_ss);

RE_ss=tauw_ss*w_ss*L_ss+tauy_ss*Y_ss+2*taus_ss*Ph_ss*Yh_ss+Pz_ss*Z_ss;

sur_ss=(RE_ss-G_ss)/(Y_ss+Ph_ss*Yh_ss);

rely_ss=(taus_ss*Ph_ss*Yh_ss+Pz_ss*Z_ss)/RE_ss;



