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

# Production Theory

A firm's production is governed by a (probably unknowable) function
that relates the quantity produced of a good (or a set of goods) to the
quantities of the various inputs that the firm employs:
        x = f(K, L, ...),
where x is the quantity produced per period, K is the number of 
units of capital services per period, and L is the number of units
of labor per period. The ellipses indicate that other inputs may be
involved in production. This expression treats units of K and L as
homogeneous. We retain this simplifying assumption for the 
remainder of the analysis.

## An Illustrative Function

The function used to illustrate some issues related to production
is this: 
    x = a*K^e * L^(1 - e)  +  b*K^(2*e) * L^(2*(1 - e))  +
          c*K^(3*e) * L^(3*(1 - e)). 
Like the Cobb-Douglas function that it generalizes, this function
is homothetic. It is not, however, homogeneous, as we see below.

f(K,L,a,b,c,e):=a*K^e*L^(1-e)+ 
      b*K^(2*e)*L^(2*(1-e))+c*K^(3*e)*L^(3*(1-e));

Later analysis will depends on the average product of capital, 
the marginal product of capital, and the counterparts for labor.
These expressions appear below.

APK(K,L,a,b,c,e):=''(f(K,L,a,b,c,e)/K);
APL(K,L,a,b,c,e):=''(f(K,L,a,b,c,e)/L);  
MPK(K,L,a,b,c,e):=''(diff(f(K,L,a,b,c,e),K));
MPL(K,L,a,b,c,e):=''(diff(f(K,L,a,b,c,e),L));

Before developing the graphs, we determine the values of
K and L that lead to maximum output, given a set of 
parameters. The parameter values appear below.

The fact that such critical values exist implies that 
returns to scale needs to be considered. We do that before
returning to the behavior of the expressions above.

[a0, b0, c0, e0] : [2, 10, - 0.2, 0.3]$

## Returns to Scale

This section examines the effect of changing the employment
level of both inputs but keeping the K/L value constant. That is,
it examines the effect of changing the scale of operation.

### Ranges of Positive and Negative 
Returns to Scale

We formalize the term returns to scale and determine the employment
levels for which returns to scale are positive. We derive an 
an expression for the production function given that m = K/L, where 
m is a positive constant. Output is expressed in terms of m and L.

The next cell shows the production function when K = L*m and 
the result of solving the first derivative of this function for L.

radcan(f(m*L,L,a0,b0,c0,e0)); 
solnm:solve(diff(f(m*L,L,a0,b0,c0,e0), L)=0, m);

The graph shows the inverse relationship between
m and the value of L at which returns to scale
become negative. To maximize its profits a firm must
operate in the "Positive returns to scale" area.

wxdraw2d( xlabel= "K/L = m",ylabel="L",
 implicit(solnm[1], m, 0.01, 5, L, 0.1, 75),
 label(["Positive returns to scale", 2.5, 10]),
 label(["Negative returns to scale", 3, 50]),
 title="L, m, and returns to scale"     )$

The first command in the next cell names the implicit relationship imp.
The second command substitutes a value of m into imp and solves the
resulting expression for L, yielding the values of L that lead to maximum
x. Only the positive value makes economic sense.

imp: solnm[1];  
solnm0:float(solve(subst(m=1,imp), L));

m = 1 implies that K = L. The resulting maximum value of x
is determined below.

L0:rhs(solnm0[2])$ K0:L0$ x0:f(K0,L0,a0,b0,c0,e0)$
matrix(["L0","K0,","x0"],[L0,K0,x0]);

### Increasing and Decreasing Returns: 
Fixed Proportions

The preceding cells show how to state the production function
in order to determine the maximum output that a firm can 
produce, given a value of m, m = 1.

Here, we extend this analysis to divide the range of positive
returns to scale into the range of increasing returns and the 
range of decreasing returns. 

We use another value of m, m = 0.5. The analysis could have
been conducted in terms of m = 1.

The original function and its first derivative appear below.
The name xm1 is assigned to the original expression for
use in a graph below.

