# First Order Difference Equations
Dowling17.wxmx

load(draw)$  set_draw_defaults(line_width=2, grid = [2,2], point_type = filled_circle,
             background_color = light_gray,  draw_realpart=false)$ 
            fpprintprec:5$   ratprint:false$     kill(all)$

load ("Econ2.mac");

## Preface

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

This update corrects some syntax errors, misspelled words, and also improves the clarity 
of the discussion. We have also used the draw2d function points(xL,yL) for simplicity.

We use Maxima to solve a few of the problems (and draw some of the plots) in Ch. 17 of the 
supplemental text: Introduction to  Mathematical Economics, 3rd ed, (Schaum's Outline 
Series), by Edward T. Dowling (1992, 2001), McGraw-Hill. This modestly priced text is a 
bargain with many worked out examples. You should compare the examples worked out 
"by hand" in this text with what we do using Maxima. 

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

Ted Woollett
https://home.csulb.edu/~woollett/
woollett@charter.net
Dec. 13, 2021

## Definitions and Concepts

A difference equation expresses a relationship between a dependent variable and a lagged 
independent variable (or variables) which changes at discrete intervals of time, for example, 
    I[t] = F( Y[t - 1] )  
with F(x) being some arbitrary function, where I and Y are each measured at the end of 
each year. Investment during year 4 is some function of national income during year 3, for
example: I[4] = F (Y[3]), etc.

The order of a difference equation is determined by the greatest number of periods 
lagged. 

A first-order difference equation contains only a time lag of one period, as above.
 
A second-order difference equation involves a time lag of two periods, in addition to
 (perhaps) also a time lag of one period, etc. 

The change in y as t changes from t to t +1  is called the "first difference of y". 

The FIRST DIFFERENCE of y is written :
  Δy/Δt  =  Δy[t] = y[t + 1] - y[t]  = change in y per period over the time (t, t+1).

In Maxima, y[t] is a (hash) array and wxMaxima will display the symbol using subscript
notation ( if display2d is set to true, which is the default).

y[t];

display2d;

### Example 1

Each of the following is a difference equation of the order indicated: 

Order 2:

I[t] = a* ( Y[t-1] - Y[t-2]);

Order 1:

Q[t] = a + b*P[t-1];

Order 3:

y[t+3] - 9*y[t+2] + 2*y[t+1] + 6*y[t] = 8;

Order 1:

Δy[t] = 5*y[t];

The last equation is the same as:

y[t+1] - y[t] = 5*y[t];

which is the same as:

y[t+1] = 6*y[t];

### First Few Values, history1 (y0, fnext), find_root(expr,x,a,b)

The first few values of a first order difference equation can be found starting from an
initial value, y0 = y[t=0]. 

For example, suppose y[t] = 0.5*y[t-1]^(-1/3), a nonlinear difference equation. We take
y0 = 0.75 as the starting value.

kill(y0,y1,y2,y3,y4)$
f(x) := 0.5*x^(-1/3);
[y0 : 0.75, y1 : f(y0), y2 : f(y1), y3 : f(y2), y4 : f(y3)];

history1 (_y0, fnext) := ([y0 : _y0, y1 : fnext(y0), y2 : fnext(y1), y3 : fnext (y2), 
            y4 : fnext (y3)] )$

For fnext, we use our definition of the Maxima function f(x) just above, and use the
same starting value y0.

history1 (0.75, f);

Decreasing y0, but using the same fnext,

history1 (0.55, f);

We anticipate the discussion of phase diagram analysis for a first order difference
equation (in the last section), and proceed to make a plot of yt vs ytm1, in which
yt stands for y[t], and ytm1 stands for y[t-1]. If there exists a stable value of y[t], it should
have the property y[t] = y[t-1] and hence lie on the 45 degree line in the first quadrant,
which is drawn in the color red below.

For our example, which uses the definition of the Maxima function f(x) above,

yt : f(ytm1);

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]", key_pos = top_left,
    yrange = [0, 3],
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 1),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 1))$

The Maxima function find_root (expr, x, a, b) looks for a numerical root of expr = 0 as x
varies in the closed interval [a,b]. 

find_root (A(x) = B(x), x, a, b) looks for the root of A(x) - B(x) = 0 for x in [a,b].

Let ystar be the root of the equation y = f(y), in which the variable is y, and we are again
using our definition of the Maxima function f(x) above.

ystar : find_root (y = f(y), y, 0.4, 0.8);

The first derivative of the slope of the curve f(y), ie., the curve of yt vs ytm1, evaluated at
a point where y = f(y), tells us whether the point ystar is a stable equilibrium for which the
given difference equation tends to drive the next value of y (for the next value of t).
The "time" t is in the background of this analysis, and in the last section we illustrate the
approach to a stable point like ystar, using a (cobweb) plot. 

We have defined yt above as a specific Maxima expression depending on ytm1:

yt;

The local slope of the curve yt(ytm1) is the first derivative of yt wrt ytm1.

slope : diff (yt, ytm1);

For ytm1 > 0, the slope is thus always a negative number.
What is the local slope value for the "equilibrium point" ytm1 = ystar?
If the MAGNITUDE of the slope at ystar is less than 1, then the given difference equation
will drive yt to the equilibrium point.

ystar;

at (slope, ytm1 = ystar);

