BGG code problem

Hello,

I am attempting to write a simple model code with BGG Accelerator. The paper I am trying to shadow is Christiansen and Dib 2005 (Bank of Canada). My code is given below. It does not work. I think the problem is in the steady state section. However, I have been able to recreate the steady state in a separate matlab file using this same model code. i.e. I used this code, took out the time indices and fed steady state values in the LHS, and the RHS was able to recreate the same steady state values.

Any idea why this .mod file is not working?

// Variable Declaration
//*************************

var e c ce n lambda y k sphi i rd rl rk rn w q nw s z pi sf;

varexo ez er ed;

parameters h gamma beta Upsilon zeta alpha delta nu chi theta phi mu rho sigm kappan
eta b psi srhop srhoy rhoe rhoz
lss nss piss rdss rnss iss qss sphiss rkss rlss wss xss ssigss
yss kss nwss css lambdass cess pss wnss sss mss uss vss vess srhoss
kss_yss nss_yss kss_nwss mss_nss wnss_pss cess_yss css_yss vss_uss;

// Parameter Values
//*************************

h = 0.65; // habit formation parameter
gamma = 2; // Risk aversion parameter
beta = 0.9979; // Discout factor
Upsilon = 1.013; // Weight on leisure in utility function
zeta = 1; // Inverse of the Frisch elasticity of labour supply
alpha = 0.33; // Capital share
delta = 0.025; // Capital depreciation rate
nu = 0.9833; // Probability that an entrepreneur survives till the next period
chi = 0.25; // Adjustment cost on investment
theta = 11; // Markup parameter
phi = 0.67; // Calvo parameter
//psi = 0.05; // Risk premium elasticity
psi = 0.01;

nss = 0.250346119815; // Steady State employment from my own .m file

srhop = 1.5; // Monetary policy response to inflation
srhoy = 0.125; // Monetary policy response to output gap.

rhoe = 0.8; // Persistence of demand shock
rhoz = 0.8; // Persistence of technology shock
sige = 0.0073;
sigz = 0.009;
sigr = 0.005;

// Steady State Values:
piss = 1;
qss = 1;
sphiss =(theta-1)/theta;
rdss = 1/beta;
rnss = 1/beta;
rkss = (1/nu)-(1-delta);
rlss = (1/nu);
kss_yss = alpha*(sphiss/rkss);
kss_nss = kss_yss^(1/(1-alpha));
yss_nss = kss_nss/kss_yss;
kss = kss_nssnss;
yss = yss_nss
nss;
wss = sphiss*(1-alpha)yss_nss;
nwss = kss
(nu/beta)^(1/psi);
cess = ((1-nu)/nu)nwss;
iss = delta
kss;
lambdass = (Upsilon/wss)(1-nss)^(-zeta);
css = yss-cess-iss;
lambdass2 = (1-beta
h)css^((1-h)(1-gamma)-1);
sss = (kss/nwss)^psi;
sfss = rlss*(kss-nwss);

// Nonlinear Model:
//*************************

model;
c=(1/lambda)(e(c/c(-1)^h)^(1-gamma)-betahe(+1)(c(+1)/c^h)^(1-gamma));
w=Upsilon/lambda
(1-n)^(-zeta);
rd=1/(beta*(lambda(+1)/lambda)(1/pi));
rn=1/(beta
(lambda(+1)/lambda)(1/pi));
w=sphi
(1-alpha)(y/n);
rk=alpha
sphi*(y/k(-1));
rl=(rk(+1)+q(+1)(1-delta))/q;
sf = rl
(qk-nw);
y=z
k(-1)^alpha*(n)^(1-alpha);
nw=nu*((rk+q*(1-delta))k(-1)-sf(-1));
rl(-1)=rd
s;
s=(qk/nw)^psi;
ce=((1-nu)/nu)nw;
q=1+chi
(i/k(-1)-delta);
k=i+(1-delta)k(-1);
pi=piss
((sphi/sphiss)^((1-phi)
(1-betaphi)/phi)(pi(+1)/piss)^beta);
y=c+ce+i+(chi/2)(i/k(-1)-delta)^2k;
rn=rnss*((pi/piss)^srhop*(y/yss)^srhoy*exp(er));

//shocks:

e = rhoee(-1)+ed;
z = rhoz
z(-1)+ez;

end;

initval;

c = css;
ce = cess;
n = nss;
lambda = lambdass;
y = yss;
k = kss;
sphi = sphiss;
i = iss;
rd = rdss;
rl = rlss;
rk = rkss;
rn = rnss;
w = wss;
q = qss;
nw = nwss;
s = sss;
sf = sfss;
z = 1;
pi = 1;
e = 1;

ez=0;
er=0;
ed=0;

end;

steady(solve_algo=0);

Thanks

Oops, sorry typo above: I was able to feed ss values on the RHS of these equations and the LHS gave me back the ss.