Im trying to solve for the steady states of the following system, originally with 28 variables. However, I have never been able to obtain the steady state, possibly because I’ve never been near a steady state enough. I tried substitution and reduce the system to 14 variables (now the equations are much longer) but still no luck. I was wondering if there is anything wrong with my model or is it just too complicated to be solved? Would there be any other choice for me other than simplifying the model?

Thank you very much.

**Original model:**

```
var c1 c2 d1 d2 la l1y l2y k1 k2 w lambc lamtfp lamer c1p c2p d1p d2p lap l1yp l2yp k1p k2p wp lambcp lamtfpp lamerp r g;
parameters nu phi phiprime chi epsilon eta mu omega alphaa delta pay payprime;
nu=0.5;
chi=1;
omega=1;
alphaa=2/3;
delta=0.025;
phi=0.1;
phiprime=0.1;
epsilon=1.1;
eta=1.1;
mu=0.01222;
pay=0.1;
payprime=0.1;
model;
g=1+mu*(la(-1)^omega+lap(-1)^omega);
r=c1p/d2;
c1+c2+(d1+d2)/r+k1(1)+k2(1)+payprime*(k2(1)*g(1)-(1-delta)*k2)^2+phiprime*c1p+wp*l1yp=l1y^alphaa*k1^(1-alphaa)+l1yp^alphaa*k2^(1-alphaa)+(1-delta)*(k1+k2)+phi*d2/r+w*l2y/r;
r*(c1p+c2p)+d1p+d2p+k1p(1)+k2p(1)+pay*(k1p(1)*g(1)-(1-delta)*k1p)^2+phi*d2+w*l2y=l2y^alphaa*k1p^(1-alphaa)+l2yp^alphaa*k2p^(1-alphaa)+(1-delta)*(k1p+k2p)+r*phiprime*c1p+r*wp*l1yp;
(c1^nu+c2^nu+d1^nu+d2^nu)^(-1)*c1^(nu-1)=lambc;
(c1^nu+c2^nu+d1^nu+d2^nu)^(-1)*c2^(nu-1)=lambc;
(c1^nu+c2^nu+d1^nu+d2^nu)^(-1)*d1^(nu-1)=lambc/r;
(c1^nu+c2^nu+d1^nu+d2^nu)^(-1)*d2^(nu-1)+lambc*phi/r=lambc/r+r*lamer/d2;
chi*(la^epsilon+(l1y+l2y)^epsilon)^(eta/epsilon-1)*la^(epsilon-1)=lamtfp*mu*omega*la^(omega-1);
chi*(la^epsilon+(l1y+l2y)^epsilon)^(eta/epsilon-1)*(l1y+l2y)^(epsilon-1)=lambc*alphaa*l1y^(alphaa-1)*k1^(1-alphaa);
chi*(la^epsilon+(l1y+l2y)^epsilon)^(eta/epsilon-1)*(l1y+l2y)^(epsilon-1)=lambc*w/r;
wp=alphaa*l1yp^(alphaa-1)*k2^(1-alphaa);
lambc(1)/g(1)*((1-alphaa)*l1y(1)^alphaa*k1(1)^(-alphaa)+1-delta)=lambc;
lambc(1)/g(1)*((1-alphaa)*l1yp(1)^alphaa*k2(1)^(-alphaa)+1-delta+2*(1-delta)*payprime*(k2(2)*g(2)-(1-delta)*k2(1)))=lambc*(1+2*payprime*(k2(1)*g(1)-(1-delta)*k2));
lambc(1)/g(1)*(alphaa*l1y(1)^alphaa*k1(1)^(1-alphaa)+alphaa*l1yp(1)^alphaa*k2(1)^(1-alphaa)+payprime*(k2(2)*g(2)-(1-delta)*k2(1))^2)=lamtfp;
lambc*(d1+d2-phi*d2-w*l2y)/r/r=lamer;
(c1p^nu+c2p^nu+d1p^nu+d2p^nu)^(-1)*c1p^(nu-1)+lambcp*r*phiprime=r*lambcp-r*lamerp/c1p;
(c1p^nu+c2p^nu+d1p^nu+d2p^nu)^(-1)*c2p^(nu-1)=r*lambcp;
(c1p^nu+c2p^nu+d1p^nu+d2p^nu)^(-1)*d1p^(nu-1)=lambcp;
(c1p^nu+c2p^nu+d1p^nu+d2p^nu)^(-1)*d2p^(nu-1)=lambcp;
chi*(lap^epsilon+(l1yp+l2yp)^epsilon)^(eta/epsilon-1)*lap^(epsilon-1)=lamtfpp*mu*omega*lap^(omega-1);
chi*(lap^epsilon+(l1yp+l2yp)^epsilon)^(eta/epsilon-1)*(l1yp+l2yp)^(epsilon-1)=lambcp*wp*r;
chi*(lap^epsilon+(l1yp+l2yp)^epsilon)^(eta/epsilon-1)*(l1yp+l2yp)^(epsilon-1)=lambcp*alphaa*l2yp^(alphaa-1)*k2p^(1-alphaa);
w=alphaa*l2y^(alphaa-1)*k1p^(1-alphaa);
lambcp(1)/g(1)*((1-alphaa)*l2y(1)^alphaa*k1p(1)^(-alphaa)+1-delta+2*(1-delta)*pay*(k1p(2)*g(2)-(1-delta)*k1p(1)))=lambcp*(1+2*pay*(k1p(1)*g(1)-(1-delta)*k1p));
lambcp(1)/g(1)*((1-alphaa)*l2yp(1)^alphaa*k2p(1)^(-alphaa)+1-delta)=lambcp;
lambcp(1)/g(1)*(alphaa*l2y(1)^alphaa*k1p(1)^(1-alphaa)+alphaa*l2yp(1)^alphaa*k2p(1)^(1-alphaa)+pay*(k1p(2)*g(2)-(1-delta)*k1p(1))^2)=lamtfpp;
lambcp*(phiprime*c1p+wp*l1yp-c1p-c2p)=lamerp;
end;
initval;
c1=2.88926199417712;
c2=2.88926199417712;
d1=11.3887510905584;
d2=1.68808278057854;
la=54.3805316898628;
l1y=15.6912341135446;
l2y=239.42601433478;
k1=0.0831311719057078;
k2=0.135914733485534;
w=0.0359395974529077;
lambc=0.010665194907865;
lamtfp=0.0184754516380136;
lamer=0.026156226574055;
c1p=3.38573459383048;
c2p=4.07134273372575;
d1p=4.78174847361626;
d2p=4.78174847361626;
lap=64.6167176936905;
l1yp=180.548846969043;
l2yp=1.90965361764221;
k1p=0.204270712225458;
k2p=182.076409874019;
wp=0.0514779962433012;
lambcp=0.0101672864916624;
lamtfpp=0.0278131543470114;
lamerp=0.0165379700568582;
r=2.00554404996009;
g=8.70680223167034;
end;
steady (solve_algo=2, maxit=200000);
```