The magnitude of the slope of the yt vs ytm1 curve at the value y = ystar is less than 1,
so a solution beginning with some y0 value will converge to y = ystar. The slope at ystar
is a negative number, so successive values of y will oscillate around y = ystar in the
process of convergence.

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]",  yrange = [0, 3],
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 1), color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 1), color = black, key = "", line_width = 1,
   explicit (ystar, ytm1, 0, ystar),  parametric (ystar, yy, yy, 0, ystar))$

### Solution of a Difference Equation, Ex. 2

The "solution" of a difference equation defines y[t] for every discrete value of t = 0,1,2,... 
and does not contain a difference expression.

Given that the initial value of y is y0 = y[0], in the example difference equation
         y[t+1] = b*y[t],
a solution is found as follows. By successive substitutions of t = 0, 1, 2, 3, etc, 
we get
   y[1] = b*y[0], 
   y[2] = b*y[1] = b*( b*y[0]) = b^2*y[0], 
   y[3] = b*y[2] = b^3*y[0],
   .....
   y[t] = b^t*y[0],  for t = 0, 1, 2, 3,....

This method is called the "iterative method." Since y[0] is constant, note the
crucial role b plays in determining values for y[t] as t increases.

## General Solution for y[t] = b*y[t-1] + a

Given a first order difference equation which is LINEAR, (i.e., all the variables are
raised to the first power only and there are no cross products) 

      y[t] = b*y[t-1]  +  a,

the general formula for a "definite solution", with y0 equal to y[0], is
     if b # 1  then      y[t] = (y0 - a/(1-b))*b^t  + a/(1-b), 
     otherwise, for b = 1, y[t] = y0 + a*t.

If no initial condition is given, the formula for a general solution can be given
in terms of some (as yet undefined) constant A as:

  if b # 1 then      y[t] = A*b^t + a/(1-b),  
  otherwise, for b = 1,  y[t] = A + a*t.

Note that b^0 = 1 for any b.

### Maxima Derivation of the General Solution

The Maxima function solve_rec (eqn, var, optional-initial-assignments)
can be used to solve for hypergeometrical solutions of a linear recurrence 
equation (eqn) with polynomial coefficients in the variable var. 

Optional are one or more "initial conditions" of the form y[0] = 1,y[3] = 2, etc, if var = y[t].

To use this function first use load ("solve_rec");  which is a "contributed package" designed
to solve recurrence relations.

load ("solve_rec");

The syntax is 
       solve_rec ( eqn, var, [init] ),
in which the brackets around init mean that one or more initial condition statements
of the form y[0] = y0, y[2] = -5, etc are optional; only eqn and var are required arguments.

kill (y0)$
soln : solve_rec (y[t] = b*y[t-1] + a, y[t], y[0] = y0);

The Maxima function grind prints out the expression in display2d = false style:

grind(soln)$

rhs(soln);

The coefficient of b^t can be picked out as follows

coeff(%, b^t);

Replacing a/(b-1) with - a/(1-b) throughout, we get the solution given above, for the case
that b # 1.

If we don't supply the optional initial condition, solve_rec produces

soln : solve_rec (y[t] = b*y[t-1] + a, y[t]);

The solve_rec algorithm returns a solution involving an arbitrary constant %k. If we define
another arbitrary constant A by:
   A = a/(b-1) + %k, 
which is the coefficiant of b^t, we get y[t] = A*b^t + a/(1-b), as given above for b # 1.

Note that solve_rec did not consider separately the special case b = 1.
By setting b = 1, we get that case as a solution from solve_rec.

solve_rec (y[t] = y[t-1] + a, y[t], y[0] = y0);

which agrees with our b = 1 solution above.

### Function DE1(a,b,t,y0) or DE1(a,b,t)

The software file Econ2.mac contains the definition of the Maxima function DE1.

DE1 (a, b, t, y0) returns the value of y[t] governed by the first order
difference equation (DE)  y[t] = b*y[t-1] + a, with y[0] = y0.

DE1 (a, b, t ) returns the general solution for y[t] governed by the first order
difference equation (DE):  y[t] = b*y[t-1] + a , in terms of an arbitrary constant
%A.

### Example 3

Consider the difference equation
  y[t] = -7*y[t-1] + 16, with y[0] = 5
in which a = 16, the base b = -7, and y0 = 5.

Let de be a Maxima expression returned by DE1.

de : DE1 (16, -7, t, 5);

