Solve cant find the answer, although I can

Hi, I have the following standard economics systema of equations to solve (if there are any economists, its a standard neoclassical model). I have to solve the following equations for a,c,k:

c + k + k*(DELTA - 1) - a*k^ALFA; == 0 (1);
1/c^SIG - (BETA*(ALFA*a*k^(ALFA - 1) - DELTA + 1))/c^SIG== 0 (2);
log(a) - RHO*log(a)== 0 (3);

From (3), a=1, from (2) k =( alpha beta / (beta (delta - 1) + 1))^(1/1-alpha) and finally from (1) we recover c.

But when using solve:

[as,cs,ks] = solve(fs(1)==0,fs(2)==0,fs(3)==0,a,c,k);

I get the error: Warning: Cannot find explicit solution. Why cant MatLab solve this system? Am I lacking some command?

THanks.

Answer 1

[John D'Errico]

Writing the equations that you have in a readable form...

(1)  c + k + k*(DELTA - 1) - a*k^ALFA == 0

(2)  1/c^SIG - (BETA*(ALFA*a*k^(ALFA - 1) - DELTA + 1))/c^SIG== 0 

(3)  log(a) - RHO*log(a)== 0

Starting with (3), I'll concede that a = 1, as long as RHO is not 1.

Then you conclude that k is computable from (2), but there is also the possibility that c=inf.

Perhaps you are ignoring some possibilities that solve was unwilling to exclude.

[Marcelo Sena]

Can I tell MatLab to ignore such cases?

[John D'Errico]

Good question. The answer is yes, but things are never as easy as you might wish.

syms a c k DELTA ALFA SIG BETA RHO
E1 = c + k + k*(DELTA - 1) - a*k^ALFA == 0;
E2 = 1/c^SIG - (BETA*(ALFA*a*k^(ALFA - 1) - DELTA + 1))/c^SIG== 0;
E3 = log(a) - RHO*log(a)== 0;

For example...

a = solve(E3,a)
Warning: The solutions are valid under the following conditions: in(exp((pi*l*(-2*i))/(ALFA - 1))*((DELTA + 1/BETA - 1)/ALFA)^(1/(ALFA - 1)),
'real'). To include parameters and conditions in the solution, specify the 'ReturnConditions' option. 
> In solve>warnIfParams at 518
  In solve at 361 
a =
1

Ok, so it looks like I can tell MATLAB not to worry so much about a.

syms a
assume(a,'real')
a = solve(E3,a)
a =
1

Next, I looked at (2), to see if I can solve for k.

solve(E2,k)
Warning: The solutions are parameterized by the symbols: l. To include parameters and conditions in the solution, specify the 'ReturnConditions'
option. 
> In solve>warnIfParams at 511
  In solve at 361 
Warning: The solutions are valid under the following conditions: 0 <= real(ALFA - 1)*(pi - imag(log((DELTA + 1/BETA - 1)/(ALFA*a))/(ALFA - 1)) +
2*pi*l*real(1/(ALFA - 1))) & 0 < real(ALFA - 1)*(pi + imag(log((DELTA + 1/BETA - 1)/(ALFA*a))/(ALFA - 1)) - 2*pi*l*real(1/(ALFA - 1))) & 1 <
real(ALFA) & in(l, 'integer') | real(ALFA - 1)*(pi - imag(log((DELTA + 1/BETA - 1)/(ALFA*a))/(ALFA - 1)) + 2*pi*l*real(1/(ALFA - 1))) <= 0 &
real(ALFA - 1)*(pi + imag(log((DELTA + 1/BETA - 1)/(ALFA*a))/(ALFA - 1)) - 2*pi*l*real(1/(ALFA - 1))) < 0 & real(ALFA) < 1 & in(l, 'integer'). To
include parameters and conditions in the solution, specify the 'ReturnConditions' option. 
> In solve>warnIfParams at 518
  In solve at 361 
ans =
exp(-(pi*l*2*i)/(ALFA - 1))*((DELTA + 1/BETA - 1)/(ALFA*a))^(1/(ALFA - 1))

Oops. The problem here is not only will I need to constrain k to be real, but I'll need to know something about the other parameters, like ALFA, BETA, DELTA, RHO, SIG, c. If they can be anything, including even complex values, solve gets upset.

So the problem is when you solve for k there, syms realizes that there are some nasty cases to worry about, possibilities you chose to ignore. So I made a few guesses here. Really, I'm being lazy, just applying wholesale constraints on the parameters.

assume(c,'real')
assume(k,'real')

assume(ALFA,'real')
assume(BETA,'real')
assume(DELTA,'real')
assume(RHO,'real')
assume(SIG,'real')

assumeAlso(0 < RHO)
assumeAlso(1 < ALFA)
assumeAlso(0 < DELTA)
assumeAlso(0 < BETA)
assumeAlso(0 < SIG)

assumeAlso(0 < c)

Lets see how this gets me.

solve(subs(E2,'a',1),k)
Warning: The solutions are parameterized by the symbols: l. To include parameters and conditions in the solution, specify the 'ReturnConditions'
option. 
> In solve>warnIfParams at 511
  In solve at 361 
Warning: The solutions are valid under the following conditions: in(exp((pi*l*(-2*i))/(ALFA - 1))*((DELTA + 1/BETA - 1)/ALFA)^(1/(ALFA - 1)),
'real') & pi*l < (ALFA - 1)*(pi/2 + pi/(2*(ALFA - 1))) & 1 <= ALFA/2 + l & BETA*DELTA + 1 < BETA & in(l, 'integer') | in(exp((pi*l*(-2*i))/(ALFA -
1))*((DELTA + 1/BETA - 1)/ALFA)^(1/(ALFA - 1)), 'real') & 2*l + 1 < ALFA & 1 <= ALFA + 2*l & BETA <= BETA*DELTA + 1 & in(l, 'integer'). To include
parameters and conditions in the solution, specify the 'ReturnConditions' option. 
> In solve>warnIfParams at 518
  In solve at 361 
ans =
exp(-(pi*l*2*i)/(ALFA - 1))*((DELTA + 1/BETA - 1)/ALFA)^(1/(ALFA - 1))

So look more closely as what we are being told. I think the problem is there are constraints on parameters like:

DELTA + 1/BETA - 1 > 0

If not, then solve forsees a problem. So, we can use assumeAlso again.

assumeAlso(DELTA + 1/BETA - 1 > 0)
solve(subs(E2,'a',1),k)
Warning: The solutions are valid under the following conditions: BETA <= BETA*DELTA + 1. To include parameters and conditions in the solution,
specify the 'ReturnConditions' option. 
> In solve>warnIfParams at 518
  In solve at 361 
ans =
((DELTA + 1/BETA - 1)/ALFA)^(1/(ALFA - 1))

So, there is one more piece of useful information to provide:

assumeAlso(BETA <= BETA*DELTA + 1)
k = solve(subs(E2,'a',1),k)
k =
((DELTA + 1/BETA - 1)/ALFA)^(1/(ALFA - 1))

Sigh.