# Optimization, Part B

load(draw)$  set_draw_defaults(line_width=2, grid = [2,2], point_type = filled_circle,
          head_type = 'nofilled, head_angle = 20, head_length = 0.5,
             background_color = light_gray,  draw_realpart=false)$     
 fpprintprec : 5$  ratprint : false$

load ("Econ1.mac");

## Preface

Dowling12B.wxmx uses Maxima to review the analysis of optimization of functions of 
economic variables which are constrained by equalities, following the second part of Ch. 12 
of Introduction to Mathematical Economics (3rd ed),  by Edward T. Dowling, 
(Schaum's Outline Series), McGraw-Hill, 2012. This text is a bargain, with many complete 
problems worked out in detail. You should compare Dowling's solutions, worked out 
"by hand", with what we do using Maxima here.

A code file Econ1.mac as available in the same section (of Economic Analysis with Maxima),
which defines many Maxima functions used in this worksheet.
Use load ("Econ1.mac");

This worksheet is one of a number of wxMaxima files available in the section
    Economic Analysis with Maxima
on my CSULB webpage.

Edwin L. (Ted) Woollett
https://home.csulb.edu/~woollett/
June 13, 2022

Part A includes unconstrained optimization examples. In Part B the use of the bordered 
Hessian matrix for second order conditions in optimization with equality constraints is 
discussed, and also the use of eigenvalues determination of sign definiteness. Also 
discussed in some detail in Part B is the subject of Input-Output Analysis in an open 
economy.

We have slightly changed some of the symbols used by Dowling in particular problems.

## References

Notes on a High School course on optimization methods with a large variety of
approaches to the problem.

https://web.stanford.edu/group/sisl/k12/optimization/#!index.md
-----------------------------------------------------
Course notes, problem sets, and Exams for Math 2640, Introduction to Optimization, 
Leeds Univ., Chris Jones, 2005

http://www1.maths.leeds.ac.uk/~cajones/math2640/MATH2640.html
--------------------------------------------------
Fundamental Methods of Mathematical Economics, Alpha C. Chiang and 
Keven Wainwright, 4th ed., 2005, McGraw-Hill
----------------------------------
Wainwright's course notes at British Columbia Institute of Technology, Burnaby, 
British Columbia, Canada
http://faculty.bcitbusiness.ca/kevinw/chiang/ChapterLectureNotes.htm
----------------------------------------------
Wainwright's 2007 Econ 331 course at Simon Fraser Univ.
http://www.sfu.ca/~wainwrig/Econ331/331.htm
-------------------------------------------------------
Martin J. Gander, Math. Dept.,  Univ. of Geneva,  Ch. 4, Optimization
https://www.unige.ch/~gander/teaching/polycopie.pdf
https://www.unige.ch/~gander/
------------------------------------------------------
Math 2070, Univ. of Sydney Optimization notes

https://www.maths.usyd.edu.au/u/UG/IM/MATH2070/r/NLoptWC.pdf
----------------------------------------------------
Wolfram summary of Mathematics methods for optimization.
https://reference.wolfram.com/language/tutorial/ConstrainedOptimizationIntroduction.html

## Optimization Constrained by Equalities [12.5]

### optimum (f, varL, h),  optimum (f, varL, [h1, h2] ), etc.

For one equality constraint h = 0, use optimum (f, varL, h).
If f depends on the economic variables x and y, varL = [x, y],
If f depends on the economic variables x, y, and z, varL = [x, y, z], etc.
For two equality constraints h1 = 0 and h2 = 0, use optimum (f, varL, [h1, h2]).
And so on for more equality constraints.

### Bordered Hessian: Bhessian (L, gList, varList), BHtest (BH, m, n)

Assume the objective function f is a function of n variables (x1, x2, ...,xn), and we have 
m constraints of  form g(x1, x2, ...,xn) = c, where c is a constant.

The bordered Hessian matrix BH can be calculated with the syntax

       BH : Bhessian (L, gList, varList)
in which L is the Lagrangian, glist is (in general) a list of expressions like [g1, g2] related
to the equality constraints in force, and varList is a list of the search variables [x,y,...].

For only one equality constraint, such as g = c, where g is a Maxima expression depending
on search variables (x,y), c is a constant, and the objective function f depends also on (x,y),
the Lagrangian function is
    L : f + λ*(c - g) OR (for example)  L : f + lam*(c - g),
and  the bordered Hessian matrix is
    BH : Bhessian (L, g, [x,y])  OR   BH : Bhessian (L, [ g ], [x,y]). 

To insert the Greek letter λ at any point, press Esc, then the letter l, and then press Enter.
You can also use some ordinary name, such as lam, or lam1, for the Lagrangian multiplier.

If the optimization is subject to two equality constraints, for example g1 = c1, g2 = c2, we
need two Lagrange multipliers to define L
    L : f + λ1*(c1 - g1) + λ2*(c2 - g2),  OR (e.g.)  L : f + lam1*(c1 - g1) + lam2*(c2 - g2),
    BH : Bhessian (L, [ g1, g2 ], [x,y]). 

For more constraints and/or more search variables, just enlarge the lists.

The Maxima function Bhessian is defined in Econ1.mac, and the above syntax produces
 a   (m + n) x (m + n)  square matrix which will involve the first derivatives of the m constraint 
functions g1, g2,..., and the second derivatives of the Lagrangian function L with respect to
 (wrt) the n search variables in varList. In the most general case this bordered Hessian 
matrix will be a function of all n search variables and all m lagrange multipliers λk.

In a general case, assume n search variables x1, x2,..., xn, and m equality constraints
g1 = c1, g2 = c2,...gm = cm.

To test for a maximum or minimum, the Hessian matrix BH should first be evaluated at a 
critical point, which is a solution of the m + n equations grad(L) = 0, where grad(L) is the 
list of the first derivatives of the Lagrangian L with respect to (wrt) the m+n "variables"  
x1, x2,..., xn, λ1, λ2,...,λm.

Let BHk stand for the "kth bordered leading principal minor" of a bordered Hessian 
matrix in which the smallest value of k = m+1, and the largest k = n.

You can calculate BHk using our Maxima function BLPM,
    BHk :  BLPM (BH, m, k) 
which returns LPM (BH, m + k), with k = m + 1 thru k = n. There are (n - m) leading principal 
minors of a bordered Hessian matrix: beginning with LPM (BH, 2 m + 1) and ending with 
LPM (BH, m + n), using our function LPM, which we discussed in Dowling12A.wxmx, Sec.
 "Leading Principal Minors of a Matrix  LPM".

But it is easiest instead to follow the path:
  1. Calculate the Lagrange expression L.
  2. Calculate the bordered Hessian matrix BH : Bhessian (L, gList, varList).
  3. Evaluate BH at a candidate critical point at which the first derivatives of L wrt the
search variables and the Lagrangian multipliers are all simultaneously equal to zero.
  4  Call BHtest (BH, m, n), in which m = number of equality constraints, n = number of
search variables, and BH corresponds to the desired critical point.

BHtest (BH, m, n) returns a list of the n - m leading principal minors of the bordered Hessian 
matrix BH, begining with LPM (BH, 2*m+1) and ending with LPM (BH, m + n) = det BH, and 
if the relevant criteria set out below are met, prints out the nature of the critical point used 
to generate BH.

Let sgn(x) = sign(x). For example, sgn(-2) = -1, sgn(6) = +1.

The number of relevant leading principal minors to check is n - m, in which
n = number of variables, m = number of equality constraints. kmin = 2*m+1, kmax = m+n,
In ascending order, these are the determinants (leading principal minors)

      LPMk ==  LPM (BH, k)  with k = kmin, kmin+1,...kmax, 
               in which kmin = 2*m+1, kmax = m+n.
If it happens that kmin = kmax (ie., n - m = 1), then there is only one determinant to find: 
LPM (BH, m+n) which is the same as determinant (BH).

However, LPM (BH, m+n) is always the same as determinant (BH).
-----------------------------------------------------------------------------------
IN SUMMARY, for n = length(varList), m = length(gList), k = 2*m+1, 2*m+2,..., m+n,

sufficient conditions (but not necessary) for a relative MAXIMUM are:
     (a)  sgn (det BH) == sgn (LPM (BH, m+n)) = (-1)^n,
     (b)  sgn(LPMk) alternates for successive k's in the range (kmin, kmax),

sufficient conditions (but not necessary) for a relative MINIMUM are:
     (a)   sgn (det BH) == sgn (LPM (BH, m+n)) = (-1)^m, 
     (b)   sgn(LPMk) remains the same as (-1)^m for all k (no sign oscillation).

### Chiang-Wainwright, Ch.12, Example 1 

Find the extremum of f = x*y subject to x + y = 6. For this example m = 1 and n =  2. 
Hence kmin =  2*m + 1 = 3, and kmax = m + n = 3, so there is only one leading principal 
minor, LPM3 = LPM (BH, 3) == determinant (BH), of the bordered Hessian matrix BH to check.

For this problem, the sign oscillation criteria play no role.

(-1)^m = (-1)^1 = -1, so we have a local minimum if LPM3 < 0.

(-1)^n = (-1)^2 = +1, so we have a local maximum if LPM3 > 0.

f : x*y;
g : x + y;

### Graphical Exploration

In the x-y plane, a solution (critical) point (x,y) must lie on the line x + y = 6. For what value
of the objective function f = x*y can this occur? We first draw the line x+y = 6 (y = 6 - x) in
blue, and then plot the locus of the points f = x*y for both f = 5 and f = 10 for orientation.
For the latter two curves, we use draw2d's function 
           implicit (F(x,y) = C, x, x1, x2, y, y1, y2),
which tries to draw all points (x,y) for which F(x,y) = C in the ranges provided.

xmax : 6$ ymax : 6$
wxdraw2d ( xlabel = "x", ylabel = "y", 
    title = "max f = x*y such that x + y = 6",
    key = " x + y = 6", explicit (6 - x, x, 0, xmax), 
    color = purple, key = " f = 10", implicit ( f = 10, x, 0, xmax, y, 0, ymax),
    color = red, key = " f = 5", implicit ( f = 5, x, 0, xmax, y, 0, ymax))$
    

This graphical exploration suggests perhaps f = 9, with x = 3, y = 3 roughly. We could
add another implicit curve corresponding to f = 9 to tie this down. Instead we use our
friend optimum (func, varList, constraint).

### optimum (f, [x,y], 6 - g)

optimum (f, [x,y], 6 - g);

Since LPM3 = 2  > 0 we have a relative maximum.

optimum creates a global list cp.

cp;

cp : cp[1];

at (f, cp);

Show the curve f = 9 together with the constraint line x + y = 6 (i.e.,  y = 6 - x).

xmax : 6$ ymax : 6$
wxdraw2d ( xlabel = "x", ylabel = "y", 
    title = "max f = x*y such that x + y = 6",
    key = " x + y = 6", explicit (6 - x, x, 0, xmax), 
    color = red, key = " f = 9", implicit ( f = 9, x, 0, xmax, y, 0, ymax),
    color = magenta,  key = "y = 3", explicit (3, x, 0, 3),
    color = purple,  key = "x = 3", parametric (3, yy, yy, 0, 3),
    color = black,  key = "", points ([ [3,3] ]))$

### Do it "by hand"...

Form the Lagrangian function 
     L =  f + λ (6 - g)  =   f + λ (6 - x - y), 
and look for the simultaneous solutions to the three first derivatives (wrt x, y, and λ) set 
equal to zero (i.e., the critical points). The jacobian function offers a fast route to a list of 
these first derivatives, although you have to use stuff[1] to get the list we want. 

[f, g];

L : f + λ*(6 - g);
gradL : jacobian ([L], [x, y, λ ] ) [1];
solns : solve (gradL, [x, y, λ]);

soln : solns[1];

The same solution as found by optimum.

at (f, soln);

The bordered Hessian matrix BH for one equality constraint g(x1, x2,,,xn) can be calculated 
using Bhessian (Lagrange-func, g, [x1,x2,...,xn]). Bhessian is defined in Econ1.mac.

BH : Bhessian (L, g, [x,y]);

For this example, BH (Bhessian) is already purely numerical, so we can skip the step of 
evaluating it at the critical values of x, y, and λ. The second derivatives of the Lagrangian 
function L (wrt x and y) appear in the 2 x 2 block in the lower right corner. The first 
derivatives of g wrt x and y appear on the borders.

With m = 1 and n = 2,  (-1)^m = -1 and (-1)^n = +1.  2*m + 1 = 3, m + n = 3, so the 
single leading principal minor we need to check is LPM3 ==  LPM (BH, 3).

LPM (BH, 3);

Since LPM3  = 2  is positive and has the same sign as (-1)^n, the critical point found
(x = 3, y = 3) is a relative maximum.

Our Maxima function BHtest (BH, m, n) returns a list of all relevant leading principal minors
of the bordered Hessian matrix BH, and if the relevant criteria set out above are met,
will print out the nature of the critical point used to generate BH.

BHtest (BH, 1, 2);

LPM (BH,3) is defined as the determinant of the submatrix of BH formed by deleting all but
the third row and all but the third column, but since BH is a 3 x 3 matrix this submatrix is 
just BH. So we could have just calculated the determinant directly.

determinant (BH);

Dowling makes use of the language: "kth bordered leading principal minor" , which can
be calculated using our Maxima function BLPM (BH, m, k), where m = number of constraints,
and k takes values from  kmin = m + 1 up to kmax = n. 
For our case, m = 1, kmin = 2, kmax = n = 2, so there is just the determinant 
BLPM (BH, 1, 2) --> LPM (BH, 3).

In this paragraph, let Hb stand for "H bar", and Hb2 for the "second bordered leading 
principal minor". Dowling uses the overbar symbol (over H) for his bordered Hessian matrix, 
but in the absence of an overbar symbol in wxMaxima I will use Hb. Then Dowling's notation 
is |Hbk| for the kth bordered leading principal minor. BLPM is defined in Econ1.mac.

Hb2 : BLPM (BH, 1, 2);

### Ch. 12,  Example 5

Optimize the function f = 4 x^2 + 3 x y + 6 y^2 subject to the constraint x + y = 56.

g : x + y;
f : 4*x^2 + 3*x*y + 6*y^2;

Here m = 1, n = 2, m+n = 3, 2*m+1 = 3, check the sign of the leading principal minor LPM3.

(-1)^m = (-1)^1 = -1, so local minimum if LPM3 < 0.

(-1)^n = (-1)^2 = +1, so local maximum if LPM3 > 0.

### Graphical Exploration

Given that x and y have to satisfy x + y = 56 (and live on the blue line below)
we are looking for the point that will yield either a minimum or a maximum value
for f. We plot some f = constant curves using draw2d's implicit function.

wxdraw2d ( xlabel = "x", ylabel = "y",
   key = " x + y = 56", explicit (56 - x, x, 0, 56), 
    color = purple, key = " f = 18,000", implicit ( f = 18000, x, 0, 56,y,0,56),
    color = red, key = " f = 12,000", implicit ( f = 12000, x, 0, 50,y,0,50),
    color = green, key = " f = 10,000", implicit ( f = 10000, x, 0, 40,y,0,40),
    color = black, key = " f = 9,000", implicit ( f = 9000, x, 0, 40,y,0,40))$
    

### optimum (f, [x,y], 56 - g)

optimum (f, [x,y], 56 - g);

Since LPM3 < 0 we have a relative minimum at (x = 36, y = 20).

optimum creates the global list cp.

cp;

soln : cp[1];

soln : at (soln, lam[1] = λ);

at (f, soln);

### Do it "by hand"....

[f, g];

L : f + λ*(56 - g);
gradL : jacobian ([L], [x, y, λ])[1];
solns : solve (gradL, [x, y, λ]);

soln : solns[1];

BH is the bordered Hessian matrix appropriate to this problem.

BH : Bhessian (L, g, [x,y] );

Since BH is already purely numerical, we again skip the step of evaluating BH at the 
critical point revealed by the first order conditions (FOC).

This is another case of m = 1 equality constraints and n = 2 variables. So there is only one
leading principal minor to check, LPM3 ==  LPM (BH,3).

BHtest (BH, 1, 2);

LPM (BH,3);

determinant (BH);

Since m = 1, (-1)^m = -1, and LPM (BH, 3) < 0, the sufficient conditions for a relative 
minimum have been met.

### Prob. 12.28, Optimization of  f(x, y, z) with one constraint

Use the bordered Hessian to check the second order conditions in Prob. 5.12(c), where
f = 4 x y z^2 was optimized subject to x + y + z = 56. This is a m = 1, n = 3 problem.
kmin = 2*m + 1 = 3, kmax = m + n = 4, check LPM(BH,3) and LPM(BH,4). 
(-1)^m = -1; local minimum if LPM4 < 0, LPM3 < 0.
(-1)^n = -1; local maximum if LPM4 < 0, LPM3 > 0.

A  "graphical exploration" is difficult for n = 3 variables. So proceed to the bordered Hessian
test.

f : 4*x*y*z^2;
g : x + y + z;
optimum (f, [x,y,z], 56 - g);

cp;

The only physical solution which has x, y, z, f all positive is the second solution in this list.

soln : cp[2];

soln : at (soln, lam[1] = λ);

at (f, soln);

LPM4 < 0, and LPM3 > 0,  so sufficient conditions for a relative maximum.

"By hand" solution steps: start with FOC (first order conditions).

[f, g];

L : f + λ*(56 - g);
gradL : jacobian ([L], [x, y, z, λ])[1];
solns : solve (gradL, [x, y, z, λ]);

The solutions for which f = 0 and λ = 0 do not involve the enforcement of the constraint, 
so the only physical solution is the second list.

soln : solns[2];

cp : rest (soln, -1);

The value of the objective function f at the given critical point is the same as the value of
the Lagrangian function L at the critical point, since the expression which multiplies λ is 
zero at the critical point.

at (f, cp);

at (L, cp);

Second order conditions (SOC).
Construct the bordered Hessian matrix using Bhessian

BH : Bhessian(L, g, [x, y, z]);

Convert to a purely numerical matrix by evaluating BH at the critical point. We can use
either soln or cp for this, since the symbol λ does not appear in this matrix.

BH : at (BH, cp);

Call BHtest (BH, m, n) for m = 1 constraint, and n = 3 variables.

BHtest (BH, 1, 3);

With n = 3, (-1)^n = -1, det BH < 0, signs alternate, so sufficient condition for a relative
maximum.

### Ch. 12, Example 6,  Cobb-Douglas Model

Extend the analysis of Example 10 in Dowling Ch. 6, Sec 9, to determine whether the
critical point found is a relative maximum or minimum.

The Cobb-Douglas model can be used to relate the quantity q of units which can be
produced as a function of the number of units of capital K and the number of units of
labor L,
        q = A*K^α*L^β
in which α > 0 is the "output elasticity of capital", 0> β < 1 is the "output elasticity of labor",
and A > 0 is an "efficiency parameter" measuring the level of technology employed.
With PK the price per unit of capital K and PL the price per unit of labor L, the cost
of production of q units is
       cost = PK*K + PL*L.

Assume the budget constraint is $108 and PK = 3$/unit of capital, PL = 4$/unit of labor,
and we want to optimize q = K^(0.4)*L^(0.5) subject to the constraint that
     3*K + 4*L = 108.

This is a m = 1, n = 2 problem.  2 m + 1 = 3, m + n = 3, so only LPM(BH,3) needs to be
checked. (-1)^m = -1 so we have a local minimum if LPM3 < 0.
(-1)^n = +1, so we have a local maximum if LPM3 > 0.

### optimum (q, [K, L], 108 - g)

solve will be more likely to find solutions if we simplify the form of q using ratsimp.

[K, L];

kill(L)$
g : 3*K + 4*L;
q : ratsimp (K^(0.4)*L^(0.5));

optimum (q, [K,L], 108 - g);

cp;

soln : cp[1], numer;

Kmax : at (K, soln);
Lmax : at (L, soln);

qmax : at (q, soln), numer;

sconcat ("q = ", qmax);

grind(%);

string(qmax);

wxdraw2d ( xlabel = "K", ylabel = "L",
   key = " 3*K + 4*L = 108", explicit ((108 - 3*K)/4, K, 0, 36),     
    color = red, key = sconcat ( " q = ", qmax), implicit ( q = qmax, K, 0, 36, L, 0, 36),
    color = brown, key = "L = 15",  explicit (Lmax, K, 0, Kmax),
    color = purple, key =  "K = 16", parametric (Kmax, yy, yy, 0, Lmax),
    color = black, key = "", points ([ [Kmax, Lmax] ]))$

### Do it "by hand"....

We include the budget constraint by defining the "Lagrange function"
    Q = q + λ ( 108 - 3 K - 4 L),
in which λ is the "Lagrange multiplier", and we are constraining the optimization
by requiring that g = 3 K + 4 L = 108. We then solve for the values of K, L, and λ for
which Q is optimized. 

We need to use a symbol different from L for the Lagrange function, since we are using L 
for one of the search variables in the objective function q and the constraint g.
Let Q be the Lagrange function.  We previously used L as the symbol
for the Lagrangian expression, hence use of kill(L) above.

This is a m = 1, n = 2 problem.  2 m + 1 = 3, m + n = 3, so only LPM(BH,3) needs to be
checked. (-1)^m = -1 so we have a local minimum if LPM3 < 0.
(-1)^n = +1, so we have a local maximum if LPM3 > 0.

[q,g];

Q : q + λ*(108 - g);
gradQ : jacobian ([Q], [K,L,λ] )[1];
solns : solve ( gradQ, [K, L, λ] );

soln : solns[1];

Let cp be the critical point (K,L) found.

cp : rest (soln, -1);

Both the objective function q and the Lagrangian Q should have the same value at the 
critical point cp, since the expression multiplying lam is zero at cp.

subst (cp, [q, Q] ), numer;

Second order conditions (SOC):

BH : Bhessian (Q, g, [K, L]);

Before checking LPM's, evaluate at critical point cp.

BH : at (BH, cp);

Let's try out BHtest using BH in the above form.

BHtest (BH, 1, 2);

float(%);

Now we convert BH to floating point numbers and then retry BHtest.

BH : float (BH);

BHtest (BH, 1, 2);

With n = 2, (-1)^n = +1, and det BH > 0, so we have met the sufficient conditions for
a relative maximum.

### Leydold-Petry Example: Cobb-Douglas

This example is from the pdf: Introduction to Maxima for Economics, Sec. 9.7.

Suppose we have to maximize a Cobb-Douglas production function
  Y = K*L^2
under  the constraint
  g (K, L) = K + L = 3.

[K, L];

Y : K*L^2;
g : K + L;

### Graphical Exploration

wxdraw2d ( xlabel = "K", ylabel = "L",
   key = "K + L = 3", explicit (3 - K, K, 0, 3), 
    color = purple, key = " Y = 5", implicit ( Y = 5, K, 0, 3, L, 0, 3),
    color = red, key = " Y = 2", implicit ( Y = 2, K, 0, 3, L, 0, 3),
    color = brown, key = " Y = 0.5", implicit ( Y = 0.5, K, 0, 3, L, 0, 3))$

### optimum (Y, [K, L], 3 - g)

optimum (Y, [K, L], 3 - g);

optimum has found a maximum at K = 1, L = 2, with Y = 4, which agrees with our graphical
exploration, and has also found a local minimum at K = 3, L = 0, Y = 0, but we are only 
interested in maximizing the production function Y at some L > 0 and Y > 0 value.

### Do it "by hand"...

Let F be the Lagrangian function.

F : Y + λ*(3 - g);
gradF : jacobian ([F], [K,L,λ] )[1];
solns : solve ( gradF, [K, L, λ] );

soln : solns[1];

at (Y, soln);

Construct the bordered Hessian matrix using the Lagrange function F.

BH : Bhessian (F, g, [K,L] );

Evaluate the bordered Hessian matrix at the critical point cp of interest.

BH : at (BH, soln);

BHtest (BH, 1, 2);

With m = 1, n = 2, 2*m + 1 = 3, m + n = 3, only one LPM to check. (-1)^n = +1, so
for a relative maximum we need LPM3 > 0.

Leydold and Petry (by error) define Y : K^2*L in line 1 of their Maxima solution, and
get the solution L = 1, K = 2 after that error.

### Francisco Feri Example, log transformation of objective func.

Francisco Feri, Chair, Dept. of Economics, Royal Holloway, Univ. of London, has
posted some pdf's on the web dealing with optimization problems.

https://intranet.royalholloway.ac.uk/economics/documents/
   pdf/courseformsandinfo/ec5555lecture6(2013)-constrainedoptimsation.pdf

Maximize f = x^3*y subject to the constraint x + y = 6.

[x, y];

f : x^3*y;
g : x + y;

### Graphical Exploration

wxdraw2d ( xlabel = "x", ylabel = "y",
   key = "x + y = 6", explicit (6 - x, x, 0, 6), 
    color = brown, key = " f = 50", implicit ( f = 50, x, 0, 6, y, 0, 6),
    color = red, key = " f = 100", implicit ( f = 100, x, 0, 6, y, 0, 6),
    color = purple, key = " f = 200", implicit ( f = 200, x, 0, 6, y, 0, 6))$  

### optimum (f, [x,y], 6 - g)

[f, g];

optimum (f, [x, y], 6 - g);

cp;

cp1 : cp[1], numer;

fmax : at (f, cp1);

wxdraw2d ( xlabel = "x", ylabel = "y",
    title = " max x^3*y s. t. x + y = 6",
   key = "x + y = 6", explicit (6 - x, x, 0, 6), 
    color = purple, key = sconcat (" f = ", fmax), implicit ( f = fmax, x, 0, 6, y, 0, 6))$

If we instead work with the natural log of f = x^3*y, and 
define    f1 = ln (x^3*y) = ln(y) + 3*ln(x),

logexpand : all$
f1 : log (f);
logexpand : true$

optimum (f1, [x, y], 6 - g);

We see that the relative maximum of log(f) is at the same critical point (4.5, 1.5) as the
maximum of f.

"BY Hand solution" using f1:

[f1, g];

L1 : f1 + λ1*(6 - g);
gradL1 : jacobian ([L1],[x, y, λ1] )[1];
solns_ln : solve (gradL1, [x, y, λ1]);

soln_ln : solns_ln[1];

BH : Bhessian (L1, g, [x,y]);

BH : at (BH, soln_ln);

BHtest (BH, 1, 2);

### Chiang-Wainwright, p. 360, Two Variable Symbolic Example

Consider a simple two-period model where a consumer's utility is a function of 
consumption in both periods. Let the consumer's utility function be
   U (x1, x2) = x1*x2,
where x1 > 0 is consumption in period 1 and x2 > 0 is the consumption in period 2. 
The consumer is also endowed with a budget B at the beginning of period 1. 
Let r > 0 be the market interest rate at which the consumer can choose to borrow or
lend across the two periods. The consumer's intertemporal budget constraint is that
x1 and the present value of x2 add up to B. Thus
   x1 + x2/(1+r)  = B.
Find the values of x1 and x2 which maximize the utility subject to this constraint.

[x1,x2,r,B,λ];

U : x1*x2;
g : x1 + x2/(1+r);
L : U + λ*(B - g);
gradL : jacobian ([L], [x1, x2, λ])[1];

### optimum (U, [x1, x2], B - g)  Symbolic solution

If we don't give Maxima information about r and B, optimum returns an "indefinite"
solution, whereas the solution is actually a relative maximum.

optimum(U, [x1, x2], B - g);

But if we give Maxima information that both r and B should be considered positive
numbers, optimum (and BHtest) find a relative maximum.  Use:    assume( r > 0, B > 0)$
Note: to remove assumption, use forget (r > 0, B > 0)$. To check on status of such
assumptions, use facts();

assume (r > 0, B > 0)$
optimum(U, [x1, x2], B - g);
forget (r>0, B >0)$
facts();

optimum creates the global list cp.

cp;

soln : cp[1], factor;

soln : at (soln, lam[1] = λ);

at (U, soln);

The constraint term doesn't contribute to L at the critical point, since B - g = 0 there.

at (L, soln);

at (B - g, soln);

"By hand" SOC:

BH : Bhessian (L, g, [x1,x2] );

With one equality constraint, m = 1, and 2 variables (x1, x2), n = 2, 2*m+1 = 3, m+n = 3
and we only need to check the sign of the determinant of BH, which is the same as LPM3.

determinant (BH);

Since (-1)^n = (-1)^2 = +1 and det BH  > 0 (with the interest rate r > 0), the critical point
found is a relative maximum.

assume (r >0)$
BHtest (BH, 1, 2);
forget (r > 0)$

### Another Two Variable Symbolic Example

Maximize f(x,y) = x^a y^b subject to  x + y = 10, where a > 0, b > 0, and f(x,y) is defined
on the set (x,y) with x >= 0, y >= 0.

Step 1. Find all stationary points (critical points) of the Lagrangian.

[x,y,a,b,λ];

f : x^a*y^b;
g : x + y;
L : f + λ*(10 - g);
gradL : jacobian ([L], [x, y, λ])[1];
soln : solve (gradL, [x,y,λ])[1];

cp : rest (soln, -1);

With a > 0 and b > 0, the solutions x = xs and y = ys are positive.

fmax : at (f, cp);

With a > 0 and b > 0, fmax > 0 at the critical point cp.

Next look at the secondary conditions needed for a relative maximum.

BH : Bhessian (L, g, [x, y]);

Evaluate BH at the critical point found.

BHcrit : at (BH, cp);

With one constraint, m = 1, and two variables, n = 2, 2*m+1 = 3, n + m = 3, we only have
 LPM3 = det(BH) to check.

determinant (BHcrit);

expand (%);

ratsimp (%);

With a, b > 0, the determinant of BH (evaluated at the critical point) is positive, and since
 (-1)^n = (-1)^2 = 1 > 0, we have a sufficient condition for a relative maximum.

### Nonlinear 3 Variables Example,  2 Equality Constraints

Maximize f = x*y*z subject to both: x^2 + y^2 = 4 and x + z = 2. The first constraint is 
nonlinear, which adds a layer of challenge when finding the critical point(s).

[x, y, z];

f : x*y*z;
g1 : x^2 + y^2;
g2 : x + z;

This is a m = 2, n = 3 problem: 2*m+1 = 5, m + n = 5. Only one leading principal minor to
check, LPM5. 
(-1)^n = -1. If LPM5 < 0, then we have a relative maximum. 
(-1)^m = +1. If LPM5 > 0, then we have a relative minimum. 
If LPM5 = 0, then we have an indefinite case.

First we optimistically try optimum (func,varL, gL):

optimum (f, [x, y, z], [4 - g1, 2 - g2]);

One solution x = 2, y = 0, z = 0, f = 0 was found, which satisfies the constraints, but we 
suspect the Maxima function solve has difficulties with this problem because of the 
nonlinearity. Let's revert to a "by hand" mode of using Maxima.

[f, g1, g2];

L : f + λ1*(4 - g1) + λ2*(2 - g2);
gradL : jacobian ([L], [x, y, z, λ1, λ2] )[1];
solns : solve (gradL, [x, y, z, λ1, λ2] );

The system of equations to solve are nonlinear, so the Maxima function solve (advertised
for solving polynomials) needs help. We can attempt  to reduce the number of variables  
and equations to achieve a solution. Give the 5 first derivatives of L (in the list gradL) 
the names e1, e2, ..., e5.

[e1,e2,e3,e4,e5] : gradL$
matrix (["e1","e2", "e3", "e4", "e5"], 
           [ e1,   e2,    e3,    e4,   e5 ]);

We can simplify by solving e2 and e3 for λ1 and λ2

λ1s : solve (e2, λ1);
λ2s : solve (e3, λ2);

Solve has assumed y is nonzero.
Let's combine these into one list we call λs.

λs : flatten ( cons ( λ1s, λ2s));

An alternative route would be: λs : [ λ1s[1], λ2s[1] ];

Now replace λ1 and λ2 in e1, preserving the name e1.

e1 : at ( e1, λs );

The equations to solve are now e1 = 0, e4 = 0, e5 = 0. We can multiply e1 by y without 
changing the equation (e1*y = 0 has the same soln as e1 = 0 if we restrict y to be nonzero).
 We let e1 be the name of this product.

e1 : y*e1, expand;

matrix ( ["e1", "e4", "e5" ], [e1, e4, e5] );

### Using eliminate (eqnL, exprL)

Both e1 and e4 now have terms with y^2 as a factor. Both e1 and e5 contain the factor z.

If we eliminate y^2 and z, using these three equations, we will get an equation in x only.
We can use the Maxima function eliminate ([eqn_1, …, eqn_n], [x_1, …, x_k]). From the
Maxima manual:

"Eliminates variables from equations (or expressions assumed equal to zero) by taking 
successive resultants. This returns a list of n - k expressions with the k variables 
x_1, …, x_k eliminated. First x_1 is eliminated yielding n - 1 expressions, 
then x_2 is eliminated, etc. If k = n then a single expression in a list is returned free of 
the variables x_1, …, x_k. In this case solve is called to solve the last resultant for the 
last variable."

eliminate ( [e1,e4, e5], [y^2, z]);

ex : %[1];

ex is a cubic equation in x. solve usually can deal with a cubic equation.

xsolns : solve (ex);

Let's assign names to these three solutions for x.

[x1s, x2s, x3s] : xsolns;

Only x2s and x3s imply positive values for x; we reject x1s (which implies a negative x).
We Let's first deal with the x3s ( x = 2) case, and solve for (y, z) for this case. 

Assuming x = 2, the equation e5 = 0 ( - z - x  + 2 = 0 ), implies z = 0. 
Then equationn e4 = 0  (- y^2 - x^2 + 4 = 0 ), implies y^2 = 0 which implies y = 0, 
consistent with eqn e1 = 0, so (x = 2, y = 0, z = 0) and f = 0, a solution we reject.

Repeating this argument using Maxima:

at (e5, x3s);

Then e5 = 0 implies z = 0, hence f  = 0.

at (e4, x3s );

Then e4 = 0 implies y = 0. We reject the x = 2 soln x3s.

Next look at the case x2s and work with floating point numbers.

x2s : float (x2s);

numer : true$

at (e5, x2s);

solve(%,z);

z2s : %[1];

e4;

at (e4, x2s);

solve (%, y^2);

fac : sqrt(rhs (%[1]));

So y = +/- 1.8016, and we choose the mathematical solution with the positive value of y
to correspond to a real world Economics solution.

y2s : y = fac;

cp2 : [x2s, y2s, z2s];

at (f, cp2);

Thus by helping out solve (using eliminate) we have uncovered a solution (critical point)
for which all the variables are positive. Next, find the values of the Lagrange multipliers.

λs;

λs :  at (λs, cp2);

soln : flatten (cons (cp2, λs));

[L, x, y, z, g1, g2];

Set up the appropriate bordered Hessian matrix:

BH : Bhessian (L, [g1, g2], [x, y, z] );

Use the replacement list soln to evaluate BH at the critical point found:

BH : at (BH, soln);

In the lower right block we have the 3x3 submatrix of the second order partial derivatives of
the Lagrangian function wrt (x,y,z), bordered by the first order partial derivatives of g1 and
g2, all evaluated at the critical point.

This is an m = 2, n = 3 problem, use BHtest (BH, m, n).

BHtest (BH, 2, 3);

With m = 2, n = 3, 2 m + 1 = 5, m + n = 5, there is only one LPM to check, LPM5 = det(BH). 
(-1)^n = (-1)^3 = -1,  and LPM5 < 0 so a relative maximum.

determinant (BH);

numer : false$

## Input-Output Analysis   [12.6] 

We roughly quote from the beginning of Sec. 5.7 of Chiang & Wainwright:

In its "static version", the input-output analysis of Professor Wassily Leontief, a Nobel
Prize winner, deals with the question: "What level of output should each of the n
industries in an economy produce, in order that it will just be sufficient to satisfy the
total demand for that product?"

In a modern economy where the production of one good requires the input of many other
goods as "intermediate goods" in the production process (steel requires coal, iron ore,
electricity, etc.), *total demand* x for product i (x[i]) will be the summation of all intermediate
demand for the *production* of the product plus the "final demand" b for the product
arising from consumers, investors, the government, and exporters, as ultimate users of
the product.

To simplify the problem, the assumptions made are
 1. each industry produces only one homogeneous commodity,
 2. each industry uses a fixed "input ratio" for the production of its output,
 3. production in each industry is subject to "constant returns to scale", so that a k-fold
       change in every input will result in an exactly  k-fold change in the output

From these assumptions, in order to produce each unit of the jth commodity, the need
for the ith commodity must be a fixed amount, which we denote by a[i,j]. Specifically,
the production of each unit of the jth commodity will require a[1,j] (amount) of the
first commodity, a[2,j] amount of the second commodity,..., and a[n,j] of the nth
commodity. The order of the indices is easy to remember: the first index refers to the
input, the second to the output, so a[i,j] indicates how much of the ith commodity is
used for the production of each unit of the jth commodity.

If we assume prices are a given, we can use "a dollar's worth" of each commodity as
its unit. Then the statement  a[3,2] = 0.35 means that 35 cents worth of the third
commodity is required as an input for producing one dollar's worth of the second
commodity.

For an n-industry economy, the input coefficients can be arranged into an n x n matrix A,
in which each column specifies the input requirements for the production of one unit (one
dollar's worth) of the output requirements for a particular industry.

Here is an example for a four industry economy:

A : genmatrix (a, 4, 4);

For this 4 x 4 matrix A, the second column states that to produce a unit (one dollar's
worth) of commodity 2, the inputs are: a[1,2] units of commodity 1, a[2,2] units of
commodity 2, etc. If no industry uses its own products as an input, then the elements in
the principal diagonal of matrix A will be zero.

### The Open Model

If the n industries constitute the entirety of the economy, then all their products would be
for the sole purpose of meeting the "input demand" of the same n industries (to be used
for further production) as against the "final demand" (such as consumer demand, not to
be used for further production). 

At the same time, all the inputs used in the economy would be in the nature of 
"intermediate inputs" (those supplied by the n industries) as against "primary inputs" 
(such as labor and capital, not industrial products). To allow for the presence of 
final demand (consumer demand) and primary inputs to production (labor, capital, etc), 
we must include in the model an "open sector" outside of the n-industry
network. Such an open sector can accommodate the activities of the consumer
households, the government sector, and even foreign countries.

In view of the presence of the open sector, the sum of the elements of each column of
the input matrix A must be less than 1. Each column sum represents the "partial
input cost" (not including the cost of the primary inputs (labor, capital,...).

Symbolically we need (using Maxima notation):  sum (a[i,j], i, 1, n) < 1, for j = 1, 2,...,n.

If a[i,j] is a "technical coefficient" expressing the value of input from industry i required 
to produce one dollars worth of product by industry j, the *total demand* for product from
industry i is x[i], which can be expressed at the sum of b[i] (final demand by consumers) 
plus the intermediate demands of the n-industry network.

x[i] = a[i,1]*x[1] + a[i,2]*x[2] + ... + a[i,n]*x[n] + b[i],

b[i] is the "final demand" (by consumers, etc) for the product of industry i.

If we move all terms proportional to x[k] to the left of the equal sign, leaving the
exogenous variables b[i] on the right hand side, we can write this set of n equations in
matrix form:

      (I - A) . X = B,  or
      M . X = B

where M = (I - A), with I standing for the unit matrix, A is the matrix of technical
coefficients, X is a matrix column vector of unknown total demands x[k], and 
B is a matrix column vector of given final demands b[k].

The matrix M =  (I - A) is called the Leontief matrix. As long as the Leontief matrix
is nonsingular, we can find its inverse M^^(-1) = invert (M), and obtain the unique
solution of the n equations using  Xs : invert(M) . B

For our n = 4 example we first construct lists xL and bL and turn each into a matrix
column vector using X :  cvec (xL), for example, where cvec is defined in Econ1.mac.
ident(4) returns a 4 x 4 unit matrix (1's on the diagonal and 0's elsewhere).

A : genmatrix (a, 4, 4);
xL : []$
bL : []$
for k thru 4 do (
    xL : cons (x[k], xL),
    bL : cons (b[k], bL))$
xL : reverse (xL);
bL : reverse (bL);
X : cvec (xL);
B : cvec (bL);
M : ident(4) - A;

meqns :  M .  X = B;

The inverse of a matrix M is returned by invert (M), and we know that
     invert (M) . M  = I
returns a unit matrix (don't try to verify this using the above definition of M).

Hence the solution matrix column vector Xs is given by
  Xs :  invert (M) . B
(again, don't try to do this symbolically now). See the numerical example below.

### A Numerical Example

Chiang and Wainwright, p. 115, have the following numerical example for a three
industry economy. A is the matrix of technical coefficients. Element (i,j) of A can be
selected in Maxima using the syntax A[i,j].

A : matrix ( [0.2, 0.3, 0.2], [0.4, 0.1, 0.2], [0.1, 0.3, 0.2] );

A[1,2];

The value of A[1,2] = 0.3 means that $0.30 worth of product 1 is required as input for 
producing $1 worth of product 2. If we add up the elements of column 2 of the matrix A
we get

sum (A[k, 2], k, 1, 3);

which is less than 1 as it should be (If it required $2 worth of input to produce $1 worth
of product 2, production would promptly cease, in a rational world). With this given data,
we assume the difference between 0.7 and 1 means that "primary input" such as the cost
of labor and capital amounts to $0.30 dollar's worth for every one dollar's worth of product 2.

To produce one billions dollar's worth of product 2 requires 0.3 billion dollar's worth of
product 1.

### Leontief(A)

The software file Econ1.mac defines a Maxima function Leontief (A) which returns
the Leontief matrix (I - A).

Leontief (A);

length (A);

Here we are creating the Leontief matrix M = (I  - A) "by hand".

M : ident (3)  - A;

length (M);

rank (M);

The rank (see the first section of Dowling12A.wxmx) of M is 3, which is equal to the
dimension of M (given by length(M)), so M is nonsingular, and the matrix inverse of 
M exists.

Minv : invert (M);

If the specific (consumer) final-demand column vector (say, the final-output target of 
a development program) happens to be (in billions of dollars) 

B : cvec ([10, 5, 6] );

length (B);

### InputOutput (M, B)

The software file Econ1.mac defines the Maxima function 
                InputOutput (M, B)
in which M = I - A is the Leontif matrix and B is the "final demand". 
InputOutput (M, B) checks some properties of M and returns the total demand matrix 
column vector solution.

InputOutput (M, B);

Here we work this problem "by hand": the total demand solution vector 
Xs (again in billions of dollars) is

Xs : Minv . B;

Xs[1];
Xs[1,1];
Xs[2,1];

We see that we can use syntax Xs[1,1] to get the first element of the total demand 
column vector matrix, Xs[2,1] to get the second element, etc. Alternatively, we can use 
Xs[1][1], Xs[2][1], etc.

Let a0(j) be the dollar amount of the primary input (like labor, capital, etc.) used in 
producing  one dollar's worth of the jth commodity. a0(j) can be calculated by 
subtracting the sum of the jth column of the input matrix A from 1. 
In general, A[k, j] returns the kth row and jth column element of the matrix A.

Here we define a Maxima function depending of j:

a0(j) := 1 - sum (A[k, j], k, 1, 3)$

then we have

[a0(1), a0(2), a0(3)];

### Pinput (A, j)

The Maxima function Pinput (A, j) (primary input j)  is defined in Econ1.mac, and does the
same job as the Maxima function a0(j) if A is the given matrix of technical coefficients, 
plus figures out what the final k value should be in the sum (using the dimension of A).

Pinput (A,1);

The Maxima lambda expression can be used to map Pinput onto a list of numbers for
a chosen value of A. See the Maxima manual index under lambda.

map (lambda ([xx], Pinput (A, xx)), [1, 2, 3] );

An important question now arises. The production of the total demand output mix 
Xs[k] (for k = 1 thru 3) must entail a definite required amount of the primary input 
(labor, capital, etc).   Would the amount *required* be consistent with what is *available* 
in the economy? 

On the basis of the values we found for a0(j) (with j = 1, 2, 3), the required 
primary input (PI) (labor, capital, etc) is the sum of the products a0(k)* Xs[k] with k = 1,2,3.

As an example, for product 1

Xs[1];
a0(1);
a0(1) * Xs[1,1];

The required primary input (PI) is:

PI : sum (a0(k)*Xs[k,1], k, 1, 3);

### PItot (A, Xs)

The Maxima function PItot (A, Xs) calculates this required primary input as a function
of the matrix of technical coefficients A and the solution vector Xs.

PItot (A, Xs);

so the specific demand (in billions) assumed in our definition of the consumer final 
demand matrix column vector B willl be feasible if and only if the available amount 
of the primary input (labor, capital,...) is at least $21 billion.

One notable feature of the previous analysis is that, as long as the input coefficients
(technical coefficients) A[i,j] = a[i,j] remain the same numbers, the matrix inverse Minv 
remains the same, even if we consider a hundred or a thousand different final-demand 
vectors - such as a spectrum of alternative development targets.

### Problem 12.33 & 12.34

Given the interindustry transaction demand table in millions of dollars below, find the
matrix of technical coefficients A. The "Sector of Origin" categories refer to inputs 
provided to "output" industries (Sector of Destination). Sector of origin category
 "primary input" refers to input of labor, capital, etc, in millions of dollars,
 required to produce the gross product by each industry (column). Dowling refers
to primary input (in his table) as "Value Added".

print ("                           Sector of Destination ")$
Mtable : matrix (["Sector of Origin", "steel","coal","iron","Auto","Final-Demand","Total-demand"],
                        ["steel",             80,     20,    110,    230,       160,                 600],
                        ["coal",             200,    50,      90,    120,       140,                 600],
                        ["iron",              220,   110,    30,       40,        0,                   400],
                        ["Auto",             60,    140,    160,   240,       400,               1000],
                        ["primary-input",40,    280,      10,   370,       " ",                    " "   ],
                        ["Gross Prod.", 600,  600,     400,  1000,     " ",                    " "   ] );

To make this table clearer, consider the column Auto production under the category Sector
of Destination. Let M stand for million.

To produce $1,000M worth of autos, the Auto producer needs as inputs,
$230M worth of steel, $120M worth of coal, $40M worth of iron, $240M worth of its own 
autos, and $370M worth of "primary inputs" (labor, capital, etc).

If we look at the "Sector of Origin" row Auto, which has the amounts of the Auto industry
inputs to various producers, Auto supplies $60M worth of autos to the steel producer, 
$140M worth of autos to the coal producer, $160M worth of autos to the iron producer, 
$240M worth of autos to itself, adding up to $600M worth of autos. If we then add 
$400M worth of autos to consumers (final-demand), we get a total demand of $1,000M 
worth of autos (one billion dollars worth).

Note that "gross production" of autos is equal to "total-demand" for autos, etc.

In the technical coefficient matrix A, element A[1,1] (row 1, column 1) is the amount of
steel needed as input to produce one dollars worth of steel. From our table we see that
it takes $80M worth of steel input to produce $600M worth of steel output. Since
80/600 = 0.1333, it takes $0.1333 worth of steel input to produce $1 worth of steel output,
so A[1,1] = 0.1333. Likewise, it takes $20M worth of steel to produce $600M worth of coal,
so A[1,2] = 20/600 = 0.033.

A : genmatrix (a, 4, 4);

At this point a[2,3] (for example) is an undefined hash array element, and we proceed
to assign values to these elements based on the matrix Mtable above.

a[2,3];

f(i,j) := float (Mtable[i+1, j+1] / Mtable[7,j+1])$

f(1,1);

for i thru 4 do
     for j thru 4 do  a[i,j] : f(i,j)$

a[1,1];

To bind the symbol A to the numerical value matrix we want, we have to force Maxima
to evaluate A using two (2) single quotes (') in front.

A :  ''A;

The amounts of primary inputs a0(j) required to produce one dollar's worth of the product
of industry j are defined by the Maxima function (for this four industry economy):

a0(j) := 1 - sum (A[k, j], k, 1, 4)$

matrix ( ["steel", "coal", "iron", "Auto"],
    [a0(1), a0(2), a0(3), a0(4)] );

If we multiply these numbers by the gross production (in millions of dollars) for the
respective industry, we get the  primary inputs (labor, capital, etc) required for each 
of the four industries.

matrix ( ["steel", "coal", "iron", "Auto"],
    [a0(1)*600, a0(2)*600, a0(3)*400, a0(4)*1000]);

To perform a consistency check on these values of the elements of the matrix A,
we need to define the column vector X whose four elements are the total demands,
and the column vector B whose four elements are the final demands.

X : cvec ([600, 600, 400, 1000] );
B : cvec ([160, 140, 0, 400]);

We compare A . X and X - B, which should be the same.

A . X;

X - B;

Since A . X = X - B, this is consistent with our original definition of the total demand 
vector X, a solution of the equation (I - A) . X = B, equivalent to:  X - A . X = B.

We can also, of course, invert the Leontief matrix M = (I - A) and the matrix product
Minv . B should reproduce the total demand vector X.

M : ident (4) - A;

rank (M);

Since the rank is the same as the dimension, M is nonsingular, and the inverse
matrix exists.

Minv : invert (M);

X;
Minv . B;

which again shows consistency.

### The Existence of Nonnegative Solutions

We roughly quote the Chiang and Wainwright, p. 116 - 117, section dealing with this
issue.

In our previous examples, the Leontief matrix M = (I - A)  happens to be nonsingular, so
solution values Xs[j] exist, and in our examples turned out to be nonnegative, as economic
sense requires. Nonnegative values for total demand are not automatic, but depend on
the Leontief matrix having specific properties, the "Hawkins-Simon condition", made clear
in a 1949 academic article by David Hawkins and Herbert A. Simon.

The theorem proved by Hawkins and Simon states:

a) Given an n x n matrix M, with M[i,j] <= 0 (i # j) (ie., with all off-diagonal elements non-
positive), and b) a matrix column vector B with all n elements nonnegative, there exists
a n dimensional solution matrix column vector Xs, with all elements nonnegative, if and
only if all of the leading principal minors of M are positive numbers.

The Leontief matrix M = (I - A) satisfies condition (a), and the final demand matrix
column vector B satisfies condition (b). Hence if M is nonsingular, and if all of the 
leading principal minors of M are positive numbers, we are guaranteed that a unique
solution for total demand Xs exists, and that all elements Xs[j] are nonnegative.

### Economic Meaning of the Hawkins-Simon Condition

We continue to follow Chiang and Wainwright (pp. 118 - 119)

Consider the two industry case. We assume a[i,j] >= 0 for all (i, j) = 1, 2.

We want to redefine A as a 2 x 2 matrix with symbolic elements a[i,j]. 
We use remarray (a) to remove our previous numerical binding of a[i,j] elements.
First, we remind the reader of some commonly used array functions.

a[1,4];

arrays;

arrayinfo (a);

listarray(a);

remarray (a);
arrays;

A : genmatrix (a, 2,2);
I : ident (2);
M : I - A;

We use our function LPM (M, k) to compute the leading principal minor of M for k = 1,2.
The function LPM is defined in the file Econ1.mac.

What does the condition LPM (M,1) > 0 require?

lpm1 : LPM (M,1);

lpm1 > 0;

This first condition requires that 1 > a[1,1], or a[1,1] < 1, or that the amount of the first
commodity used in the production one dollar's worth of the first commodity be less than
one dollar (and can be zero).

What does the condition LPM (M, 2) > 0 require?

lpm2 : LPM (M, 2);

lpm2 : expand(lpm2);

lpm2 > 0;

This second condition requires that 
       a[1,1] + a[1,2]*a[2,1] + a[2,2]*(1 - a[1,1])  <  1.
Since the 3rd term on the left hand side is nonnegative, and could be zero, we must have:
       a[1,1] + a[1,2]*a[2,1]  <  1.
In terms of economics, a[1,1] measures the *direct* use of commodity 1 as input in the
production of commodity 1, and the product a[1,2]*a[2,1] measures the *indirect* use
- it gives the amount of commod. 1 needed in producing the specific quantity of 
commod. 2 that is needed in the production of one dollar's worth of commod. 1.
Thus we require that  the amounts of commod. 1 used as both direct and indirect inputs in
producing one dollar's worth of commod. 1, must add up to less than one dollar.

Thus the Hawkins-Simon condition requires practical restrictions on the production
process details that result in a meaningful and viable process, in the sense of Economics.

### Summary of Input/Output Analysis Functions

### InputOutput (M, B)

Note that if the determinant of a matrix is not equal to zero, the matrix is said to be 
"nonsingular", and all its rows and columns are linearly independent, and the rank
of the matrix equals its dimension.

But the determinant of an n dimensional matrix (n x n) is the same as the nth
leading principal minor. If LPM (M, n) > 0, we already satisfy the condition that the 
determinant of M is not equal to zero, so M is nonsingular, and an inverse of M exists, 
and a unique solution of the matrix equation M . X = B exists. 

The Maxima function InputOutput (M, B) simply computes the n leading principal minors
of the square matrix M, and if they are all positive numbers, returns  the unique solution
matrix column vector with all nonnegative elements, otherwise returns "not viable".

The off-diagonal elements of M should be nonpositive numbers.

### Leontief (A)

Leontief (A) returns the Leontief matrix (I - A) and checks that each element of A is positive.

### Pinput (A, j)

Pinput (A, j) subtracts the sum of the jth column elements of the matrix A from 1, 
returning the primary input value (labor, capital, etc) required for the jth industry. 
A is the matrix of technical coefficients.

### PItot (A, Xs)

PItot (A, Xs) returns the total primary input (labor, capital, etc) required, given the matrix 
of technical coefficients A and the total demands solution vector Xs of the equation 
(I - A) . Xs = B, in which B is the matrix column vector of final demands.

## Eigenvalue Tests for Sign Definiteness [12.7]

To this point, the sign definiteness of a Hessian and a quadratic form have been tested by 
using the leading  principal minors. Sign definiteness can also be tested by using the 
characteristic roots of a matrix.  Given a square matrix A, if it is possible to find a vector 
V # 0 and a scalar c such that  A . V = c*V.

The scalar c is called the characteristic root, latent root, or eigenvalue of the matrix A, and 
the vector V is called the characteristic vector, or eigenvector of the matrix A. 

To  find the eigenvalues and eigenvectors of a given matrix A, we note that, using the identity
 matrix I, with the property: I . V = V,   A . V = c*V is equivalent to:
      A . V = c* (I . V)
which can be rearranged as
    (A - c*I) . V  = 0,
where (A - c*I) is called the "characteristic matrix of A". 
 
When using the shorthand way of writing this matrix equation: (A - c*I) . V  = 0, we are
glossing over the notational convenience of writing the matrix on the right hand side of
the equation as 0. Taken literally, this is mathematical nonsense. What is meant (by the
context) is a column matrix of zeros, such as

zeromatrix (4, 1);

if A is a 4 x 4 square matrix. We can also use

cvec ([0, 0, 0, 0]);

With this convenience convention understood, we continue to use the usual matrix 
equation notation M = 0 (caveat emptor).

Since by assumption, V # 0, the characteristic matrix must be singular, and thus its 
determinant must vanish. Setting determinant (A - c*I) equal to zero allows for the 
determination of the eigenvalues and eigenvectors of the matrix A.

If all characteristic roots (c) are positive then A is positive definite,
    else if all c’s are negative then A is negative definite, 
    else if all c’s are nonnegative and at least one c = 0 then A is positive semidefinite,
    else if all c’s are nonpositive and at least one c = 0 then A is negative semidefinite,
    else some c’s are positive and others negative, and A is sign indefinite. 

### Example 8

Let A be the 2 x 2 matrix

A : matrix ( [-6, 3], [3, -6] );

Let I be the 2 x 2 unit matrix

I : ident (2);

To find the characteristic roots (eigenvalues) of the matrix A, the determinant of the
characteristic matrix (A - c*I) is set equal to zero.

c;

eqn : determinant (A - c*I) = 0;

solns : solve( eqn, c);

[c1, c2] : map ('rhs, solns);

Testing for sign definiteness, since both characteristic roots (eigenvalues) are negative,
the matrix A is negative definite.

### Example 9

Continuing with Example 8, the first root c1 = -9 in now used to find the characteristic
vector (eigenvector) of the given matrix A.

[v11, v12];

V1 : cvec ( [v11, v12] );

We replace c by c1 (as defined above in Example 8), and V by V1 in the matrix 
equation (A - c I) . V = 0, we get a matrix equation which implies two ordinary equations.

(A - c1*I) . V1 =  zeromatrix (2,1);

The two equations are identical and require v12 = - v11, so add the normalization condition
transpose (V1) . V1 = 1, which is satisfied if v11 = +/- 1/sqrt(2)

V1 : at (V1, v12 = - v11);

eqn : transpose (V1) . V1 = 1;

solve (eqn, v11);

We can choose either of the two solutions, choose the second.

V1 : cvec ([ 1/sqrt(2), - 1/sqrt(2) ]);

check normalization

transpose (V1) . V1;

A . V1 = c1*V1;

Next, look for the eigenvector V2, using c = c2 as the corresponding eigenvalue of A.

V2 : cvec ( [v21, v22] );

(A - c2*I) . V2 =  zeromatrix (2,1);

Both equations imply we need v22 = v21, and adding the normalization condition gives

V2 : cvec ([ 1/sqrt(2),  1/sqrt(2) ]);

transpose (V2) . V2;

A . V2 = c2*V2;

### eigenvalues (M)

Maxima has a function eigenvalues (M) which returns a list of two lists.
The first list is a list of eigenvalues found, the second list is a list of
multiplicities of the listed eigenvalues.

[eivals, mult] : eigenvalues (A);

eivals;

[c1, c2] : eivals;

c1;

### eigenvectors (M)

Maxima has a function eigenvectors (M) which returns a list of two lists.
The first list is the same as is returned using eigenvalues (M) and has two lists: 
1. a list of eigenvalues found, 2. a list of the multiplicities.
The second list returned by eigenvectors contains a list of eigenvectors (in a list form) 
for each eigenvalue found.

[EI, EV] : eigenvectors (A);

EI;

[eivals, mult] : EI;

eivals;

[c1, c2] : eivals;

c1;

EV;

EV[1];

EV[1][1];

[v11, v12] : EV[1][1];

v11;

EV[2];

EV[2][1];

[v21, v22] : EV[2][1];

At this point eigenvectors(A) has given us un-normalized eigenvector component in a
list form. We can use (for example) our cvec function to convert to a column matrix as 
we ultimately need.

V1 : cvec ([v11, v12]);

V1 is not yet normalized.

transpose (V1) . V1;

V2 : cvec ([v21, v22]);

transpose (V2) . V2;

### Problem 12.43

Use eigenvalues to determine sign definiteness for the matrix

A : matrix ( [4, 6, 3], [0, 2, 5], [0, 1, 3] );

[eivals, mult] : eigenvalues (A);

[c1, c2, c3] : float (eivals);

With all three eigenvalues positive, A is positive definite.

### Problem 12.46 (a)

Find the characteristic roots (eigenvalues) for the matrix A:

A : matrix ( [6, 6], [6, -3] );

[eivals, mult] : eigenvalues (A);

With one eigenvalue negative, and the other positive, A is sign indefinite.

## Binding and Nonbinding Constraints in Economics

### Reference

https://intranet.royalholloway.ac.uk/economics/documents/pdf/
                courseformsandinfo/ec5555lecture6(2013)-constrainedoptimsation.pdf

A constraint is binding if at the optimum the constraint function holds with
equality (sometimes called an equality constraint) giving a boundary solution
somewhere on the constraint itself.

Otherwise the constraint is non-binding or slack (sometimes called an inequality
constraint).

If the constraint is binding we can use the Lagrangian technique.
 
Often we can use our economic understanding to tell us if a constraint is binding.

– Example: a non-satiated consumer will always spend all her income so
the budget constraint will be satisfied with equality.

But in general we do not know whether a constraint will be binding ( = , > or < ).

In this case we use a technique which is related to the Lagrangian, but which is
slightly more general called linear programming, or in the case of non-linear
inequality constraints, non-linear programming or Kuhn-Tucker programming
after its main inventors. 

### Example 1

A firm chooses output x to maximize a profit function π = -x^2 + 10x - 6 .
Because of a labor shortage, the firm cannot produce an output higher than x = 4.

What are the objective and constraint functions?

The objective function: π = -x^2 + 10*x - 6
The constraint: x ≤ 4 or 0 ≤ 4 – x

x;

Maxima will not allow an assignment to be made to the symbol π. So use Pr to
represent profit π,

Pr: -x^2 + 10*x - 6;
soln : solve ( diff (Pr, x), x);

at (Pr, soln);

We can determine if this is a local maximum by looking at the sign of the second derivative.

at (diff (Pr, x, 2), soln);

So the profit Pr will be maximum if the output x = 5, but our constraint is that x <= 4,
so the constraint is binding. The optimum solution, taking into account the constraint, is
x = 4, a solution "on the boundary" of the constraint x <= 4 (as close to the constraint-free 
optimum as possible).

wxdraw2d (xlabel = "x", yrange = [-5, 25], 
     title = "Profit Curve and Binding Constraint", 
    key = "Profit",  explicit (Pr, x, 0, 12), 
    color = black, key = "", explicit (0, x, 0, 12),
     line_width = 1, parametric (5, yy, yy, -5, 19), color = red,
    line_width = 2, key = "constraint", parametric (4, yy, yy, -5, 25) ),
     wxplot_size = [680, 680]$

[at (Pr, x = 4),  at (Pr, x = 5)];

So the given constraint x <= 4 is binding in this problem. A hypothetical  constraint
 x <= 6 would NOT be binding, since the constraint-free optimum (x = 5) is accessible 
without violating the x <= 6 constraint.