To define an ordinary Maxima function y(t) at this point, we need to precede the expression
de with TWO single quotes (') to force evaluation.

y(t) := ''de;

grind(y(t))$

Check this definition at t = 0 and t = 1.

map ('y , [0,1]);

which is equivalent to

[y(0), y(1)];

y[1]'s value should be given in terms of this ordinary Maxima function y(t) using
y(1) = -7*y(0) + 16, using the given difference equation which must be satisfied for
discrete values of t = 0, 1, 2, ...

-7*5 + 16;

which is the correct value of y[1].

Without supplying the fourth arg y0 when calling DE1, we get:

de : DE1 (16, -7, t);

which involves the arbitrary constant %A.

## Stability Conditions: General Characterization

The exponential expression b^t will generate seven qualitatively different time paths 
depending on the nature of the base b. The general solution
         y[t] = A*b^t + c
will be dynamically stable only if the "complementary function" A*b^t approaches 0 as
t approaches infinity. This then all depends on the general nature of b, as illustrated in 
the following  plots.

If the magnitude of b is greater than 1 ( |b| > 1), the time path will explode and move
farther and farther away from equilibrium. If the magnitude of b is less than 1 (|b| < 1), 
the time path will be damped and will move toward equilibrium.

If b is negative (b < 0), the time path will oscillate between positive and negative values. 
If b is positive  (b > 0),  the time path will be non-oscillating.

Here we can make some plots for different values of b, in the simple case y[t] = b^t,
implied by the  difference equation y[t] = b*y[t-1], with y[0] = 1.

We use the draw2d function points (xL, yL) to plot the evolution in time.

tL : makelist (t, t, 0, 6);
yL : makelist (1.5^t, t, 0, 6);

Note that at the top of this file we defined a default point type as a filled circle for
draw2d.

wxdraw2d(yrange = [0,12], xrange = [0, 7], xlabel = "t",
    ylabel = "y", key_pos = top_left, key = "y = 1.5^t", 
    points (tL, yL))$

The jump in the value of b^t gets larger each successive discrete value of t we use. This 
is an example of an "explosive instability."

As a constrast, we show convergence to zero with a value of b in the range 0 < b < 1, 
say b = 0.5, as an example. For t = 0, 0.5^0  = 1.

yL : makelist (0.5^t, t, 0, 6);

wxdraw2d(yrange = [-0.5, 1.5], xrange = [0, 7], xlabel = "t",
    ylabel = "y",  key = "y = 0.5^t", points (tL, yL), key = "", color = black,
    line_width = 1, explicit (0,t,0,7))$

For -1 < b < 0 we get oscillating convergence of b^t.
Here is an example with b = - 0.5.

yL : makelist ((-0.5)^t, t, 0, 6);

wxdraw2d(yrange = [-1, 1.5], xrange = [0, 7], xlabel = "t",
    ylabel = "y(", key = "y = (-0.5)^t", points (tL, yL), key = "",
    line_width = 1, color = black, explicit (0, t, 0, 7))$

With b < -1 we get oscillating explosive instability.
Here is an example with b = - 1.5.

yL : makelist ( (-1.5)^t, t, 0, 6);

wxdraw2d(yrange = [-10, 12], xrange = [0, 7], xlabel = "t",
    ylabel = "y", key_pos = top_left, key = "y = (- 1.5)^t", 
    points (tL, yL), key = "",  line_width = 1, color = black, explicit (0, t, 0, 7))$

The borderline case b = -1 shows oscillating behavior with no convergence and no
growth trend.

yL : makelist ( (-1)^t, t, 0, 6);

wxdraw2d(yrange = [-2, 2], xrange = [0, 7], xlabel = "t",
    ylabel = "y",  key_pos = top_left, key = "y = (-1)^t", 
    points (tL, yL), key = "",line_width = 1, color = black, explicit (0, t, 0, 7))$

The other borderline case is b = 1, for which b^t is a flat line at y = 1.

yL : makelist (1^t, t, 0, 6);

wxdraw2d(yrange = [-2, 2], xrange = [0, 7], xlabel = "t",
    ylabel = "y",  key_pos = top_left, key = "y = 1^t", 
    points (tL, yL), key = "",
    line_width = 1, color = black, explicit (0, t, 0, 7))$

## Lagged Income Determination Model

In the simple income determination model of Dowling Ch. 2, Sec. 2.3, there were no lags. 
Now assume  that  consumption is a function of the previous period’s income, and national 
income Y is the sum of consumption plus investment, and that investment per year is 
constant.

/* these are not assignment statements!  */
C[t] = C[0] + c*Y[t-1];
Y[t] = C[t] + I[t];
I[t] = I[0];

Writing Y[t] in standard form,  Y[t] = c*Y[t-1] + C[0] + I[0], and letting C0 stand for C[0], 
I0 stand for I[0], by comparison with our standard form: y[t] = b*y[t-1] + a, we identify b = c, 
a = C0 + I0. Let Y0 stand for Y[0], the optional initial value arg in our Maxima function
DE1 defined in Econ2.mac.

DE1 (C0 + I0, c, t, Y0);

With the "propensity to consume" c > 0 and c < 1, this solution converges to (I0 + C0)/(1-c)
with no oscillation, based on our previous examples.

As a numerical example, assume C0 = 200, I0 = 100, c = 0.9, Y0 = 4500.

de : DE1 (300, 0.9, t, 4500);

Use two single quotes (') in front of the expression de to allow us to define a Maxima
function Y(t) which can be mapped onto a list of times.

Y(t) := ''de;

tL : makelist (t, t, 0, 50)$
YL : map ('Y, tL)$

wxdraw2d(yrange = [2900, 4500], xrange = [0, 50], xlabel = "t",
    ylabel = "Y[t]", key = sconcat ("Y[t] = ", Y(t)), points (tL, YL),
    color = black, line_width = 1, key = "",  explicit (3000, t, 0, 50))$

## The Cobweb Model

For many products, such as agricultural commodities, which are planted a year before 
marketing,  current supply depends on last year’s prices. This poses interesting stability 
questions. Let Qd[t] be the quantity demanded function in terms of prices and Qs[t] be
the quantity supplied function depending on prices the previous year. The equilibrium
price during year t (t = 0, 1, 2, ...) is then found by equating Qd[t] and Qs[t].

Qd[t] : c + b*P[t];
Qs[t] : g + h*P[t-1];
solns : solve (Qd[t] = Qs[t], P[t]);

solns : expand (solns);

Use DE1(a,b,r,y0)  to find the solution for P[t], assuming P0 is the initial value P[0].
Here "a" = (g-c)/b, "b" = h/b.

DE1 ((g-c)/b, h/b, t, P0);

Evidently the equilibrium price is Pe = (g-c)/(b-h), which appears inside the first term as
well as the last term. Let's show an equation with  P[t] expressed in terms of Pe.

P[t] =  (P0 - Pe)*(h/b)^t + Pe;

de : rhs(%);

Now bind the symbol Pe to its symbolic value:

Pe : (g - c)/(b - h);

To get de fully evaluated, incorporating the above assignment of the symbol Pe,
we use two single quotes (') in front of de.

de : ''de;

If |h| > |b|, then |h/b| > 1, and the time path P[t] explodes.

If |h| = |b|, then |h/b| = 1 and the time path P[t] oscillates without convergence.

With an ordinary "negative demand function" (b < 0) and "positive supply function" (h > 0),
h/b < 0, and the time path will oscillate. 

If |h| < |b|, then we get convergence to Pe in this ordinary case since |h/b| < 1.

Quoting Dowling here:
"In short, when Q = f(P) in supply-and-demand analysis, as is common in mathematics, 
the supply curve must be flatter than the demand curve for stability. See Example 7 and 
Problems 17.21 to 17.25.  But if P = f(Q), as is typical in economics, the reverse is true;
The demand curve must be flatter (or more elastic), than the supply curve if the model is 
to be stable."

As a numerical example, (we declare t an integer here after finding this version of Maxima
has some trouble with combining factors each raised to the power t.

declare (t, integer)$
deEx : subst ([c = 86, b = -0.8, g = -10, h = 0.2], de);

1.25*0.2;

Maxima is having trouble combining factors with the same exponent (t) here.

Since -1.25*0.2 = -0.25, our result is
P[t] = (P0 - 96)*(-0.25)^t + 96

With the base "b" equal to -0.25, we have oscillatory convergence to the equilibrium
price Pe = 96 in this numerical example.

### Visualizing the Cobweb

We first ignore time lags and draw Q(P) curves (lines) for both demand (D) and
supply (S), in which the supply curve is "flatter" than the demand curve.

Qd : 86 - 0.8*P;
Qs : -10 + 0.2*P;

The coordinates of the supply-demand crossing are (Pe,Qe) in this Q(P) plot.

Pe : at (P, solve(Qd = Qs, P));

Qe : at (Qs, P = Pe);

Here we are following the mathematics convention with Q (the dependent variable) on the 
vertical axis and P (the independent variable) on the horizontal axis. Note that our model 
parameters imply that Qs = 0 when P = 50, meaning that suppliers need a price of at 
least 50 to supply any of the good being modeled.

wxdraw2d(xlabel = "P", ylabel = "Q", yrange = [0, 20], xrange = [40, 140],
   key_pos = top_left, key = "D",explicit (Qd, P, 40, 140), color = red,
    key = "S", explicit (Qs, P, 40,140))$

Suppose the initial ((t = 0) market price P0 = 120 and let the supply curve predict what the
quantity supplied will be in year t = 1. Let Q1 stand for Qs[1].

P0 : 120$

Q1 : at (Qs, P = P0);

Add a vector to the plot, from point (120,0) to point (120,14), recalling the vector syntax
used by draw2d:   vector ([x,y], [dx, dy]). For this task, x = 120, y = 0, (the tail of the vector),
and dx = 0,  dy = 14. (dx,dy) descibe the orientation and length of the vector, whereas
(x,y) describe the tail of the vector. Thus we use vector ([P0, 0], [0, Q1]).

wxdraw2d(xlabel = "P", ylabel = "Q", yrange = [0, 20], xrange = [40, 140],
   key_pos = top_left, key = "D",explicit (Qd, P, 40, 140), color = red,
    key = "S", explicit (Qs, P, 40,140),key = "", color = black, line_width = 1,
    vector ([P0, 0], [0, Q1]))$

Let the demand curve determine the price P1 offered for quantity Q1 = 14

soln : solve (Qd = 14, P);

P1 : at (P, soln);

Add a vector from (P0,Q1) to (P1, Q1), ie., from (120,14) to (90,14).

wxdraw2d(xlabel = "P", ylabel = "Q", yrange = [0, 20], xrange = [40, 140],
   key_pos = top_left, key = "D",explicit (Qd, P, 40, 140), color = red,
    key = "S", explicit (Qs, P, 40,140),key = "", color = black, line_width = 1,
    vector ([P0, 0], [0, Q1]), vector ([P0, Q1], [P1 - P0, 0]))$

The quantity supplied for year 2, Q2, is then determined by the supply curve and the
price P1 = 90 offered in year 1.

Q2 : at (Qs, P = P1);

Add a vector from (P1, Q1) = (90,14) to (P1, Q2) = (90, 8).

wxdraw2d(xlabel = "P", ylabel = "Q", yrange = [0, 20], xrange = [40, 140],
   key_pos = top_left, key = "D",explicit (Qd, P, 40, 140), color = red,
    key = "S", explicit (Qs, P, 40,140),key = "", color = black, line_width = 1,
    vector ([P0, 0], [0, Q1]), vector ([P0, Q1], [P1 - P0, 0]),
    vector ([P1,Q1], [0, Q2 - Q1]))$

Let the demand curve determine the price P2 offered for quantity Q2 = 8.

soln : solve (Qd = Q2, P);

P2 : at (P, soln), numer;

Add a vector from (P1,Q2) = (90,8) to (P2,Q2) = (97.5, 8).

wxdraw2d(xlabel = "P", ylabel = "Q", yrange = [0, 20], xrange = [40, 140],
   key_pos = top_left, key = "D",explicit (Qd, P, 40, 140), color = red,
    key = "S", explicit (Qs, P, 40,140),key = "", color = black, line_width = 1,
    vector ([P0, 0], [0, Q1]), vector ([P0, Q1], [P1 - P0, 0]),
    vector ([P1,Q1], [0, Q2 - Q1]), vector ([P1,Q2], [P2 - P1, 0]))$

The quantity supplied for year 3, Q3, is then determined by the supply curve and the
price P2 = 97.5 offered in year 2.

Q3 : at (Qs, P = P2);

Add a vector from (P2,Q2) = (97.5,8) to (P2,Q3) = (97.5, 9.5).

wxdraw2d(xlabel = "P", ylabel = "Q", yrange = [0, 20], xrange = [40, 140],    
   key_pos = top_left, key = "D",explicit (Qd, P, 40, 140), color = red,
    key = "S", explicit (Qs, P, 40,140),key = "", color = black, line_width = 1,
    head_length = 1.2,   vector ([P0, 0], [0, Q1]), vector ([P0, Q1], [P1 - P0, 0]),
    vector ([P1,Q1], [0, Q2 - Q1]), vector ([P1,Q2], [P2 - P1, 0]),
    vector ([P2, Q2], [0, Q3 - Q2])), wxplot_size = [680,680];

A cobweb pattern is gradually drawn as the vectors converge on the equilibrium value of
market price (96).

## The Harrod Model

The Harrod model attempts to explain the dynamics of growth in the economy. It assumes
 national saving are some fraction of the national income:
      S[t] = s Y[t],  
where s is a constant equal to both the marginal propensity to save (MPS) and the average 
propensity to save (APS). With both S and Y expressed in dollars, s is dimensionless.

The Harrod model also assumes the "acceleration principle", i.e.,  investment is proportional 
to the rate of change of national income over time, in the sense:
       I[t] = a (Y[t] - Y[t-1]),
where a is a constant equal to both the marginal and average capital-output ratios.
With both I and Y expressed in dollars, a is dimensionless.

In equilibrium, national investment equals national savings:
       I[t] = S[t],
and hence
   a (Y[t] - Y[t-1]) = s Y[t], or
   (a - s) Y[t] = a Y[t-1].

a is normally larger than 1 and s is positive but less than 1, so (a-s) # 0, and we can 
normally divide by (a - s) to get
      Y[t] = (a/(a-s)) Y[t-1].
Comparing with our standard form y[t] = b*y[t-1] + a, the base b = a/(a-s) and there is no
additional constant term. Find Y[t] using DE1, using Y0 to stand for  Y[0].

DE1(0, a/(a - s), t, Y0);

grind(%)$

So our solution is Y[t] = (a/(a-s))^t*Y0.
The base b of the exponent t is b = a/(a-s), in which 0 < s < 1 and a > 1 normally, so
a/(a-s) is positive but greater than 1, so we have non-oscillatory explosive behavior of
national income. Income will expand indefinitely with no bounds.

### Example 8: Warranted Rate of Growth

Continuing with the Harrod model, the "rate of growth between successive periods"  is
a pure number defined as
       G = (Y[1] - Y[0]) / Y[0]

Y[t] := (a/(a - s))^t*Y0;
G : (Y[1] - Y0)/Y0, ratsimp;

We can absorb the minus sign out front of G with the following trick:

G : (-num(G))/(-denom(G));

The "warranted rate of growth" Gw is the rate of growth which is needed to have an 
equilibrium between savings and investment each year. Since we found Y[t] and then the 
rate of growth G assuming that equilibrium,
    Gw = s/(a-s).

If s = 0.12, for example, and a = 2.12, the "warranted rate of growth is"

at (G, [s = 0.12, a = 2.12]);

ie., 6 per cent.

## Other Economic Applications

### 17.28: Compound Interest Formula from a Difference Equation

Derive the formula for the value P[t] of an initial amount of money P0 deposited at i 
interest (decimal interest rate) for t years when compounded annually. 

When interest is compounded annually, the principal P increases according to
   P[t+1] = P[t] + i P[t] = (1 + i) P[t],
which is mathematically equivalent to (let t --> t-1 everywhere)
   P[t] = (1 + i) P[t-1],
a form which encourages us to use our difference equation function DE1(a,b,t,y0), 
in a case b = (1+i) # 1 since i > 0, and a = 0, and we let P0 stand for P[0].

kill(P0)$
de : DE1(0, 1+i, t, P0);

Value of principal after t years is then 
   P[t] = P0 (1 + i)^t.
We used this in Dowling Ch. 8, Sec. 6.

### Prob 17.29: Market Price vs Prev. Yr. Level of Inventory

Assume
     Qd[t] = c + z P[t],
     Qs[t] = g + h P[t].
Also assume that the market price is determined not only by the usual market clearing
mechanisms (the pursuit of Qs = Qd) but also by the level of inventory (Qs - Qd) available 
the previous year:
        P[t+1] = P[t] - a ( Qs[t] - Qd[t] ), for t = 0, 1, 2, ....

We also make the reasonable assumption that an increase in inventory will tend to reduce
the price, and likewise, a decrease in inventory will tend to increase the market price.
So we assume a > 0.

Let I[t] = Qs[t] - Qd[t], then P[t+1] = P[t] - a I[t]. A larger I[t] value produces a smaller P[t+1]
provided a > 0.

Substituting for Qs and Qd, solving for P[t+1] in terms of P[t], and then replacing t+1 --> t
and t --> (t-1), we get
  P[t] = B P[t-1] + A,
with B = 1 + a (z - h),  A = a (c - g).

Let P0 stand for P[0], and find P[t] using DE1.

de : DE1 (A, B, t, P0);

Evidently, Pe = A/(1 - B.

kill(Pe)$
de1 : ratsubst (Pe, A/(1-B), de);

coeff(de1, B^t);

de2 : %*B^t + Pe;

We given this equation a name for later convenience (Prob. 17.30).

soln : P[t] = de2;

This is our expression for P[t]. Stability depends on the value of B = 1 + a (z - h), and
 hence on the values of a, z, and h.
We assume a > 0 (as above), and under normal conditions, z < 0, h > 0, so a(z-h) < 0, 
  which in turn implies  B < 1.

With B < 1, we can distinguish five cases:

1.)    0 < B < 1,  P[t] converges to Pe = A/(1-B) and is nonoscillatory; this case requires
               |a*( z - h)| < 1.
2.)    B = 0, 0^t = 0, P[t] = Pe = constant, and evidently Pe must equal P0; this case
            requires a*(z - h) = -1.
3.)    -1 < B < 0, P[t] converges to Pe = A/(1-B) and is oscillatory; this case requires
                -2 < a*(z-h) < -1.
4.)     B = -1, uniform oscillation between P0 and (2*Pe - P0); this case requires a*(z-h) = -2.
5.)     B < -1, oscillatory explosion in value; this case has a*(z-h) < -2.

### Prob. 17.30: Numerical Example of 17.29

As a numerical example, take (in 17.29 above) c = 120, z = -0.5, g = -30, h = 0.3,  a = 0.2,
P0 = 200.

Let soln1 be soln (above) in a form in which the explicit dependence on all the model
parameters is visible. We first replace Pe in terms of A and B. Next we replace A in 
terms of a,c,g. Finally, we replace B in terms of a,z,h.

soln1 : subst( [Pe = A/(1 - B), A = a*(c - g), B = 1 + a*(z - h)], soln);

case1730 is a list of replacement values for all the model constants in our example
(Prob. 17.30).

case1730 : [c = 120, z = -0.5, g = -30, h = 0.3, a = 0.2, P0 = 200];

We then find P[t] given these model parameter assumptions:

soln1730 : at (soln1, case1730);

grind(%)$

at (soln1730, t = 0);

With the base B = 0.84, 0 < B < 1, P[t] converges to Pe = 187.5 without oscillation.

rhs(soln1730);

We need to precede % with two single quotes (') to force evaluation to define
a Maxima function f(t) with the right hand side of soln1730.

f(t) := ''%;

We now have a defined Maxima function we can use like any other function in Maxima.
We create a list of times t (years) using makelist and a list of prices P[t], calling the first
list tL and the second list PL.

tL : makelist (t, t, 0, 20);
PL : map ('f, tL);

wxdraw2d (xlabel = "t", ylabel = "P[t]", yrange = [185, 205], xrange =  [0,20],    
    points ( tL, PL ), color = black, line_width = 1, key = "Pe = 187.5",
    explicit (187.5, t, 0, 20))$

## Phase Diagram Analysis for 1st Order Difference Eqns.

Potentially stable equilibrium points in a plot of y[t] vs. y[t-1] are points on the 45 degree line
y[t] = y[t-1]. If the local "slope" of the curve of y[t] has a magnitude less than 1 (at a 
crossing point), that point is a point of stable equilibrium.

### Example 10

Make a phase plot of the nonlinear difference equation:

    y[t] = (y[t-1])^(1/2), which is equivalent to the difference equation
    y[t] = y[t-1]^(0.5).

Potential equilibrium points must satisfy y[t] = y[t-1].
If, in the given de,  we set both y[t] and y[t-1] equal to y,
  y = y^(1/2), or
  y - y^(1/2) = 0.
This is satisfied for y = 0 and y = 1.   

solve cannot deal with nonlinear equations.

solve (y  = y^(1/2), y);

A plot of f(y) vs. y shows where f(y) = 0.

wxdraw2d( xlabel = "y", explicit (y - y^(1/2), y, -1,2), 
    color = black, line_width = 1, explicit (0,y,-1,2))$

find_root can deal numerically with nonlinear equations.

find_root (y = sqrt(y), y, 0.1, 2);

However, if we try using findroot (y = sqrt(y), y, -0.1,0.2) we get an error since
if y < 0, sqrt(y) is a complex number.

We can find the limit  of the expression y - sqrt(y) as we approach y = 0 from values
of y > 0.

limit (y - sqrt(y), y, 0, plus);

To draw a phase diagram for a first order difference equation, it is convenient to let the 
symbol yt stand for y[t] and the symbol ytm1 stand for y[t-1].

In terms of these symbols, our difference equation is yt = ytm1^0.5 = sqrt (ytm1).
And slope = ∂(yt)/∂(ytm1) =  diff (yt, ytm1).

We first draw the yt vs. ytm1 curve, and also the 45 degree line along which yt = ytm1 (the
locus of possible equilibrium points).

The first line in the following cell assigns the symbol yt to the expression sqrt(ytm1).
If we use the optional draw2d setting: proportional_axes = 'xy, a 45 deg line will look like 
a 45 deg line.

yt : sqrt(ytm1);
wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]", key_pos = top_left,
    proportional_axes = 'xy,
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2.5),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2.5), color = black, key = "", line_width = 1,
    parametric (1,yy,yy,0,1))$

Let ytp (yt-prime) stand for the local slope, as a function of ytm1:

ytp : diff (yt, ytm1);

The numerical value of the slope of the phase curve at the point ytm1 = 1 is:

slope : at (ytp, ytm1 = 1);

The slope at ytm1 = 1 is 1/2 > 0, so no oscillation; the slope magnitude < 1, so y = 1 is
a point of stable equilibrium.

For a graphical check of stability, let ytm1 = 0.25 be yt[0]. We can then find succeeding 
values of yt from the governing model difference equation: yt = sqrt(ytm1). We can show 
the evolution of succeeding values of y[t], for t = 1,2,3,..., on our phase diagram.

Let yj = y[j] in this sidebar.

(1): starting at the point (ytm1 = y0 = y[0], yt = 0), draw a vertical line to the intersection 
with the yt vs ytm1 curve (which determines y1):     (y0, 0) --> (y0, y1),

(2): draw the horizontal line (y0,y1) --> (y1,y1) (ends on the 45 deg line), 

(3): draw a vertical line to intersect the yt vs. ytm1 curve (which determines y2):
        (y1, y1) --> (y1, y2),

(4): draw the horizontal line (y1, y2) --> (y2, y2) (ends on the 45 deg line), 

(5): draw the vertical line to intersect the yt vs. ytm1 curve (which determines y3):
        (y2, y2) --> (y2, y3),

and so on.

In the following cell, we use the given model difference equation to define the history
of y[t]; the values of y0, y1,...,y4, assuming y0 = 0.25.

kill (y0,y1,y2,y3,y4)$
[y0 : 0.25, y1 : sqrt (y0), y2 : sqrt(y1), y3 : sqrt(y2), y4 : sqrt(y3)];

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]", key_pos = top_left,
    proportional_axes = 'xy,
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2.5),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2.5), color = black, key = "", line_width = 1,
    parametric (1,yy,yy,0,1), head_length = 0.03, vector ([y0, 0], [0, y1]),
    vector ([y0,y1], [y1-y0, 0]), vector ([y1,y1], [0, y2-y1]), vector ([y1,y2],[y2-y1,0]),
    vector ([y2,y2], [0, y3-y2]), vector ([y2,y3], [y3-y2,0]))$

