RBC code simulation error

Hello guys,

I am trying to simulate RBC model. I have the dynare code below:
var c k y i w Rent n r a;
varexo e;
parameters alpha beta delta sigma teta rho sigmae;
beta=0.99;
alpha=0.36;
sigma=2;
delta=0.025;
sigmae = 0.01;
rho=0.95;
teta=1;
model;
1/exp(c(-1))=beta*((1/(exp©))((aalphaexp(k)^(alpha-1)exp(n)^(1-alpha))+(1-delta)));
teta/(1-exp(n))=(1/exp©)
(exp(k)^(alpha)a((1-alpha)exp(n)^(-alpha)));
exp(k)=(a
exp(k(-1))^(alpha)exp(n(-1))^(1-alpha))-exp©+((1-delta)exp(k(-1)));
exp(y)=(a
exp(k(-1))^(alpha)exp(n)^(1-alpha));
exp(i)=exp(y)-exp©;
1/exp(c(-1))=beta
(1/exp©)
(1+exp®);
exp(w)=(a
exp(k)^(alpha)((1-alpha)exp(n)^(-alpha)));
exp(Rent)=(a
alpha
exp(k)^(alpha-1)exp(n)^(1-alpha));
a=rho
a(-1)+e;
end;
initval;
w =log(1);
Rent =log(0.1);
r =log(0.002);
n =log(0.6);
i =log(0.6);
y =log(1.5);
c =log(1.5);
k =log(16);
end;
shocks;
var e=sigmae^2;
end;
steady;
stoch_simul(hp_filter=1600,order=1,irf=40);

I am getting the following error:

Error using lnsrch1 (line 53)
Some element of Newton direction isn’t finite. Jacobian maybe
singular or there is a problem with initial values

Could you help me to find the error in my code.

Thanks

There are a bunch of errors. Your TFP in the equations should be exp(a) as otherwise output in steady state is 0 as a has steady state 0. Moreover, your timing is wrong. For example, your Euler equation is not forward looking.

Thanks jpfeifer for your help, I considered your corrections and it worked. Do you think the following code is accepted as standard RBC model simulation.

var c k y i w Rent n r a;
varexo e;
parameters alpha beta delta sigma teta rho sigmae;
beta=0.99;
alpha=0.36;
sigma=2;
delta=0.025;
sigmae = 0.01;
rho=0.95;
teta=1;
model;
1/exp©=beta*((1/(exp(c(+1))))*((exp(a(+1))alphaexp(k(+1))^(alpha-1)exp(n(+1))^(1-alpha))+(1-delta)));
teta/(1-exp(n))=(1/exp©)
(exp(k)^(alpha)exp(a)((1-alpha)*exp(n)^(-alpha)));
exp(k)=(exp(a)*exp(k(-1))^(alpha)*exp(n(-1))^(1-alpha))-exp©+((1-delta)*exp(k(-1)));
exp(y)=(exp(a)exp(k(-1))^(alpha)exp(n)^(1-alpha));
exp(i)=exp(y)-exp©;
1/exp(c(-1))=beta
(1/exp©)
(1+exp®);
exp(w)=(exp(a)exp(k)^(alpha)((1-alpha)*exp(n)^(-alpha)));
exp(Rent)=(exp(a)alphaexp(k)^(alpha-1)exp(n)^(1-alpha));
exp(a)=rho
exp(a(-1))+e;
end;
initval;
w =1;
Rent =1;
r =1;
n =1;
i =1;
y =1;
c =1;
k =1;
a=0;
end;
shocks;
var e=sigmae^2;
end;
steady;
stoch_simul(hp_filter=1600,order=1,irf=40);

Many thanks,

The timing for capital is wrong in some equations and correct in others.

Here is my code, I spent one two days by trying to find the solution. Someone can 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);