m1:0.5$ xm1:f(m1*L,L,a0,b0,c0,e0);
diff(f(m1*L,L,a0,b0,c0,e0),L);

As with m = 1, we determine the values of K and L that lead
to maximum x, and the value of x. These values appear below.

matrix(["L values", "L","K","x"], 
            [solnimpm1: float(solve(subst(m1,m,imp), L)), 
             L1:rhs(solnimpm1[2]),  
             K1:m1*L1,  
             x1:f(K1,L1,a0,b0,c0,e0)]);

The second derivative with respect to L divides the range
of positive returns to scale into its two component parts.
For L < 20.519, the function exhibits increasing returns
to scale; for L > 20.519, it exhibits decreasing but positive
returns to scale; and for L = 20.519, it exhibits constant
returns to scale.

solve(diff(xm1,L,2),L)$    L2:float(rhs(%[1]))$  
K2:m1*L2$   x2:f(K2,L2,a0,b0,c0,e0)$
matrix(["K","L","x"] ,[K1,L1,x1], [K2, L2, x2]  );

The graph below depicts these results. Other than the
color and dimension changes, this graph differs from 
the one in the text in two respects. First, it assigns the
expression for the production fuction to the key for
the explicit( ) command. Second, it places the value of
m1 in the xlabel. These changes involve the concat( )
and string( ) commands. Note the key="" used to 
turn off the key option.

Exercise: Remove key="" to see the effect.

wxdraw2d(
 xlabel = concat("L, when K = L*",string(m1)),
ylabel="x",  yrange=[0,1.1*x1], key=string(xm1),
 explicit(xm1,L,0.1,50),points_joined=true, 
 line_type=dots, color=black,   key="",
 points([ [L2,0],[L2,x2] ]),points([[L1,0],[L1,x1]]),
 label(["Diminishing Returns", L2+.5*(L1-L2),x1/4]), 
 label( ["Increasing Returns", L2/2, x2/4]),
 dimensions=[540,480])$

The table below shows a set of values consistent with 
K = m1*L. These values are added to the surface plot, 
illustrating how the production function with a constant input ratio 
relates to the production function. 

Llist: makelist(i*L2/5.0,i,0,6)$ /*  <==List of L values */
Klist:m1*Llist$                        /*  <== List of K values*/          
xlist : f(Klist,Llist,a0,b0,c0,e0)$  /* <== List of x values*/
addrow(matrix(
 ["K","L","x"]),transpose(matrix(Klist,Llist,xlist)));

The graph below represents the production function in a 
three-dimension graph. The lists appear  as a set of points. 
The title used in the book is replaced with two key entries.
The zrange is extended slightly in order to move the key
away from the graph.

x1;

wxdraw3d(view=[60,330],zrange=[0,1.1*x1],xlabel="K",
  ylabel="L",zlabel="x",xtics={0,m1*L2},color=green,
  ytics={0,L2},ztics={0,x0}, key="Product Function",
  explicit(f(K,L,a0,b0,c0,e0),  K,.1,1.2*K2,L,.1,1.2*L2), 
  color=blue, key="x when K = m1*L",points(Klist,Llist,
  xlist),dimensions=[480,480] )$

## Production with Variable Proportions

The graph below extends the ones above in two ways. First, it 
extends the range of employment to include values at which 
increasing in one or both of the inputs causes output to decrease.

The graph also shows contour lines that emphasize the fact that
a firm can substitute among inputs. The contour lines shows 
combinations of K and L that generate x = 3390, approximately.
The segment of the line nearest the origin (the southwestern
segment) is the relevant part. Being on the other segment would
generate the same output but at a higher cost.

We examine the implications of substitutions below. First, we 
consider what happens when only a single input's employment
can be varied. 

Suggestion: Copy the cell below and paste it into another location.
Remove the wx from wxdraw3d. Execute the cell and rotate the 
resulting graph to see which view you find most informative.