An initial value  y0 = 0.25 heads closer and closer to y = 1 as the value of t increases.

If the initial value starts out at, say y0 = 2, we can again follow the progression as
t =  1, 2, 3, ...

Here we redefine the successive y[t] values: y0, y1, y2,...
We can then reuse the above graphics instructions without change to show stability
from the high side.

kill (y0,y1,y2,y3,y4)$
[y0 : 2, y1 : sqrt (y0), y2 : sqrt(y1), y3 : sqrt(y2), y4 : sqrt(y3)], numer;

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]", key_pos = top_left,
    proportional_axes = 'xy,
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2.5),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2.5), color = black, key = "", line_width = 1,
    parametric (1,yy,yy,0,1), head_length = 0.03, vector ([y0, 0], [0, y1]),
    vector ([y0,y1], [y1-y0, 0]), vector ([y1,y1], [0, y2-y1]), vector ([y1,y2],[y2-y1,0]),
    vector ([y2,y2], [0, y3-y2]), vector ([y2,y3], [y3-y2,0]))$

We again see convergence to the equilibrium value y = 1, this time from the high side.

### Prob. 17.31

Construct a phase diagram for the nonlinear difference equation
   y[t] = y[t-1]^3
and check for stability.

Let's start with the phase diagram plot.

kill(yt, ytm1)$
yt : ytm1^3;
wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]", key_pos = top_left,
    proportional_axes = 'xy,  yrange = [0, 2],
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2), color = black, key = "", line_width = 1,
   explicit (1,ytm1,0,1), parametric (1,yy,yy,0,1)), wxplot_size = [680, 680];

