Solve() not working for Expenditure Minimization Problem in Sage

asked 2020-07-15 22:17:52 +0200

EconJohn
111 ●2 ●6 ●19

updated 2020-07-16 21:19:12 +0200

Though I've figured out how to go through a standard utility maximization problem analytically in sage Im having difficulty solving its dual problem, expenditure minimization. The code I have is the following:

x1, x2, l, p1, p2, a, b,  R= var('x1, x2, l, p1, p2, a, b,  R')
U = x1^a*x2^b
m = p1*x1+p2*x2;
L = m+ l * (R-U);
dLdx = L.diff(x1);
dLdy = L.diff(x2);
dLdl = L.diff(l);
solve([dLdx == 0, dLdy == 0, dLdl == 0], x1, x2, l)

The output I get however is:

[-b*l*x1^a*x2^(b - 1) + p2 == 0, -a*l*x1^(a - 1)*x2^b + p1 == 0, -x1^a*x2^b + R == 0]

which are only the first order conditions and not solving for $x_1$ and $x_2$.

This code is exactly like my utility maximization code here however I only moved where m and U are.

why is this the case?

Note: if I were to consider the following code:

x1, x2, l, p1, p2, a, b,  R= var('x1, x2, l, p1, p2,  R')
U = x1^1/3*x2^2/3
m = p1*x1+p2*x2;
L = m+ l * (R-U);
dLdx = L.diff(x1);
dLdy = L.diff(x2);
dLdl = L.diff(l);
solve([dLdx == 0, dLdy == 0, dLdl == 0], x1, x2, l)

I obtain the following analytical solution:

[[x1 == 1/2*18^(1/3)*R/(R*p1/p2)^(2/3), x2 == 18^(1/3)*(R*p1/p2)^(1/3), l==1/2*18^(1/3)*p1/(R*p1/p2)^(2/3)], 
[x1 == 1/4*18^(1/3)*R*(-I*sqrt(3) - 1)/(R*p1/p2)^(2/3), x2 == -1/2*18^(1/3)* 
(R*p1/p2)^(1/3)*(I*sqrt(3) + 1), l == 1/4*18^(1/3)*p1*(-I*sqrt(3) - 1)/(R*p1/p2)^(2/3)], 
[x1 == -1/4*18^(1/3)*R*(-I*sqrt(3) + 1)/(R*p1/p2)^(2/3), x2 == 1/2*18^(1/3)*(R*p1/p2)^(1/3)*(I*sqrt(3) - 1), l == -1/4*18^(1/3)*p1*(- I*sqrt(3) + 1)/(R*p1/p2)^(2/3)]]

My guess is that there is a bug with regards to the rules of multiplying undefined exponents.

edit retag flag offensive close merge delete

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2020-07-17 02:44:24 +0200

EconJohn
111 ●2 ●6 ●19

updated 2020-07-17 02:45:09 +0200

Thanks to Juanjo's answer to my question here I've managed to get a script which gives us an analytical solution for our hicksian demand equations. This is:

x1, x2, l, p1, p2, a, b,  R= var('x1, x2, l, p1, p2, a, b,  R')
U = x1^a*x2^b 
m = p1*x1+p2*x2;
L = m+ l * (R-U);
dLdx = L.diff(x1);
dLdy = L.diff(x2);
dLdl = L.diff(l);
step1=solve([dLdx == 0, dLdy == 0, dLdl == 0], p1, p2, R)
step2=solve((p1/p2).subs(step1)==p1/p2, x1)
good2step3 = solve((U/R).subs(step2).log().log_expand()==0, x2)
good1step4=solve(step2[0].subs(good2step3).canonicalize_radical(),x1)

where good1step4 is our hicksiand demand for $x_1$ and good2step3 is our hicksian demand for $x_2$.

Its interesting how you need a different procedure when your constraint is non-linear.

edit flag offensive delete link

add a comment

answered 2020-07-16 12:19:58 +0200

Emmanuel Charpentier
7057 ●6 ●46 ●137

updated 2020-07-16 17:28:31 +0200

It cannot be done "automagically" (neither by Maxima's solver (default) nor Sympy's), but "manually" solving each equation at a time, high school-style, works.

Elements of the problem :

1, x2, l, p1, p2, a, b,  R= var('x1, x2, l, p1, p2, a, b,  R')
U = x1^a*x2^b
m = p1*x1+p2*x2
L = m+ l * (R-U)
dLdx = L.diff(x1)
dLdy = L.diff(x2)
dLdl = L.diff(l)

Equations to solve :

E1=dLdx==0
E2=dLdy==0
E3=dLdl==0

Solve one-by-one, thus eliminating all unknowns (with a little manual help):

S1=((E1-p1).log().log_expand()-I*pi).solve(x1)
S2=((E2.subs(S1)-p2).log().log_expand()-I*pi).solve(x2)
S3=((E3.subs(S1).subs(S2)-R).log().log_expand()-I*pi).solve(l)

Substitute each unknown in partial solutions (and clean up a bit) :

Sol={S3[0].lhs():S3[0].rhs().log().log_expand().factor().exp().log_expand()}
Sol.update({S2[0].lhs():S2[0].subs(Sol).rhs().log().log_expand().factor().exp()})
Sol.update({S1[0].lhs():S1[0].subs(Sol).rhs().log().log_expand().factor().exp()})

So we get :

sage: Sol
{l: e^(-(a*log(R) + b*log(R) + a*log(a) + b*log(b) - a*log(p1) - b*log(p2) - log(R))/(a + b)),
 x2: e^(-(a*log(a) - a*log(b) - a*log(p1) + a*log(p2) - log(R))/(a + b)),
 x1: e^((b*log(a) - b*log(b) - b*log(p1) + b*log(p2) + log(R))/(a + b))}

If you happen to have access to Mathematica, you can check your solution against it :

MSol=mathematica.Solve([E1,E2,E3],[x1,x2,l])
sage: MSol
{{x1 -> E^((b*Log[a] - b*Log[b] - b*Log[p1] + b*Log[p2] + Log[R])/(a + b)), 
  x2 -> E^((-(a*Log[a]) + a*Log[b] + a*Log[p1] - a*Log[p2] + Log[R])/
     (a + b)), l -> E^((-(a*Log[a]) - b*Log[b] + a*Log[p1] + b*Log[p2] + 
      Log[R] - a*Log[R] - b*Log[R])/(a + b))}}

Which can be translated as :

MST={u[1]:u[2].sage() for u in MSol[1]}
sage: MST
{x1: e^((b*log(a) - b*log(b) - b*log(p1) + b*log(p2) + log(R))/(a + b)),
 x2: e^(-(a*log(a) - a*log(b) - a*log(p1) + a*log(p2) - log(R))/(a + b)),
 l: e^(-(a*log(R) + b*log(R) + a*log(a) + b*log(b) - a*log(p1) - b*log(p2) - log(R))/(a + b))}

Same difference...

HTH,

edit flag offensive delete link

Comments

I'm having a bit of difficulty reading the unformatted code in the middle there can you just type set it?

EconJohn ( 2020-07-16 16:32:30 +0200 )edit

My bad (a forgotten newline...). The markup used by ask.sagemath.org is not the most streamlined on the market. I'd rather use Github-flavored markdown (for which pandoc or org-mode offer convenient solutions...).

Emmanuel Charpentier ( 2020-07-16 17:34:32 +0200 )edit

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Solve() not working for Expenditure Minimization Problem in Sage

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Solve() not working for Expenditure Minimization Problem in Sage edit

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Solve() not working for Expenditure Minimization Problem in Sage