are inconsistent with r=1/beta from the Euler equation. Also, the last equation implies a unit root in the price level, meaning there are infinitely many steady states.

Sir，I have corrected the last part of the code, but there is a problem in the steady-state solution part, so that the result is still not obtained.In particular, there is a contradiction when manually solving the equilibrium solution, but the same method can be obtained and operated normally from other journals, phi variables and D variables.

var K I W C L r Q Yi A Z y Pi P phi RZ D RD tau pi B SLF NA SR;
varexo epsilon epsilon2 epsilonR;
parameters delta theta beta mu rho rho_r phi_pi rho_1 rho_phi theta_r Re rho_s r_ss gamma_R gamma_pi gamma_y pi_ss y_ss SR_ss;
delta = 0.06;
theta = 0.5;
beta = 0.99;
mu = 0.5;
rho = 6;
rho_r = 1.5;
phi_pi = 0.9;
rho_phi = 0.9;
theta_r = 0.1;
rho_1 = 0.68;
Re = 1.004;
gamma_R = 0.8;
gamma_pi = 1.5;
gamma_y = 0.5;
rho_s = 0.9;
r_ss = 0.0315;
pi_ss = 0.9263;
y_ss = 0.0737;
SR_ss = 0.0466;
model;
K(+1) = I + (1 - delta) * K;
W = ((1-theta) / theta) * C / (1 - L);
r = 1/beta * C(+1)/C;
C = C(-1) * beta * (1 - delta + Q);
Yi = A * K^mu * L^(1-mu);
Z = W * L + I;
L = (1 - mu) * Yi/W;
Q = mu * A * K^(mu-1) * L^(1-mu);
y = (Yi^((rho-1)/rho))^(rho/(rho - 1));
Yi = (Pi/P)^(-rho) * y;
P = (Pi^(1-rho))^(1/(1-rho));
B = (1 - theta_r) * D + SLF - Z;
D = (1 - phi) * D(-1);
phi = B/D;
log(phi) = rho_phi * log(phi(-1)) + epsilon;
RD = theta_r * Re + (1 - theta_r) * (1 - phi) * r;
RZ = r + tau/phi;
SLF = B + theta_r * D - NA;
r = r_ss * (r(-1)/r_ss)^(gamma_R) * ((pi/pi_ss)^gamma_pi * (y/y_ss)^gamma_y)^(1-gamma_R) * SR;
log(SR) = (1 - rho_s) * log(SR_ss) + rho_s * log(SR(-1)) + epsilonR;
pi = P/P(-1);
y = C + I + Z * tau/phi;
log(A) = rho_1 * log(A(-1)) + epsilon2;
end;
initval;
y = 1;
Yi = 1;
Pi = 1;
C = 29/62;
P = 1;
Q = 31/100;
W = 25/31;
r = 1.25;
A = 1;
K = 50/31;
I = 3/31;
Z = W * L + I;
D = 1;
phi = 1;
RD = theta_r * Re + (1 - theta_r) * (1 - phi) * r;
RZ = r + tau/phi;
SLF = B + theta_r * D - NA;
tau = (RZ - r) * phi;
pi =1;
epsilon = 0;
epsilon2 = 0;
epsilonR = 0;
end;
steady;
check;
shocks;
var epsilon;
periods 1:9;
values 0.05;
end;
check;
model_diagnostics;
stoch_simul(order = 1,irf = 40);