We see points of possible stable equilibrium: y = 0 and y = 1. 
The slope at (0,0) is 0 and the slope at (1,1) is 3 (right?). So we infer that (1,1) is not a
stable equilibrium point, and (0,0) is a stable equilibrium point.

Visual confirmation with the plot: start with y0 = 0.75, and watch what happens.

kill(y0,y1,y2,y3,y4,f)$
f(x) := x^3;
[y0 : 0.75, y1 : f(y0), y2 : f(y1), y3 : f(y2), y4 : f(y3)];

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]",
    yrange = [0,2], key_pos = top_left,
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2), color = black, key = "", line_width = 1,
    explicit (1, ytm1,0,1), parametric (1,yy,yy,0,1), head_length = 0.03, 
    vector ([y0, 0], [0, y1]), vector ([y0,y1], [y1-y0, 0]), vector ([y1,y1], [0, y2-y1]), 
    vector ([y1,y2],[y2-y1,0]),vector ([y2,y2], [0, y3-y2]), vector ([y2,y3], [y3-y2,0]))$

The phase curve takes y[t] away from y[0] = 0.75 in decreasing y direction towards y = 0.
Next we start at y0 = 1.1, on the high side of y = 1.

[y0 : 1.1, y1 : f(y0), y2 : f(y1), y3 : f(y2), y4 : f(y3)];

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]",
    yrange = [0,3], xrange = [0, 2.5], key_pos = top_left,
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2.5),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2.5), color = black, key = "", line_width = 1,
    parametric (1,yy,yy,0,1), head_length = 0.03, vector ([y0, 0], [0, y1]),
    vector ([y0,y1], [y1-y0, 0]), vector ([y1,y1], [0, y2-y1]), vector ([y1,y2],[y2-y1,0]),
    vector ([y2,y2], [0, y3-y2]), vector ([y2,y3], [y3-y2,0]))$

