Solve for steady state

Here is my code. I need a help.

%%% Welfare-Maximizing Monetary and Fiscal Policy Rules in Closed Economy%%%
var c gc ig y b w l r rk theta gk k mc tauk tauc taul t pi pstar n j ip def;
varexo e_gc e_ig e_theta;

parameters beta, alpha, omega, gamma, psi, delta, r_ss, tauc_ss, tauk_ss, l_ss, taul_ss, b_ss, y_ss, theta_ss, gc_ss, ig_ss, pi_ss, eta, nu, rhogc, rhoig, rhotheta, rhopi, rhoy, rhodef, phikb, phikgc, phikig, phiktheta,
phicb, phicgc, phicig, phictheta, philb, philgc, philig, philtheta;

beta = .99;
alpha = .24;
omega = .38;
gamma = .39;
delta = 0.025;
r_ss = 1.1;
tauc_ss = 0.05;
tauk_ss = 0.15;
l_ss = .33;
taul_ss = 0.26;
b_ss = 0.5;
theta_ss = 1;
y_ss = 1.12923626;
gc_ss = 0.1179;
ig_ss = 0.0996;
pi_ss = 1;
eta = 0.75;
nu = 6;
rhogc = 0.95;
rhoig = 0.95;
rhotheta = 0.95;
rhopi = 1.53;
rhoy = 0;
rhodef = 0;
phikb = 8.56;
phikgc = 0;
phikig = 0;
phiktheta = 0;
phicb = 8.56;
phicgc = 0;
phicig = 0;
phictheta = 0;
philb = 8.56;
philgc = 0;
philig = 0;
philtheta = 0;
b_ss = 0;
gk_ss = (ig_ss/delta);
mup = (nu/(nu-1));
mc_ss = (1/mup);
r_ss = (1/beta);
pstar_ss = ((1-etapi_ss^(nu-1))/(1-eta))^(1/(1-nu));
k_ss = (alpha
mc_sstheta_ss(gk_ss^omega)(l_ss^(1-alpha)))/((r_ss - (deltatauk_ss) - (1-delta))/(1-tauk_ss));
rk_ss = alphamc_sstheta_ss*(gk_ss^omega)(l_ss^(1-alpha))(k_ss^(alpha-1));
ip_ss = deltak_ss;
w_ss = (1-alpha)mc_sstheta_ss
(gk_ss^omega)(k_ss^alpha)(l_ss^(-alpha));
n_ss = (1/(1-etabeta))theta_ss(gk_ss^omega)(k_ss^alpha)(l_ss^(1-alpha))mc_ss;
j_ss = (1/(1-eta
beta))theta_ss(gk_ss^omega)
(k_ss^alpha)(l_ss^(1-alpha))pstar_ss;
j_ss = (nu/(nu-1))n_ss;
y_ss = theta_ss
(gk_ss^omega)
(k_ss^alpha)
(l_ss^(1-alpha));
c_ss = y_ss - ip_ss - ig_ss - gc_ss;
psi = (1/(c_ss + gammagc_ss))((1-taul_ss)/(1+tauc_ss))w_ss;
t_ss = gc_ss + ig_ss;
t_ss = tauc_ss
c_ss + tauk_ssrk_ssk_ss - deltatauk_ssk_ss + taul_ssw_ssl_ss;
def_ss = b_ss;

model;

