load(draw)$ 
set_draw_defaults(fill_color=gray40,
draw_realpart=false,font_size=10)$
fpprintprec:5$ ratprint:false$ kill(all)$

# Demand Theory: Constraints 
               and Optimization

This chapter introduces the budget constraint and 
combines the budget constraint with the preferences developed in the 
preceding chapter to provide a model of rational consumer
choice. 

## Budgets and Budget Constraints

## Indifference Curves and Choice

## Lagrangian Multipliers

Langrangian multipliers offer an alternative to graphical analysis 
of optimization. To apply this method,  augment the objective
function by adding a term (or terms) containing the constraint(s), 
stated as an implicit function (implicit functions), multiplied by the
term(s) M (M1, M2 ...), the Lagrange Multiplier(s).

The next cell illustrates the approach using one constraint and 
a budget line as the constraint. Now utility is a function of four
goods: x, y, z, and w. The function L, which consists of the 
utility function and the constraint, which equals zero, achieves 
its extreme value at the same set of values at which u(x,y,z,w)
achieves its extreme. 

The depends( ) command tells Maxima that u is a function of 
x, y, z, and w.

kill(all)$   depends(u,[x,y,z,w])$
L : u(x, y, z, w) + M*(m - px*x - py*y - pz*z - pw*w);

TeX:
The five partial derivatives used to define the first-order condition 
for optimization appear below. Note that the partial derivaive Lm
returns the budget line, i. e., the constraint.

[Lx,Ly,Lz,Lw,LM]:[diff(L,x),diff(L,y),diff(L,z),
                           diff(L,w),LM:diff(L,M)]$   
transpose( matrix(
 ["Lx","Ly =","Lz =","Lw =","Lm ="],[Lx,Ly,Lz,Lw,LM]));

The cell below shows the expressions for the CES utility function of two 
goods.

kill(all)$  
matrix(
  [u(x,y,a,b):=(a*x^-b+(1-a)*y^-b)^-(1/b),
   L :  u(x, y, a, b) + M*(m - px*x - py*y)]);
 transpose(matrix(
   [Lx(x, y, a, b, px, py, m, M) := ''(diff(L, x)),
    Ly(x, y, a, b, px, py, m, M) := ''(diff(L, y)),
    LM(x, y, a, b, px, py, m, M) := ''(diff(L, M)) ] ) ) ;

### Illustration of the Lagrangian Multiplier

Using mnewton to find the (approximately) optimal values, we illustrate 
the fact that M is the marginal utility of income. The cell below defines
parameters and solves the system.

[a0,b0]:[1/3,4/5]$[px0,py0,m0]:[50,12,1389.4]$
load("mnewton")$     soln0 : float(mnewton(
[Lx(x,y,a0,b0,px0,py0,m0,M),Ly(x,y,a0,b0,px0,py0,m0,M),
 LM(x,y,a0,b0,px0,py0,m0,M)],[x,y,M],[20,20,0.01]))[1];

The table below shows these values, along with the initial utility level, u = 30.
Given that M = 0.0216, adding $1 to m should increase u to approximately
u = 30.0216.

[x0,y0,M0]:[rhs(soln0[1]),rhs(soln0[2]),rhs(soln0[3])]$
matrix(["x0","y0","M0","u0"],[x0,y0,M0,u0:u(x0,y0,a0,b0)]);

The cell below confirms this result.

soln1:float(mnewton([ Lx(x,y,a0,b0,px0,py0,m0+1,M),  
Ly(x,y,a0,b0,px0,py0,m0+1,M),LM(x,y,a0,b0,px0,py0,
  m0+1,M)],[x,y,M],[20,20,0.01]))[1]$  [x1,y1]:
 [rhs(soln1[1]), rhs(soln1[2])]$  u1:u(x1,y1,a0,b0)$
u_change: u(x1,y1,a0,b0) -u(x0, y0,a0,b0)$  matrix(
 ["x1","y1","u1","u__change"],[x1,y1,u1,u_change]);

A $1 income decrease reduces u by 0.0216 units.

soln2:float(mnewton([ Lx(x,y,a0,b0,px0,py0,m0+1,M),  
Ly(x,y,a0,b0,px0,py0,m0+1,M),LM(x,y,a0,b0,px0,py0,
  m0-1,M)],[x,y,M],[20,20,0.01]))[1]$  [x2,y2]:
 [rhs(soln2[1]), rhs(soln2[2])]$  u2:u(x2,y2,a0,b0)$