Again the phase diagram sends y away from y[0] = 1, so the only stable equilibrium point
is y = 0.

### Prob. 17.32

Construct a phase diagram for the nonlinear difference equation
      y[t] = y[t-1]^(-1/4)
and check for stability.

The only possible equilibrium point is (1,1), and the yt vs. ytm1 curve has a local slope
of -1/4 = -0.25 at that point, whose magnitude is less than 1, indicating a true equilibrium
point. The negative sign of the local slope at point (1,1) means that the history of y[t] will
oscillate around y = 1 as convergence proceeds. 

We draw a phase diagram plot for this case, and add a line tangent to the curve at point
(1,1) whose slope is -1/4.

kill(yt, ytm1)$
yt : ytm1^(-1/4);
wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]", yrange = [0, 2], 
    proportional_axes = 'xy,  key = "phase curve y[t]", key_pos = top_center,
    explicit (yt, ytm1, 0, 2),color = red, key = "y[t] = y[t-1]", explicit (ytm1, ytm1, 0, 2), 
    color = dark_turquoise, key = "tangent",  explicit (-x/4 + 5/4, x, 0.5, 1.5),
    color = black, key = "", line_width = 1, explicit (1, ytm1, 0,1),
    parametric (1,yy,yy,0,1)), wxplot_size = [680, 680];