wxdraw3d(view = [55,330],dimensions=[480,480], 
 contour=base,contour_levels={.9*x0},
 xlabel="K",ylabel="L",zlabel="x",xtics={K0/4,3*K0}, 
 ytics = {L0/4,  1.5*K0}, ztics = {0,x0},color="green", 
 explicit(f(K,L,a0,b0,c0,e0),K,.1,3.0*K0,L,.1,1.5*L0))$

## Varying a Single Input

Fix K = K0, which was defined earlier. Then the production function
becomes a function of L, the only input whose employment
level can be varied:

matrix(["K's fixed value", "The production function "],
  [K0,f(K0,L,a0,b0,c0,e0)]);

The commands below do the following:
* Recall the values of L0 and K0
* Generate an object that shows x as a function of K with L's 
value held constant (Ktotalcurve) anbd as a function of L with
K's value held constant (Ltotalcurve). These are cross-sections
of the 3-dimensional production function over the range of 
positive returns to scale.
* Derive the per-unit curves, average and marginal, from the
two total product curves.
* Draw these curves.
To avoid clutter, the keys are suppressed. Total product
curves are black, average product curves are blue, and 
marginal product curvew are red.

print("L0 =", L0, " and K0 =", K0)$
Ktotalcurve:gr2d(xlabel="K", ylabel="x",yrange=[0, 1.1*x0], 
 ytics=2000,xtics=10,line_width=2,color=black,
 explicit(f(K, L0,a0,b0,c0,e0),K, .1,1.1*K0))$
Kunitcurves:gr2d(yrange=[-50,200],xlabel="K",xtics=10,
  ylabel="x per unit of K", xaxis=true,ytics = 100,
  line_width=2,color=blue,
  explicit(APK(K, L0, a0, b0, c0, e0), K, .1, 1.1*K0),
  color=red,
  explicit(MPK(K, L0, a0, b0, c0, e0), K, .1, 1.1*K0) )$
Ltotalcurve:gr2d(  yrange = [0, 1.1*x0],ytics=2000,
  xtics=10,xlabel = "L", ylabel="x",line_width=2,color=black,
  explicit(f(K0, L, a0, b0, c0, e0), L, .1, 1.1*L0) )$
Lunitcurves: gr2d( xaxis = true, yrange=[-50, 200],
  xlabel="L",ylabel="x per unit of L",xtics=10,ytics=100,
  color=blue,line_width=2,
  explicit(APL(K0, L, a0, b0, c0, e0), L, .1, 1.1*L0),
  color=red,
  explicit(MPL(K0, L, a0, b0, c0, e0), L, .1, 1.1*L0))$
wxdraw(dimensions=[540,480],Ktotalcurve,  Ltotalcurve, 
    Kunitcurves, Lunitcurves, columns =2)$

### Interactions Between Inputs

The next cell creates two total product curves, two average
product curves, and two marginal product curves for input L.
Each curve corresponds to a value of input K.

totals: gr2d(user_preamble="set key bottom right",
 yrange=[0,1.1*x0],title="Total Product of L for 2 Values of K",
 xlabel = "L", ylabel="x", ytics=2000,
 line_width=2,key="x when K = K0/4",
 explicit(f(K0/4,L,a0,b0,c0,e0),L,.1,2*L0),color=red,
 key="x when K = K0/2",explicit(f(K0/2,L,a0,b0,c0,e0),L,.1,2*L0))$
averages: gr2d(   user_preamble="set key bottom left",
  xlabel = "L", ylabel="x per unit of L", ytics=50,
  title="Average Product of L for 2 Values of K",
  line_width=2, key = "APL when K = K0/4",
  explicit(APL(K0/4, L, a0, b0, c0, e0), L, .1, 2*L0),
  key = "APL when K = K0/2", color=red,
  explicit(APL(K0/2, L, a0, b0, c0, e0), L, .1, 2*L0))$
