Can't find a steady state, please help!

[miles]

Please help with my code below, I keep getting the following error:

??? Error using ==> print_info at 57
Impossible to find the steady state. Either the model doesn't have a
steady state, there are an infinity of steady states, or the guess values
are too far from the solution

Error in ==> steady at 92
print_info(info,options_.noprint);

Error in ==> NKHabit1 at 252
steady;

Error in ==> dynare at 120
evalin('base',fname) ;

My model code as follows:

var y c I g a k_hat k n u nu pi pi_hush A_hat D_hat mc w R lambda mu theta m r q y_star;

varexo epsilon_a epsilon_theta epsilon_g epsilon_r epsilon_y;

parameters alpha epsilon phi beta gamma si delta sigma_o tau omega rho_a rho_g rho_r rho_theta phi_pi phi_y v theta_star pi_star rho_y;

alpha=0.3;
epsilon=6;
phi=0.5;
beta=0.99;
gamma=0.75;
si=1;
delta=1.4040;
sigma_o=0.025;
tau=0.5;
pi_star=0.005;
omega=0.2;
rho_a=0.97;
rho_g=0.85;
rho_r=0.85;
rho_theta=0.85;
phi_pi=1.5;
phi_y=0.5;
theta_star=1.9358;
v=1;
rho_y=0.95;

model;
exp(y)=exp(c)+exp(I)+exp(g);
exp(y)=(exp(a)*(exp(k_hat)^alpha)*exp(n)^(1-alpha))/exp(nu);
exp(k_hat)=exp(u)*exp(k);
exp(pi)=(((1-phi)*(exp(pi_hush)^(1-epsilon))+phi))^(1/(1-epsilon));
exp(nu)=(1-phi)*((exp(pi_hush)^(-epsilon))*(exp(pi)^epsilon))+phi*(exp(pi)^epsilon)*exp(nu(-1));
exp(pi_hush)=exp(pi)*(epsilon/(epsilon-1))*exp(A_hat)/exp(D_hat);
exp(A_hat)=exp(lambda)*exp(y)*exp(mc)+phi*beta*(exp(pi(+1))^epsilon)*exp(A_hat(+1));
exp(D_hat)=exp(lambda)*exp(y)+phi*beta*(exp(pi(+1))^(epsilon-1))*exp(D_hat(+1));
exp(w)=exp(mc)*(1-alpha)*((exp(k_hat)/(exp(n)))^alpha);
exp(R)=exp(mc)*alpha*((exp(k_hat)/(exp(n)))^(alpha-1));
exp(lambda)=(1/(exp(c)-gamma*exp(c(-1))))-((beta*gamma)/(exp(c(+1))-gamma*exp(c)));
exp(lambda)*exp(w)=theta*((1-exp(n))^(-si));
exp(lambda)*exp(R)=exp(mu)*delta*sigma_o*(exp(u)^(delta-1));
exp(lambda)=beta*exp(lambda(+1))*(1+exp(r(+1)))*(exp(pi(+1))^(-1));
exp(mu)=beta*(exp(lambda(+1)*exp(R(+1))*exp(u(+1))+exp(mu(+1))*(1-sigma_o*(exp(u(+1))^delta))));
exp(m)^(-v)=exp(lambda(+1))-beta*exp(lambda(+1))*(exp(pi(+1))^(-1));
exp(lambda)=exp(mu)*(1-((tau/2)*(((exp(I)/exp(I(-1)))-1)^2))-(tau*((exp(I)/exp(I(-1)))-1)*exp(I)/exp(I(-1))))+(beta*exp(mu(+1))*tau*((exp(I(+1))/exp(I))-1)*((exp(I(+1))/exp(I))^2));
exp(k(+1))=((1-((tau/2)*((exp(I)/exp(I(-1)))-1)^2))*exp(I))+(1-sigma_o*(exp(u)^delta))*exp(k);
exp(q)=exp(mu)/exp(lambda);
a=rho_a*(a(-1))+epsilon_a;
theta=(1-rho_theta)*theta_star+rho_theta*theta(-1)+epsilon_theta;
g=(1-rho_g)*(omega+y_star)+rho_g*g(-1)+epsilon_g;
exp(r(+1))=rho_r*exp(r)+(1-rho_r)*phi_pi*(exp(pi)-pi_star)+((1-rho_r)*phi_y*((exp(y)/exp(y(-1)))-1))+epsilon_r;
y_star=rho_y*(y_star(-1))+epsilon_y;
end;

initval;
y=0;
c=0.5;
I=0.17;
g=0.18;
a=1;
k_hat=7;
k=7;
n=0.3;
u=1;
nu=1;
pi=0.0005;
pi_hush=0.0266;
A_hat=7;
D_hat=9;
mc=0.8;
w=1.5;
R=0.035;
lambda=1.87;
mu=1.87;
theta=1.93;
m=35.6;
r=0.015;
q=1;
end;
steady;

shocks;
var epsilon_a; stderr 0.01;
var epsilon_theta; stderr 0.01;
var epsilon_g; stderr 0.01;
var epsilon_r; stderr 0.01;
var epsilon_y; stderr 0;
end;

stoch_simul(periods=1000, order=1, irf=20)y c I g a k_hat k n u nu pi pi_hush A_hat D_hat mc w R lambda mu theta m r q;

Thank you so mcuh!

[mghazias]

There is no way to help with this,
You need paper and pencil + time, to solve the model and find the steady state values, or at least close guess to SS values.

[selima]

excuse what do you mean by paper and pencil method to find steady state ? thank you for the answer.

[mutichai]

Let go back to see your results again.
I think that before this line

"Error using ==> print_info at 57
Impossible to find the steady state. Either the model doesn't have a
steady state, there are an infinity of steady states, or the guess values
are too far from the solution"

Matlab will give you the residual value of each equation.
To get steady state values, you need to provide initial values that make residual of each equation close to zero.
I suggest you to fix the equation that the residual value is still large.
Hope this help.