**After substitution:**

```
var c1 d2 la l1y l2y k1 k2 c1p d2p lap l1yp l2yp k1p k2p;
parameters nu phi phiprime chi epsilon eta mu omega alphaa delta pay payprime;
nu=0.5;
chi=1;
omega=1;
alphaa=2/3;
delta=0.025;
phi=0.1;
phiprime=0.1;
epsilon=1.1;
eta=1.1;
mu=0.01222;
pay=0.1;
payprime=0.1;
model;
2*c1+((c1/((c1p/d2)^(1/(nu-1))))+d2)/(c1p/d2)+k1+k2+payprime*(k2*(1+mu*(la^omega+lap^omega))-(1-delta)*k2)^2+phiprime*c1p+(alphaa*l1yp^(alphaa-1)*k2^(1-alphaa))*l1yp=l1y^alphaa*k1^(1-alphaa)+l1yp^alphaa*k2^(1-alphaa)+(1-delta)*(k1+k2)+phi*d2/(c1p/d2)+(alphaa*l2y^(alphaa-1)*k1p^(1-alphaa))*l2y/(c1p/d2);
(c1p/d2)*(c1p+(d2p*(c1p/d2)^(1/(nu-1))))+2*d2p+k1p+k2p+pay*(k1p*(1+mu*(la^omega+lap^omega))-(1-delta)*k1p)^2+phi*d2+(alphaa*l2y^(alphaa-1)*k1p^(1-alphaa))*l2y=l2y^alphaa*k1p^(1-alphaa)+l2yp^alphaa*k2p^(1-alphaa)+(1-delta)*(k1p+k2p)+(c1p/d2)*phiprime*c1p+(c1p/d2)*(alphaa*l1yp^(alphaa-1)*k2^(1-alphaa))*l1yp;
d2^(nu-1)+c1^(nu-1)*phi/(c1p/d2)=c1^(nu-1)/(c1p/d2)+(c1p/d2)*(c1^(nu-1)*((c1/((c1p/d2)^(1/(nu-1))))+d2-phi*d2-(alphaa*l2y^(alphaa-1)*k1p^(1-alphaa))*l2y)/(c1p/d2)/(c1p/d2))/d2;
chi*(la^epsilon+(l1y+l2y)^epsilon)^(eta/epsilon-1)*la^(epsilon-1)=(((2*c1^nu+(c1/((c1p/d2)^(1/(nu-1))))^nu+d2^nu)^(-1)*c1^(nu-1))/(1+mu*(la^omega+lap^omega))*(alphaa*l1y^alphaa*k1^(1-alphaa)+alphaa*l1yp^alphaa*k2^(1-alphaa)+payprime*(k2*(1+mu*(la^omega+lap^omega))-(1-delta)*k2)^2))*mu*omega*la^(omega-1);
chi*(la^epsilon+(l1y+l2y)^epsilon)^(eta/epsilon-1)*(l1y+l2y)^(epsilon-1)=((2*c1^nu+(c1/((c1p/d2)^(1/(nu-1))))^nu+d2^nu)^(-1)*c1^(nu-1))*alphaa*l1y^(alphaa-1)*k1^(1-alphaa);
chi*(la^epsilon+(l1y+l2y)^epsilon)^(eta/epsilon-1)*(l1y+l2y)^(epsilon-1)=((2*c1^nu+(c1/((c1p/d2)^(1/(nu-1))))^nu+d2^nu)^(-1)*c1^(nu-1))*(alphaa*l2y^(alphaa-1)*k1p^(1-alphaa))/(c1p/d2);
((1-alphaa)*l1y^alphaa*k1^(-alphaa)+1-delta)=(1+mu*(la^omega+lap^omega));
((1-alphaa)*l1yp^alphaa*k2^(-alphaa)+1-delta+2*(1-delta)*payprime*(k2*(1+mu*(la^omega+lap^omega))-(1-delta)*k2))=(1+2*payprime*(k2*(1+mu*(la^omega+lap^omega))-(1-delta)*k2))*(1+mu*(la^omega+lap^omega));
c1p^(nu-1)+d2p^(nu-1)*(c1p/d2)*phiprime=c1p/d2*d2p^(nu-1)-c1p/d2*(d2p^(nu-1)*(phiprime*c1p+(alphaa*l1yp^(alphaa-1)*k2^(1-alphaa))*l1yp-c1p-(d2p*(c1p/d2)^(1/(nu-1)))))/c1p;
chi*(lap^epsilon+(l1yp+l2yp)^epsilon)^(eta/epsilon-1)*lap^(epsilon-1)=(((c1p^nu+(d2p*(c1p/d2)^(1/(nu-1)))^nu+2*d2p^nu)^(-1)*d2p^(nu-1))/(1+mu*(la^omega+lap^omega))*(alphaa*l2y^alphaa*k1p^(1-alphaa)+alphaa*l2yp^alphaa*k2p^(1-alphaa)+pay*(k1p*(1+mu*(la^omega+lap^omega))-(1-delta)*k1p)^2))*mu*omega*lap^(omega-1);
chi*(lap^epsilon+(l1yp+l2yp)^epsilon)^(eta/epsilon-1)*(l1yp+l2yp)^(epsilon-1)=((c1p^nu+(d2p*(c1p/d2)^(1/(nu-1)))^nu+2*d2p^nu)^(-1)*d2p^(nu-1))*(alphaa*l1yp^(alphaa-1)*k2^(1-alphaa))*c1p/d2;
chi*(lap^epsilon+(l1yp+l2yp)^epsilon)^(eta/epsilon-1)*(l1yp+l2yp)^(epsilon-1)=((c1p^nu+(d2p*(c1p/d2)^(1/(nu-1)))^nu+2*d2p^nu)^(-1)*d2p^(nu-1))*alphaa*l2yp^(alphaa-1)*k2p^(1-alphaa);
((1-alphaa)*l2y^alphaa*k1p^(-alphaa)+1-delta+2*(1-delta)*pay*(k1p*(1+mu*(la^omega+lap^omega))-(1-delta)*k1p))=(1+2*pay*(k1p*(1+mu*(la^omega+lap^omega))-(1-delta)*k1p))*(1+mu*(la^omega+lap^omega));
((1-alphaa)*l2yp^alphaa*k2p^(-alphaa)+1-delta)=(1+mu*(la^omega+lap^omega));
end;
initval;
c1=0.5243;
d2=14.6371;
la=16.7470;
l1y=1.8784;
l2y=14.8029;
k1=6.0658;
k2=20.9940;
c1p=6.3620;
d2p=4.6059;
lap=13.6277;
l1yp=17.8523;
l2yp=8.0043;
k1p=22.0770;
k2p=28.1082;
end;
steady (solve_algo=2, maxit=200000);
```