u_change: u(x2,y2,a0,b0) -u(x0, y0,a0,b0)$  matrix(
 ["x2","y2","u2","u__change"],[x2,y2,u1,u_change]);

More generally, a Lagrangian multiplier states the change in the value 
of the objective function per one-unit change in the constraint. 

### Multiple Constraints

With Lagrangian approach generalizes, the analysis in that it can 
to any number of goods. Also, additional constraints can be added. 

L : u(x,y,a,b)+M1*(m-px*x-py*y)+M2*(T -tx*x-ty*y);

The Lagrangian equation now contains four variables: x, y, and the two 
multipliers. The first-order conditions are defined by setting the four
partial derivatives equal to zero.

transpose(matrix(["Lx =","Ly =","LM1 =", "LM2 ="],
 [diff(L,x),diff(L,y),diff(L,M1),diff(L,M2)]));

### Expenditure Minimization

The constrained maximization problem stated thus far (maximize
utility subject to a budget constraint) can be rephrased in terms
of its dual: minimize the expenditures required to provide
a specified utility level. The Lagrangian equation W below
states this problem for the CES utility function. The Lagrangian
multiplier is now labeled mu.

W:m-px*x-py*y+%mu*(u0-u(x,y,a,b));

The partial derivatives required to determine the expenditures-
minimizing combination of x and y appear below. 

Wx:diff(W,x)$ Wy:diff(W,y)$ Wmu:diff(W,%mu)$
transpose(matrix(
  ["Wx =","Wy =","Wmu ="],[Wx,Wy,Wmu]));

The multiplier from the dual problem is the change in expenditure 
implied by a one-unit change in the target utility level, so mu = 1/M.

Setting Wx and Wy equal to zero and solving allows an expression
for the optimal value of y/x, named yxratio below. The budget 
constraint can be written as m = px*xopt + py*yxratio*xopt, where
xopt is the optimal quantity of x, to be determined. Then yxratio
implies yopt, the optimal quantity of y.

matrix([yxratio:((px/py)*(1-a)/a)^(1/(1+b)),
            xopt:m/(px+py*yxratio),
            yopt:xopt*yxratio] );

Substituting these values into the consumer's
utility function produces an expression for utility
given the constraints and optimizing behavior, the 
resulting expression, the indirect utility function.

Utility is now stated in terms of the original utility function's 
parameter, along with the prices and income level.
The x and y values do not appear directly (hence, "indirect.")

v(px, py, m, a, b) :=''(u(xopt, yopt, a, b))$
v0: float(v(px0, py0, m0, a0, b0))$
v1:float(v(px0,py0,m0+1,a0,b0))$  chg:v1-v0$
transpose(matrix(
  ["v =","Initial Utility, v =","New Utility, v =",
  "change ="], [u(xopt,yopt,a,b),v0, v1, chg]));

We use this formulation to examine the effects of two taxes,
one on the consumption of good x and the other a lump-sum tax
that is unrelated to the consumer's consumption pattern.

x0:float(subst([m=m0,px=px0,py=py0,a=a0,b=b0],(xopt)))$
/* Replace px0 with px0*1.5 -- impose a 50% ad valorem tax. */
x1:float(subst([m=m0,px=px0*1.5,py=py0,a=a0,b=b0],(xopt)))$           
taxRev : %*0.5*px0$      
x2: float(subst([m =m0-taxRev,px=px0,py=py0,a=a0,b=b0],xopt))$
matrix(["x /no tax","x /ad valorem tax","x /lump sum tax",
"Tax revenue"], [x0, x1, x2, taxRev]);

At the initial values of the system's parameters, xopt = 15.616. 
A 50% tax on x reduces x to 11.221 units and generates 280.53 
in revenue per time period. 
 
A per-period lump-sum tax of this amount results in xopt = 12.463.
The lump-sum tax reduces the consumer's utility by less than the 
ad valorem tax.

v0: float(v(px0*1.5, py0, m0, a0, b0))$
v1: float(v(px0, py0, m0-taxRev, a0, b0))$
matrix(["Utility/ad valorem tax", 
"Utility/lump sum tax"], [v0, v1]);

## Second-Order Conditions

This section develops the first-order and second-order conditions and
applies the analysis to the utility maximization model.

### Functions of One Variable

The expression in the next cell is a quadratic function of x.
It illustrates the nature of the second-order for a function
of one variable. (Some functions have multiple local extreme
values.)

