clear all;
close all;

%parameters
beta=0.99;
sigma=2;
phi=2;
delta=0.025;
alpha=1/3;
epsilon_p=11;
theta_p=0.75;
epsilon_wnr=11;
epsilon_wr=11;
theta_wr=0.75;
theta_wnr=0.75;
eta=1;
rho_r=0.8;
rhor_pi=1.7;
rhor_y=0.45;
rho_g=0.8;
rho_t=0.8;
rhot_g=0.5;
rhot_b=0.1;
rho_A=0.8;
omega=0.5;
gss=0.05;
BoverY=0.4;
%Steady state of variables
piss = 1;
pistarss = 1;
dss = 1;
Rbss = 1/beta-1;
Rkss = 1/beta-(1-delta);
mcss = (epsilon_p-1)/epsilon_p;
Ass = 1;
KNratio  = (Rkss/(mcss*alpha*Ass))^(1/(alpha-1));%K/N
w_rss = mcss*(1-alpha)*Ass*((KNratio)^alpha);
w_nrss = w_rss;
w_rstarss = w_rss;
w_nrstarss = w_nrss;
N_nrss = (((w_nrss)^(1-sigma))*((epsilon_wnr-1)/epsilon_wnr))^(1/(phi+sigma));
C_nrss = w_nrss*N_nrss;
DD = (1-gss)*Ass*(KNratio^alpha)-delta*KNratio;%auxiliary parameters for simplicity
save afile 
%solve the equations of C_rss and N_rss
xx = fsolve(@(x)somevarsolve(x,DD,omega,N_nrss,C_nrss,w_rss,epsilon_wr,sigma,phi),[1;1]);
C_rss =xx(1);
N_rss = xx(2);
Nss = omega*N_rss+(1-omega)*N_nrss
Kss= KNratio*Nss;
Yss = Ass*(Kss^alpha)*(Nss^(1-alpha));
Gss = gss*Yss;
bss = BoverY*Yss;
tauss = bss*Rbss+gss*Yss;
iss = delta*Kss;
lamba_rss = C_rss^(-sigma);
lamba_nrss = C_nrss^(-sigma);
X1ss = ((C_rss)^(-sigma))*mcss*Yss/(1-theta_p*beta);
X2ss = ((C_rss)^(-sigma))*Yss/(1-theta_p*beta);
H1_rss = N_rss^(1+phi)/(1-theta_wr*beta);
H2_rss = lamba_rss*N_rss/(1-theta_wr*beta);
H1_nrss = N_nrss^(1+phi)/(1-theta_wnr*beta);
H2_nrss = lamba_nrss*N_nrss/(1-theta_wnr*beta);
Css = omega*C_rss+(1-omega)*C_nrss;
save steadystate %save the steadystate value of variables
dynare work