We can visualize the convergence in the case y[0] = y0 = 0.75.

kill(y0,y1,y2,y3,y4,f)$
f(x) := x^(-0.25);
[y0 : 0.75, y1 : f(y0), y2 : f(y1), y3 : f(y2), y4 : f(y3)];

We see that successive values of y (y1, y2, y3,...) oscillate between values larger than 1 
and values smaller than 1, but converging on y = 1.

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]",
    yrange = [0.5,1.5], key_pos = bottom_right, proportional_axes = 'xy,
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2), color = black, key = "", line_width = 1,
    parametric (1,yy,yy,0.5,1), head_length = 0.03, vector ([y0, 0], [0, y1]),
    vector ([y0,y1], [y1-y0, 0]), vector ([y1,y1], [0, y2-y1]), vector ([y1,y2],[y2-y1,0]),
    vector ([y2,y2], [0, y3-y2]), vector ([y2,y3], [y3-y2,0])), wxplot_size = [680,680];

If we start with a y0 larger than 1, say 1.25,

kill(y0,y1,y2,y3,y4,f)$
f(x) := x^(-0.25);
[y0 : 1.25, y1 : f(y0), y2 : f(y1), y3 : f(y2), y4 : f(y3)];

we again get oscillating convergence to y = 1.

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]",
    yrange = [0.5,1.5], key_pos = bottom_right, 
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2), color = black, key = "", line_width = 1,
    parametric (1,yy,yy,0.5,1), head_length = 0.03, vector ([y0, 0], [0, y1]),
    vector ([y0,y1], [y1-y0, 0]), vector ([y1,y1], [0, y2-y1]), vector ([y1,y2],[y2-y1,0]),
    vector ([y2,y2], [0, y3-y2]), vector ([y2,y3], [y3-y2,0])), wxplot_size = [680,680];