kill(all)$ f(x,a,b,c):=a+b*x+c*x^2$
dfdx:diff(f(x,a,b,c),x)$   
xsoln:rhs((solve(dfdx, x ))[1])$
d2fdx2:diff(f(x,a,b,c),x,2)$   
matrix(["f(x)","df/dx","x,f=0","d2f/dx2"],
       [f(x,a,b,c),dfdx,xsoln,d2fdx2]);

The next cell specifies parameters for two parabolas, both with with 
extreme values when x > 0. For each, we determine the value of x 
that yields an extreme value, and the extreme value.

[a0,b0,c0]:[500,5,-0.5]$  [a1,b1,c1]:[150,-1.5,1.0]$
xsoln0:subst([ a=a0,b=b0,c=c0],xsoln)$    xsoln1:
  subst([a=a1,b=b1,c=c1],xsoln)$  xmax:max(xsoln0,xsoln1)$
ymax:max(f(xsoln0,a0,b0,c0),f(xsoln1,a1,b1,c1))$
ymin:min(f(xsoln0,a0,b0,c0),f(xsoln1,a1,b1,c1))$
matrix(
  ["","a","b","c","x, dy/dx = 0","Extreme Value of y"],
  ["","--","--","--","----","-----"],
  ["Initial",a0,b0,c0,xsoln0,f(xsoln0,a0,b0,c0)],
  ["New",a1,b1,c1,xsoln1,f(xsoln1,a1,b1,c1)] );

TeX: The graph confirms that the extreme (maximum or minimum)
values of the function occur at the x values calculated above. It also 
shows that when c < 0 the extreme value is a maximum, and when
c > 0 the extreme value is a minimum.

wxdraw2d(user_preamble="set key bottom",xlabel="x",
 ylabel="f(x)",yrange=[.5*ymin,1.1*ymax], key=concat(
"Maximum: c=",string(c0)), explicit(if x < xsoln0 
  then f(xsoln0,a0,b0,c0)     else 0,x,-10,xsoln0),
 line_type=dots, key=concat("Minimum:c =",string(c1)),
 explicit(if x < xsoln1 then f(xsoln1,a1,b1,c1) 
   else 0,x,-10,xsoln1), line_type = solid, key =
  string(f(x,a0,b0, c0)), line_width=2,explicit(
  f(x,a0,b0,c0),x,-10,5*xmax),color=red,key = 
  string(f(x,a1,b1,c1)),line_width=3, 
 explicit(f(x,a1,b1,c1),x,-10, 5*xmax),
 dimensions=[640,480])$

The illustration below shows a simple function that has
two extremes over the range considered. At the maximum
x, dz/dx = 0 and the the second derivative is negative
(-1, to be precise, see d2zdx2). Similar reasoning applies
to the x value that generates the minimum z value.

We use this simple example because the initial function
and its derivatives are in the same units (more precisely,
they are unit-free). For most economic examples, placing
the initial function and derivatives in the same graph
is a bad idea.

kill(Z)$  Z:sin(x)$ dZdx: diff(Z,x)$ d2Zdx2:diff(Z,x,2)$
print("Z = ", Z, ", Zx = dZ/dx =", dZdx,", and Zxx =", d2Zdx2)$
wxdraw2d(xaxis=true,user_preamble="set key left",
 xlabel="x",ylabel="Z and dZ / dZ", 
 key="Z", explicit(Z,x,-%pi,%pi),
 color=red, key="dZdx", explicit(dZdx,x,-%pi,%pi),
 line_type=dots,points_joined=true, key="",
 color=black, points([[-%pi/2,-1],[-%pi/2,1]  ]),
 points([[%pi/2,-1],[%pi/2,1]  ])
)$

### Functions of More than One Variable

Extend the analysis to a function u = u(x, y, z, w). To achieve an 
extreme value, it must small changes  in any of the variables must 
not affect u's value. The text develops the nature of the second-
order conditions, that determine whether an extreme value is a 
local minimum or a local maximum.

These condtions involve a Hessian matrix like the one below.
A Hessian matrix,is a a square matrix that consists of all 
second-order derivatives of the function being examined.

H1: matrix([u[xx],u[xy],u[xz]],
    [u[yx],u[yy],u[yz]],[u[zx],u[zy],u[zz]]);

Before addressing the second-order conditions, we review 
two important aspects of a square matrix: the determinant
and the principal minors.

The Determinant