exp© = exp(y) - exp(ip) - exp(ig) - exp(gc);
exp(k) = exp(ip) + exp(k(-1));
exp(gk) = exp(ig) + exp(gk(-1));
beta*(((exp© + gammaexp(gc))/(exp(c(+1)) + gammaexp(gc(+1))))((1 + exp(tauc))/(1 + exp(tauc(+1))))(1/exp(pi(+1)))exp®) = 1;
beta
(((exp© + gammaexp(gc))/(exp(c(+1)) + gammaexp(gc(+1))))((1 + exp(tauc))/(1 + exp(tauc(+1))))(((1-exp(tauk(+1)))exp(rk(+1)) + (deltaexp(tauk(+1))) + (1-delta)))) = 1;
psi = (1/(exp© + gammaexp(gc))) ((1 - exp(taul))/(1 + exp(tauc)))exp(w);
exp(y) = exp(theta)
(exp(gk)^omega)
(exp(k)^alpha)
(exp(l)^(1-alpha));
exp(rk) = alphaexp(mc)exp(y)/exp(k);
exp(w) = (1-alpha)exp(mc)(exp(y)/exp(l));
exp(j) = (nu/(nu-1))exp(n);
exp(pstar) = ((1 - eta
(exp(pstar)^(nu - 1)))/(1 - eta))^(1/(1-nu));
exp(n) = exp(y)exp(mc) + (etabeta
((exp© + gamma
exp(gc))/(exp(c(+1)) + gammaexp(gc(+1))))(exp(pi(+1)^nu)exp(n(+1))));
exp(j) = exp(y)exp(pstar) + (etabeta
((exp© + gammaexp(gc))/(exp(c(+1)) + gammaexp(gc(+1))))(exp(pstar)/exp(pstar(+1)))(exp(pi(+1)^(nu-1))exp(j(+1))));
exp(gc) + exp(ig) + exp(r(-1))b(-1) = b + exp(t);
exp(t) = exp(tauc)exp© + exp(tauk)exp(rk)exp(k(-1)) - deltaexp(tauk)exp(k(-1)) + exp(w)exp(l)exp(taul);
exp(gc - STEADY_STATE(gc)) = exp(rhogc
(gc(-1) - STEADY_STATE(gc)) + e_gc);
exp(ig - STEADY_STATE(ig)) = exp(rhogc
(ig(-1) - STEADY_STATE(ig)) + e_ig);
exp(theta - STEADY_STATE(theta)) = exp(rhotheta
(theta(-1) - STEADY_STATE(theta)) + e_theta);
def = b - b(-1);
exp(tauc - STEADY_STATE(tauc)) = phicb
(b - STEADY_STATE(b)) + exp((phicgc
(gc - STEADY_STATE(gc))) + phicig
(ig - STEADY_STATE(ig)) + (phictheta
(theta - STEADY_STATE(theta))));
exp(tauk - STEADY_STATE(tauk)) = phikb*(b - STEADY_STATE(b)) + exp((phikgc*(gc - STEADY_STATE(gc))) + phikig*(ig - STEADY_STATE(ig)) + (phiktheta*(theta - STEADY_STATE(theta))));
exp(taul - STEADY_STATE(taul)) = philb*(b - STEADY_STATE(b)) + exp((philgc*(gc - STEADY_STATE(gc))) + philig*(ig - STEADY_STATE(ig)) + (philtheta*(theta - STEADY_STATE(theta))));
exp(r - STEADY_STATE®) = exp((rhopi*(pi - STEADY_STATE(pi))) + (rhoy*(y - STEADY_STATE(y)))) + (rhodef*(def - STEADY_STATE(def)));

end;

initval;

c = log(c_ss);
gc = log(gc_ss);
ig = log(ig_ss);
y = log(y_ss);
b = b_ss;
w = log(w_ss);
l = log(l_ss);
r = log(r_ss);
rk = log(rk_ss);
theta = log(theta_ss);
gk = log(gk_ss);
k = log(k_ss);
mc = log(mc_ss);
tauk = log(tauk_ss);
tauc = log(tauc_ss);
taul = log(taul_ss);
pi = log(pi_ss);
pstar = log(pstar_ss);
n = log(n_ss);
j = log(j_ss);
ip = log(ip_ss);
def = def_ss;

end;

steady;
check;
shocks;
var e_gc = 0.01;
var e_ig = 0.01;
var e_theta = 0.01;

end;

options_.printing=1;
stoch_simul(order=2, irf=25) c y l ig gc pi gk k; %, hp_filter = 1600);

Put resid(1); before steady to see the residuals. You need to provide better starting values/correct steady state values.

Thanks. Let me try again.

Thanks, it works perfectly.