Problem solving for steady state

[diatrochoi]

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:

Code: Select all

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:

Code: Select all

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);

[jpfeifer]

You should try to analytically compute the steady state. This does not simply involve plugging one equation into another one. You need to find a handle on how to sequentially compute the endogenous variables.

[diatrochoi]

Thanks for the reply. I have tried doing so by writing a separate steady state solving file but so far no luck yet. I will try one more time. However, if the model is too complicated for each variable to be solved explicitly, is it correct to say that my only choice then would be to simplify the model's assumption and reduce the model to an easier one?

[jpfeifer]

I have not yet seen a model where the steady state computation cannot be broken down into one or two equations that can be solved numerically in the end. Take a look at the NK_baseline.mod in the Dynare example folder for one case dealing with a rather complicated model.

[diatrochoi]

I have tried but am still able to reduce the model to the point where a variable can be solved explicitly in terms of just the parameters as in the example provided by dynare. If I cannot express a variable in terms of just the parameters, does it mean there is no way to solve the model analytically? I have attached a pdf writeup of a model. If it is possible, could you give me some hints on how I could solve this analytically? Thank you very much and I really appreciate your help.

[jpfeifer]

Have you tried calibrating labor and then solving the steady state. That makes the life often easier.

[diatrochoi]

Do you mean setting all the labors as parameters over a range, then randomly select their values within that range, then solve for the other variables and finally check the square of residuals in the equations that haven't been used? That is what I am doing as of now. However, it is not working yet and I feel like I'm just throwing darts randomly. Would there be another way? Thanks

[jpfeifer]

No, I mean the often used approach of setting labor in steady state to 1/3 of the available time and then later on finding the labor disutility parameter that makes this value for labor the steady state value. In this case you can often rather easily solve for the other variables

[diatrochoi]

In my model the amount of labor that agent can supply is not fixed. There are also 3 types of labor, i.e. R&D labor, production labor dedicated to product 1 and production labor dedicated to product 2. So do you mean that I should set the total quantity of labor fixed and then divide that by 3 to arrive at the calibration value for each type of labor?

[jpfeifer]

Please read my statement carefully: it said you could try to fix labor to a value in steady state and then find the parameters that are compatible with this setting. This is often easier than fixing a parameter and then finding labor that is consistent with it. That often simplifies things, because with labor given in a first step, you can often directly compute capital from the Euler equation.

[paidjo]

Dr. Pfeifer,

I have a question related to this problem. So in my code, I fix steady-state labor l_ss = 1/3 then solve for bunch of steady state value of variables using it. I know that I need to calibrate the disutility of labor parameter (for example, CHIW). But when I do this, i.e., setting CHIW = lambda_ss*w_ss*(1-l_ss) from the FOC of labor, the code won't work and dynare will update my steady-state value of l, which is then different than what I initially set. The code works when I set both l_ss = 1/3 and CHIW = 0.5, for example. Can you tell me what I did wrong? Thank you.

[jpfeifer]

Usually this implies you did make a mistake in your computations. If the parameter is correctly updated, the steady state should have the correct value for labor.