The determinant of a square matrix is a value that results from 
a process of cross-multiplication of the matrix's terms. The
first output line below lists a 2 x 2 matrix and its determinant.
The second output line shows the determinant of H1.

[matrix([a,b],[c,d]),determinant(matrix([a,b],[c,d]))]; 
determinant(H1);

Principal Minors

To determine whether second-order conditions are met, examine
a Hessian matrix's principal minors, which are sub-determinants 
of the initial matrix . An n x n  matrix  has n principal minors. 
The principal minors for the 2 x 2 matrix above appear below. 

kill(A)$ A:matrix([a,b],[c,d]); matrix(
  ["1st principal minor"," 2nd principal minor "],
  [determinant(minor(A,2,2)), determinant(A)]);

The Hessian matrix H1 has three principal minors.

The second-order condition for an n x n matrix depends 
on whether a maximum value or a minimum value of the function
is indicated, as indicated in the text. 

We apply this rule to a two equations that have two independent
variables, so that graphs illustrate the conditions. Then
we consider an example in which a consumer chooses levels of three
goods subject to an income constraint.

The first example has neither a maximum nor a minimum
value. Rather, the (single) extreme value of this function is at a 
minimum  when approached from one direction and a maximum
when approached from another direction.  

f(x, y) := (x -1)^2 -(y-2)^2$   fx : diff(f(x, y), x)$ 
fy:diff(f(x,y),y)$  fxx:diff(f(x,y),x,2)$  fxy:diff(fx,y)$
fyx : diff(fy, x)$    fyy:diff(f(x, y), y, 2)$    
matrix(["f(x,y)", "fx", "fy"], [f(x, y), fx, fy]);  
matrix(["fxx","fxy","fyx","fyy"],[ fxx,fxy,fyx,fyy]);

The graph illustrates the nature of the saddle point.

wxdraw3d(view = [60, 300], xlabel ="x",
explicit(f(x, y), x, -10, 10, y, -10, 10) ,
dimensions = [640,480])$

The cell below animates the view of this saddle point.

with_slider_draw3d(iList, makelist(15*i,i,0,8),
view = [60, iList], xlabel ="x",
explicit(f(x, y), x, -10, 10, y, -10, 10),
dimensions=[640,480] )$

TeX: The behavior described above suggests that this function should
fail the second-order conditions that are listed in the table above.
It does: that |H21| = 2 violates the condition required for 
having found a maximum value, and that |H22|  = -4 violates the 
condition required for having found a minimum value. 

[fx, fy], x=1, y=3$ H2 : matrix([fxx, fxy], [fyx, fyy]);
H21: determinant(H2[1,1])$  H22: determinant(H2)$
matrix(["|H21|", "|H22|"], [H21, H22]);

TeX: Next consider an expression that results in a minimum
value for the function at its extreme point.

g(x, y):=x^2 +y^2-2*x-6*y+14$  gx:diff(g(x,y),x)$    
gy:diff(g(x,y), y)$  gxx: diff(g(x,y),x,2)$   gxy:diff(gx,y)$
gyx : diff(gy, x)$    gyy:diff(g(x,y),y,2)$  
matrix(["g(x,y)","gx","gy"], [g(x,y),gx,gy]);
matrix(["gxx","gxy","gyx","gyy"], [gxx,gxy,gyx,gyy]);