marginals: gr2d(yrange=[-50, 200],xaxis=true,
 user_preamble="set key top right", xlabel = "L", 
 ylabel="x per unit of L", line_width=2, 
 title="Marginal Product of L for 2 Values of K",
  key = "MPL when K = K0/4",
    explicit(MPL(K0/4, L, a0, b0, c0, e0), L, .1, 2*L0),
  key = "MPL when K = K0/2", color=red,
  explicit(MPL(K0/2, L, a0, b0, c0, e0), L, .1, 2*L0) )$
wxdraw(totals, averages,marginals, dimensions = [540,640])$

## Isoquants and Substitution

The contour line shown above is an isoquant, the production
counterpart to an indifference curve.  An isoquant passes through 
any combination of L and K and shows the other combinations 
of L and K that would yield the same quantity. 

As with indifference curves, the slope of the isoquant is 
important. The next cell derives a functional expression of
the slope, which is named the marginal rate of technical 
substituton.

Casual examination of the output shows that the function is
homothetic. A given ratio of K to L produces the same slope,
regardless of the levels of the two inputs.

/*mrts(K,L,e) := ''(
    radcan(-diff(f(K, L,a,b,c,e),L)/diff(f(K,L,a,b,c,e),K)));
The quote-quote operator is replaced with define ( )
*/
define(mrts(K,L,e),
      radcan(-diff(f(K, L,a,b,c,e),L)/diff(f(K,L,a,b,c,e),K)));

define(mrts(K,L,e),
      radcan(-diff(f(K, L,a,b,c,e),L)/diff(f(K,L,a,b,c,e),K)));

### Elasticity of Substitution

As noted previously, this function generalizes the Cobb-Douglas
function. Like that function, this one has unitary elasticity of
substitution between K and L.

[k:mrts*(e/(1-e)), diff(k,mrts)*(mrts/k)];

define(mrts(K,L,e),
      radcan(-diff(f(K, L,a,b,c,e),L)/diff(f(K,L,a,b,c,e),K)));

### Effect of L on mrts

Given that mrts(K,L,e) = [(e-1)/e]*K/L}, 
the mrts must approach zero as L becomes large, as we
move along any isoquant, because increasing L implies
reducing $K$, so the ratio decreases.

If K is constant, then increasing L causes movement to a higher 
isoquant, to a point at which mrts has decreased. The graph below 
compares the slopes at the intersections of the three isoquants 
and the dotted line, where K = 1.5.

### Illustration

The following input group defines the three output levels 
for which isoquants are drawn below. It also determines
the value of mrts when K/L = m1.

[L1,L2,L3]:[2,4,6]$    [K1,K2,K3]:m1*[L1,L2,L3]$  
[x1,x2,x3]:[f(m1*L1,L1,a0,b0,c0,e0),
 f(m1*L2, L2,a0,b0,c0,e0), f(m1*L3,L3,a0,b0,c0,e0)]$ 
slope:mrts(m1*L1,L1,e0)$  matrix(["slope =",slope]); 
matrix(["L","K","x"],[L1,K1,x1],[L2,K2,x2],[L3,K3,x3]);

wxdraw2d(xlabel="L, when K = m1*L",ylabel="K",
 yrange=[0,2*K3], key="Isoquants and tangents",
 implicit(f(K,L,a0,b0,c0,e0)=x1,L,.1,1.5*L3,K,.1,1.5*K3),
 key = "",explicit(K1-L1*slope + slope*L, L, L1-1, 1.5*L1) ,
 implicit(f(K,L,a0,b0,c0,e0)=x2,L,.1,1.5*L3,K,.1,1.5*K3), 
 explicit(K2-L2*slope+slope*L,L,L2-1, 1.25*L2),
 implicit(f(K,L,a0,b0,c0,e0)=x3,L,.1,1.5*L3,K,.1,2*K3),
 explicit(K3-L3*slope + slope*L, L, L3-2, 1.5*L3),color=red,
 line_width=2,key="K = m1*L",explicit(m1*L,L, 0, 1.5*L3),
 color=black, key = "K = 1.5", explicit(1.5,L,0,1.5*L3))$