### Prob. 17.33

Construct a phase diagram for the nonlinear difference equation
      y[t] = y[t-1]^(-3/2)
and check for stability.

The only possible equilibrium point is (1,1) and the local slope at that point is -1.5,
with a magnitude greater than 1. Hence we have explosive (oscillating) instability.

kill(yt, ytm1)$
yt : ytm1^(-3/2);
wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]", key_pos = top_right,
    yrange = [0, 2.5], proportional_axes = 'xy,
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 2),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 2), color = black, key = "", line_width = 1,
    explicit (1, ytm1,0,2), parametric (1,yy,yy,0,1)), wxplot_size = [680, 680];

For a visual look at the explosive instability, we assume y[0] =  y0 = 0.75.

kill(y0,y1,y2,y3,y4,f)$
f(x) := x^(-1.5);
[y0 : 0.75, y1 : f(y0), y2 : f(y1), y3 : f(y2), y4 : f(y3)];

We see y1, y2, y3, ... values getting further and further away from y = 1, with oscillation 
between values larger than and smaller than y = 1.

wxdraw2d (xlabel = "ytm1 = y[t-1]", ylabel = "yt = y[t]",
    yrange = [0,4], key_pos = top_center, 
  key = "phase curve y[t]", explicit (yt, ytm1, 0, 4),color = red, key = "y[t] = y[t-1]",
    explicit (ytm1, ytm1, 0, 4), color = black, key = "", line_width = 1,
    explicit (1, ytm1,0,1), parametric (1,yy,yy,0,1), head_length = 0.03, 
    vector ([y0, 0], [0, y1]), vector ([y0,y1], [y1-y0, 0]), vector ([y1,y1], [0, y2-y1]), 
    vector ([y1,y2],[y2-y1,0]),vector ([y2,y2], [0, y3-y2]), 
    vector ([y2,y3], [y3-y2,0])), wxplot_size = [680,680];