Both gxx and gyy are positive, so both gx and gy are negative for 
small values of x and y and positive for larger values. Approaching (1, 3)
from either direction causes g(x, y to decrease. The graph below
confirms that the figure is bowl-shaped and has a minimum as its
extreme value.

wxdraw3d(view=[70,310],xlabel="x", enhanced3d=true,
 explicit(g(x,y),x,0,2,y,2,4)  )$

For comparison, look at the same graph with enhanced3d=false.

wxdraw3d(view=[70,310],xlabel="x", enhanced3d=false,
 explicit(g(x,y),x,0,2,y,2,4)  )$

TeX: Finally, the Hessian H3 satisfies the conditions for a minimum value
to have been identified. The first-order conditions are met: $gx = gy = 0$.
Also, the second order conditions are met: $|H31| > 0$ and $|H32 > 0|$.

subst([x=1,y=3],[gx, gy]);  H3:matrix([gxx,gxy], [gyx,gyy]);
H31: determinant(H3[1,1])$  H32: determinant(H3)$
matrix(["|H31|", "|H32|"], [H31, H32]);

Application to Optimization

We illustrate a constrained optimization problem using a 
Cobb-Douglas utility function of three goods and a 
corresponding budget constraint.  This cell shows the original
functions, the implied Lagrangian equation, and the first partial
derivatives of the Lagrangian equation.

G:40-(2*x+3*y+4*z)$   u:x^.2*y^.5*z^.3$    L:u+M*(G)$
[Lx,Ly,Lz,LM]:[diff(L,x),diff(L,y),diff(L,z),diff(L,M)]$
[Lxx,Lyy,Lzz]:[diff(L,x,2),diff(L,y,2), diff(L,z,2)]$
[Lxy, Lxz, Lyz] :[diff(Lx, y),diff(Lx, z), diff(Lz,y)]$
[Gx, Gy, Gz] : [diff(G,x), diff(G, y), diff(G, z)]$
transpose(matrix(["Utility Function","---",
"Budget Constraint","---","Lagrangian","---",
"Lx","Ly","Lz"],[u,"---",G,"---", L,"---",Lx,Ly,Lz]));

The "bordered" Hessian matrix below contains the 
second-order partial derivatives of the utility function.

The last row and the last column, which form the border, 
are the derivatives of the constraint. 

HB: matrix([Lxx,Lxy,Lxz,Gx],[Lxy,Lyy,Lyz,Gy],
   [Lxz,Lyz,Lzz,Gz],[Gx,Gy,Gz,0] );

Next, mnewton( ) solves for x, y, and z. The resulting utility level 
is reported. 

load(mnewton)$  soln:mnewton([Lx,Ly,Lz,LM],[x,y,z,M],
  [.1,.1,.1,.1])$   [x0,y0,z0,u0]:[rhs(soln[1][1]),
  rhs(soln[1][2]), rhs(soln[1][3]),subst(soln[1],u)]$
  matrix(["x0","y0","z0","u0"],[x0,y0,z0,u0]);

Confirm that the consumer's expenditures are as follows: 20% 
on x, 50% on y, and 30% on z. 

/*Recall that m=40, px=2, py=3, and pz=4*/
[2*x0/40, 3*y0/40, 4*z0/40];

Evaluate the terms in the bordered Hessian, using x0, y0, and z0.

HB0: subst(soln[1],HB);

Like extreme values in the absence of a constraint, the second-order
conditions relate to the principal minors of the bordered Hessian 
matrix.  With the bordered Hessian, however, only border-preserving 
principal minors are considered. 

First, consider the order-preserving principal minor with r = 1
($n - r = 2,
because $n = 3$ in this illustration). As always, the principal minor for 
$r = 1$ is a four-element matrix with a zero in the lower right-hand corner
and identical terms as the off-diagonal elements. The determinant must 
be negative, as it is in this case: 0 - (-4)^2 = -16, as the cell below
shows.

Matrix HB0 is repeated here, along with the principal minor 
described above. 

[HB0, matrix([HB0[3][3],HB0[3][4]],[HB0[4][3],HB0[4][4]]) ];

The matrix( ) command is being used in two different ways. The 
outermost command is used to build a table. The interior commands
are used analytically. 

matrix(['determinant(matrix([HB0[3][3],HB0[3][4]],
  [HB0[4][3],HB0[4][4]])),"=",determinant(matrix(
  [HB0[3][3], HB0[3][4]],[HB0[4][3],HB0[4][4]  ]))]);

The second order condition to establish a maximum is this: Given that 
the first-order conditions are satisfied, the border-preserving 
principal minors of order r must have the sign (-1)^r, for 
r = 2, ... , n. That is, the border-preserving principal minors must 
alternate in sign, beginning with a positive second-order minor.

In this example, n = 3, so we consider r = 2 (n - r = 1) and r = 3 
(n - r = 0), the latter being the determinant of the original matrix. 
Thus, we delete the two determinants below, with the attached number 
indicating the value of r. The first term is positive, 
and the third term is negative, as they should be.

[minor(HB0,1,1), 
     [HB2:determinant(minor(HB0,1,1)),
      HB3 : determinant(HB0)]];

If computation alone is required, then the hessian( ) command provides the 
information directly. 

HB:hessian(L,[x,y,z,M]); 
HB0:subst(soln,HB);

The values satisfy the first-order conditions yield maximum
utility.

We could have stated the problem as an expenditures minimization 
problem. The resulting values of x, y, and z would have been the 
same, but the second-order condition would have been different. The 
second order condition to establish a minimum is this: Given 
that the first-order conditions are satisfied, the border-preserving 
principal minors of all orders must be negative.