# Dowling18C.wxmx: 
Inflation and Unemployment Models

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

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

In Dowling18C.wxmx, we use Maxima to discuss inflation and unemployment models.
Our treatment of continuous time models follows Chiang and Wainwright, Ch. 16, Sec. 5.
Our treament of discrete time models follows Chaing and Wainwright, Ch. 18, Sec. 3.
We have added sections dealing with differential equations with a variable term and
higher order ordinary differential equations (using desolve) following C&W 16.6,7.

We have changed some of the symbols used in particular problems. An approximate
pdf translation (using Microsoft print to pdf) is available as Dowling18Cfit.pdf. That pdf file
can be searched using Ctrl-F.

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

## References

1. Fundamental Methods of Mathematical Economics, Alpha C. Chiang and Kevin 
Wainwright, 4th ed., 2005, McGraw-Hill (Ch. 16 & 18)

2. Phillips Curve Discussion and plots:
https://courses.lumenlearning.com/boundless-economics/chapter/
                      the-relationship-between-inflation-and-unemployment/

3. Short Run Phillips Curve Video, 15 min
https://www.youtube.com/watch?v=vsbj2ppq7IQ

4. Long Run Phillips Curve Video, 10 min.
https://www.youtube.com/watch?v=GycD6uitKrA

5. Ch 30: Money Growth and Inflation
https://www.youtube.com/watch?v=1mXrsUkABWM

6. Ch. 35: The Short Run Tradeoff Between Inflation and Unemployment, 32 min video
by Jonathan Keisler
https://www.youtube.com/watch?v=VYr-GPdVYvg

7. Ch. 8, The Natural Rate of Unemployment and the Phillips Curve
Macroeconomics, Olivier Blanchard, Fifth Edition, 2011, Prentice Hall

8. Ch. 35, The Short-Run Tradeoff between Inflation and Unemployment, Principles of
Economics, Sixth Ed., N. Gregory Mankiw

## Inflation and Unemployment in Continuous Time

We review the continuous time modeling of the interaction between inflation and 
unemployment, as discussed in Chiang and Wainwright, Ch. 16, Sec. 5.

### The Phillips Relation

In the journal Economica (Nov. 1958), A.W. Phillips published a paper based on data relating
the rate of growth of money wages and unemployment in the U.K. over the period 1861-1957.

In Phillips' original formulation,  W(t) is the wage at time t, w = (1/W)*dW(t)/dt is defined as 
the "rate of growth of money wage" at time t,  U(t) is the "rate of unemployment", and
Phillips proposed the relation w = f(U) such that dw/dU < 0, which says that an increase
in unemployment normally means a decrease in the rate of growth of money wages for 
production workers, and conversely, a decrease in unemployment normally means an 
increase in the rate of growth of wages.

The unemployment rate U represents the number of unemployed as a percentage of the 
labor force. Labor force data are restricted to people 16 years of age and older, who 
currently reside in 1 of the 50 states or the District of Columbia, who do not reside in 
institutions (e.g., penal and mental facilities, homes for the aged), and who are not on active
duty in the Armed Forces. [units: percent; https://fred.stlouisfed.org/series/UNRATE]

The Phillips relation was then extended to finding a function that connects "the rate of
inflation" p = (1/P)*dP/dt (P is the price level), to the rate of unemployment U. This was 
based on the assumption that if w is positive, reflecting a growing money-wage cost, 
inflationary pressures on prices would appear, so the rate of inflation p should also be 
the same function of unemployment U: p = f(U) such that dp/dU < 0 (an increase in 
unemployment normally means a decrease in the rate of inflation).

As a math diversion, note that d(ln(x))/dt = (1/x)*dx/dt:

depends (x,t);

diff (log(x), t);

So we can write p = d(ln(P))/dt, and w = d(ln(W))/dt, as well.

The inflationary pressure of a positive value of w can be offset by an increase in labor
productivity T. Prices can remain stable in the face of increased wages paid by producers
provided each worker is able to produce goods at a greater rate. If money wage grows faster
than productivity, then inflationary pressure on prices should normally occur.

The simplest assumption is a linear relation between p, w, and T:
        p(t) = w(t)  - T,
and then the assumption that f(U) itself is some linear function of U:
        w(t) = α - β*U(t),
with α and β both positive constant parameters, so dw/dU = - β < 0,  we get a linear relation
between p and U:
       p(t) = α - T - β*U(t),       (α, β > 0).

### The Expectations-Augmented Phillips Relation

A more recent model of the relation between w and U is the "expectations-augmented 
version" of the Phillips relation:
          w(t) = f(U(t)) + g*π(t),  ( 0 < g <= 1),
where π(t) is the expected rate of inflation. As emphasized by Milton Friedman, if an 
inflationary trend has been present long enough, workers factor some future inflation 
expectations into wage demands, so w should be an increasing function of π. Our equation 
relating  rate of inflation p, unemployment U, and inflation expectation π is now
         p(t) = α - T - β*U(t) + g*π(t),      ( 0 < g <= 1),       (α, β > 0).            (1)

The "adaptive expectations hypothesis" is the equation
         dπ/dt = j*(p -  π),    (0 < j <= 1),                                                          (2)
which simply says that if the present rate of inflation is greater than the present expected
rate of inflation, the expected rate of inflation should increase with time. Likewise, if p falls
short of π, then π is revised in the downward direction. The driving force is the difference
between the actual and the expected rate of inflation.

### The Feedback from Inflation to Unemployment U(t)

The pair of equations (1) and (2) should be solved together. After taking the the labor
productivity T to be an exogenous constant, we have three time dependent variables:
(p, U, π) with only two equations. We must choose one of these three to also be exogenous
to make further progress. However if we make an assumption about the factors that directly 
influence unemployment and its dependence on time t, one approach is to take into account
only monetary policy and its influence on the rate of change of unemployment.

Let M(t) be the nominal money balance (the amount of money in the economy) and 
let m(t) = (1/M)*dM/dt be the "rate of growth of the money balance". A simple model which
relates U, p, and m is:
      dU/dt = -k*(m - p),   (k > 0).                                                                    (3)
This relation says that if the money growth rate is larger than the inflation rate, then
unemployment will decrease with time, and conversely, if the inflation rate is greater than
the money growth rate, unemployment will tend to increase.

The difference (m - p)  is the rate of growth of "real money". Eqn. (3) asserts an "interaction"
between the rate of inflation p and the "time path of unemployment U(t)", dU/dt.

### The Time Path of π, the Expected Rate of Inflation, croots(A,B)

Starting with the three equations (1), (2), (3), we can condense the model using substitution.
We first substitute (1) into (2) to get
    dπ/dt = j(α - T - β U) - j(1-g) π.                                                                     (4)
To be able to use (3) we differentiate (4) wrt time t (recall T = const.)
    d²π/dt² =  - j β dU/dt - j(1-g) dπ/dt.                                                                (5)
Using (3) for dU/dt in (5) results in
    d²π/dt² = j β k*(m - p) - j(1-g) dπ/dt.                                                               (6)
Solve eqn (2 ) for p in terms of π:
    p = (1/j) dπ/dt + π.                                                                                         (7)
Use (7) for p in (6);

    d²π/dt² + A*dπ/dt + B*π = C,                                                                          (8)

where
    A = β k + j(1-g),   B = j β k,   C = j β k m.                                                        (9)

We have written a linear second order ordinary differential equation for π (here with constant
positive coefficients, A > 0, B > 0, C > 0) in standard form, as we did in Dowling18A.wxmx.

The particular integral of (8) is π = constant, or B*π = C, or πe = C/B = m, so in this
model the intertemporal equilibrium value of the expected rate of inflation π equals the rate
of growth of nominal money.

We adopt a notation based on Dowling's treatment in his Ch. 18 and on our work in 
Dowling18A.wxmx, Dowling18Afit.pdf. The characteristic roots are the solutions of
the equation r^2 + A*r + B = 0,  and are called r1 and r2.

A Maxima function croots (A,B) is defined in Econ2.mac which returns a list [r1,r2].
We follow Chiang and Wainwright's conventions. Here we call croots with symbolic
symbols A and B (unbound to values):

fundef (croots);

croots (A,B);

Returning to our ode for π, Eq. (8) above, we have distinct real roots if A^2 > 4*B, 
double real roots if A^2 = 4*B, and complex roots if 4*B > A^2. There are no general reasons
for any of these relations to always be in effect, so in general we can have all three types 
of solutions.

### Case 1, Distinct Real Roots r1 and r2, A^2 > 4*B, y = π

Recall that both A and B are positive.

Case 1, distinct real roots, A^2 > 4*B. Since B > 0, (A^2 - 4*B) < A^2, and the square root
is less than A, hence r1 is negative, as is r2, which implies a dynamically stable equilibrium.
The time path of π(t) approaches  πe = m for "large" positive values of t.

Here is an example with A = 12, B = 1, C = 5, πe = C/B = 5. We let y stand for π.
We use the Maxima functions ode2 and ic2, as introduced in Dowling18A.wxmx.

We first use croots (A,B) to verify that both roots are negative, and hence a stable solution
will be found. [Recall that for this case, A^2 > 4*B, the solution of the complementary
equation is yc = %A1*exp(r1*t) + %A2*exp(r2*t).]

croots (12, 1), numer;

Using Lode2(y, t, type, A, B, C), defined in Econ2.mac, with type = real, double, or complex,
(here type = real) to get a general solution with two arbitrary constants %k1 and %k2:

gsoln : Lode2(y, t, real, 12,1,5);

For our first definite solution, we assume that at t = 0, y = 15 and y' = -2, and use ic2.

dsoln1 : ic2(gsoln,t = 0,y = 15, 'diff (y,t) = -2);

yex1 : rhs (dsoln1);

wxdraw2d (xlabel = "t", ylabel = "y", yrange = [0, 20],  explicit (yex1, t, 0, 50),
    color = red, line_width = 1, key = "ye = 5",  explicit (5,t, 0,50))$

For our second definite solution, we assume that at t = 0, y = 2 and y' = 0.2.

dsoln2 : ic2(gsoln,t = 0,y = 2, 'diff (y,t) = 0.2);

yex2 : rhs (dsoln2);

Show both solutions on the same plot.

wxdraw2d (xlabel = "t", ylabel = "y", yrange = [0, 20],explicit (yex1, t, 0, 50),
    color = brown, explicit (yex2, t, 0, 50),
    color = red, line_width = 1, key = "ye = 5",  explicit (5,t, 0,50))$

As Dowling discusses, the complementary solution to the linear second order ordinary
differential equation with real characteristic roots r1 and r2 has the general form
  k1*exp(r1*t)  + k2*exp(r2*t) 
which is why we need both r1 and r2 to be negative (in this real roots case) to have 
convergence to πe = 5. Let us construct such a solution using this form as a starting point,
and looking for the same solution as we found above with yex2 (at t = 0, y = 2, y' = 0.2).
This is an alternative to using ode2 or Lode2 for the two real distinct roots case.

[r1, r2] : croots (12, 1), numer;

We have to use %k1 and %k2 as the arbitrary constants in gsoln in order to get ic2 to not 
return an empty bracket [].

gsoln : y =  5 + %k1*exp (r1*t) + %k2*exp (r2*t);

dsoln : ic2(gsoln, t = 0, y = 2, 'diff(y, t) = 0.2), numer;

yex : rhs (dsoln);

wxdraw2d (xlabel = "t", ylabel = "yex", yrange = [0, 10],
    color = brown, explicit (yex, t, 0, 50),
    color = red, line_width = 1, key = "yex = 5", explicit (5, t, 0, 50))$

### Case 2: Double Real Root, A^2 = 4*B, y = π

If A^2 = 4*B, we have one real root r = -A/2.

croots (A,B);

ratsubst (4*B, A^2, %);

If A > 0, the single real root is negative, and we have convergence to a stable solution πe,
and conversely, if A < 0, the single real root is positive and we have an unstable evolution
in time of π(t).

### Case 3: Complex Roots, 4*B > A^2.

In this case, [r1,r2] = [g + %i*h, g - %i*h], with g = -A/2, h = (1/2)*sqrt(4*B - A^2) > 0.
The solution of the complementary equation is
  yc = exp(g*t)*[sum of cosine and sine terms], or
  yc = exp(-A/2)*(%k1*sin(h*t) + %k2*cos(h*t) ).

We can use our function Lode2(y,t,type,A,B,C) where type (here) is complex. We introduced
Lode2 in Dowling18A.wxmx.

Lode2 (y,t, complex, A, B, C);

In this case, the real part of the complex roots is g = -A/2, (the square root term is pure 
imaginary, the -A/2 term is the real part of each root) and, as in Case 2, we have a ,
stable converging solution if A > 0, otherwise an unstable solution.

Chiang and Wainwright (p. 536) consider a complex root case, starting with the basic
equations ( in the first eqn, we took α - T = 1/6, which says that T < α):
  p = 1/6 - 3 U + π,
  dπ/dt = (3/4)*(p - π),
  dU/dt = -(1/2)*(m - p).

Then the constants relevant for the values of A and B are β = 3, g = 1, j = 3/4, and k = 1/2.
Then A = 3/2 > 0 hence convergent, B = 9/8, C = (9/8)*m, A^2 = 9/4, 4*B = 9/2 > A^2, 
hence complex.

We can use our Lode2 function, defined in Econ2.mac, and used in Dowling18A.wxmx,
with the syntax Lode2 (y, t, type, A, B, C), with type = real, double, or complex (here
complex). We first use croots(A,B) with A = 3/2, B = 9/8:

croots (3/2, 9/8);

ypiSoln : Lode2(y, t, complex, 3/2, 9/8, 9*m/8);

ypi : rhs (ypiSoln);

Solving the second of the three starting equations for p, we get
p = (4/3)*dπ/dt + π.

p : (4/3)*diff (ypi,t) + ypi, ratsimp, expand;

Like the expected rate of inflation π, the actual rate of inflation p has a fluctuating time path
converging to the equilibrium value m.

Solving the first of our starting equations for U:
  U = (1/3)*(π - p) + 1/18.

U : (ypi - p)/3 + 1/18;

The time path of "rate of unemployment" U(t) is a damped fluctuating expression which 
eventually reaches Ue = 1/18.

### Long Run Unemployment Example

We see from this complex root example a prediction that the rate of unemployment reaches
an equilibrium value which does not depend on the monetary policy parameter m.

Quoting p. 537 of Chiang and Wainwright:

"Because the intertemporal equilibrium values of π and p are both equal to the monetary-
policy parameter m (the rate of growth of nominal money M: m = (1/M)*dM/dt), the value of
m provides the axis around which the time paths of π and p fluctuate. If a change occurs
in m, a new equilibrium value of π and p will immediately replace the old one, and whatever
values the π and p variables happen to take at the moment of the monetary-policy change
will become the initial values from which the new π and p paths emanate."

"In contrast, the intertemporal equilibrium Ue does not depend on m. In our model, U
converges to 1/18 regardless of the rate of growth of nominal money, and hence regardless
of the equilibrium rate of inflation p. This constant equilibrium value of U is referred to as
the "natural rate of unemployment". The fact that the natural rate of unemployment is
consistent with any equilibrium rate of inflation p can be represented in the (U, p) plane
by a vertical straight line (U = const) parallel to the vertical p axis. That vertical line is known
as the "long-run Phillips curve".

In the U.S. the natural rate of unemployment (NRU) is about 5%. A "natural rate of 
unemployment" greater than zero exists due to the natural churn in job expectations, short
term employment needs, and the search for more rewarding jobs.

Let's draw the long run Phillips Curve.

wxdraw2d (xlabel  = "Unemployment Rate (%)", xrange = [1, 9],yrange = [-3,9],
   ylabel = "Inflation Rate (%)", title = "Long Run Phillips Curve", 
    parametric (5, yy, yy, -3, 9), line_width = 1, color = red, explicit (2, x, 1, 5), 
    explicit (6, x, 1, 5), color = black, explicit(0,x,1,9),points ([[5,2],[5,6]]))$

### The Phillips Curve in Context

We quote a discussion from:
https://courses.lumenlearning.com/boundless-economics/chapter/
                                  the-relationship-between-inflation-and-unemployment/

"Key Points

"The relationship between inflation rates and unemployment rates is inverse. Graphically, 
this means the short-run Phillips curve is L-shaped."

"A.W. Phillips published his observations about the inverse correlation between wage 
changes and unemployment in Great Britain in 1958. This relationship was found to hold 
true for other industrial countries, as well."

"From 1861 until the late 1960’s, the Phillips curve predicted rates of inflation and rates of 
unemployment. However, from the 1970’s and 1980’s onward, rates of inflation and 
unemployment differed from the Phillips curve’s prediction. The relationship between 
the two variables became unstable."

"Key Terms
Phillips curve: A graph that shows the inverse relationship between the rate of 
unemployment and the rate of inflation in an economy.

stagflation: Inflation accompanied by stagnant growth, unemployment, or recession."

"The Phillips curve relates the rate of inflation with the rate of unemployment. The 
Phillips curve argues that unemployment and inflation are inversely related: as levels of 
unemployment decrease, inflation increases. The relationship, however, is not linear. 
Graphically, the short-run Phillips curve traces an L-shape when the unemployment 
rate is on the x-axis and the inflation rate is on the y-axis."

### Short Run Phillips Curve

We can construct our own plot of a "theoretical Phillips curve", which is designed to be
similar to the example shown in the online webpage.

kill(a,b)$
pp(x) := -a + b/cosh(x) - x;

psoln : solve([pp(3) = 6, pp(5) = 2], [a,b]), numer;

kill(U)$
IR : at(pp(U), psoln[1]);

[xmin, xmax,ymin,ymax] : [1,10,-5,10]$
wxdraw2d (xlabel = "Unemployment Rate (%)",
     ylabel = "Inflation Rate (%)", xrange = [xmin,xmax], yrange = [ymin,ymax],
    title = "Short Run Phillips Curve",
    explicit (IR,U,xmin,xmax), color = red, line_width = 1, 
    explicit (6, U, xmin, 3),explicit (2, U, xmin, 5),
    color = black, explicit(0,U, xmin, xmax), parametric(3,yy,yy, ymin, 6),
    parametric(5,yy,yy, ymin, 2), points ([ [3,6], [5,2] ]))$

"Short Run Phillips Curve: The Phillips curve shows the inverse trade-off between inflation 
and unemployment. As one increases, the other must decrease. In this example plot, an 
economy can either experience 3% unemployment at the cost of 6%  inflation, or 
increase unemployment to 5% to bring down the inflation levels to 2%."

"History of the Phillips Curve"

"The early idea for the Phillips curve was proposed in 1958 by economist A.W. Phillips. 
In his original paper, Phillips tracked wage changes and unemployment changes in 
Great Britain from 1861 to 1957, and found that there was a stable, inverse 
relationship between wages and unemployment. This correlation between wage changes 
and unemployment seemed to hold for Great Britain and for other industrial countries. In 
1960, economists Paul Samuelson and Robert Solow expanded this work to reflect the
relationship between inflation and unemployment. Because wages are the largest 
components of prices, inflation (rather than wage changes) could be inversely linked to 
unemployment."

"The theory of the Phillips curve seemed stable and predictable. Data from the 1960’s
modeled the trade-off between unemployment and inflation fairly well. The Phillips curve 
offered potential economic policy outcomes: fiscal and monetary policy could be used to 
achieve full employment at the cost of higher price levels, or to lower inflation at the cost 
of lowered employment. However, when governments attempted to use the Phillips curve 
to control unemployment and inflation, the relationship fell apart. Data from the 1970’s and 
onward did not follow the trend of the classic Phillips curve. For many years, both the 
rate of inflation and the rate of unemployment were higher than the Phillips curve would
have predicted, a phenomenon known as “stagflation. ” Ultimately, the Phillips curve was
proved to be unstable, and therefore, not usable for policy purposes."

The Relationship Between the Phillips Curve and AD-AS Curve

Changes in aggregate demand cause movements along the Phillips curve, all other 
variables held constant. Aggregate demand (AD) is the total demand for final goods and 
services in the economy at a given time and price level.

"The rate of unemployment and rate of inflation found in the Phillips curve correspond to 
the real GDP and price level of aggregate demand.

Changes in aggregate demand translate as movements along the Phillips curve.

If there is an increase in aggregate demand, such as what is experienced during 
demand-pull inflation, there will be an upward movement along the Phillips curve. 
As aggregate demand increases, real GDP and price level increase, which lowers the
unemployment rate and increases inflation."

"The Phillips curve shows the inverse trade-off between rates of inflation and rates of 
unemployment. If unemployment is high, inflation will be low; if unemployment is low, 
inflation will be high."

Points on the Phillips Curve must correspond to points on the AD-AS Curve. The Phillips 
curve is the relationship between inflation, which affects the price level aspect of aggregate 
demand, and unemployment, which is related to aggregate supply (AS). Consequently, 

"To see the connection more clearly, consider the example illustrated by" the two plots we
construct below. "Let’s assume that aggregate supply, AS, is stationary, and that aggregate
demand starts with the curve, AD1. There is an initial equilibrium price level and real GDP 
output at point A. Now, imagine there are increases in aggregate demand, causing the 
curve to shift right to curves AD2 through AD4. As aggregate demand increases, 
unemployment decreases as more workers are hired, real GDP output increases, and the 
price level increases; this situation describes a demand-pull inflation scenario."

"Phillips Curve and Aggregate Demand: As aggregate demand increases from AD1 to AD4, 
the price level and real GDP increases. This translates to corresponding movements along 
the Phillips curve as inflation increases and unemployment decreases."

Let's make an AD-AS plot with fixed supply and shifting demand.

### Plot of Shifting Demand with Fixed Supply

If we model the fixed supply curve as a straight line (for simplicity) P = a*Q + b, with P the
price level and Q the corresponding quantity supplied, we look for the parameter values a
and b such that when P = 100, Q = 5, and when P = 110, Q = 10.

solve ([100 = 5*a + b, 110 = 10*a + b],[a,b]);

Our supply curve is then described by the linear relation Ps = 2*Q + 90. which has a slope
dPs/dx = 2. For convenience let's use a downward sloping demand curve which is also a
straight line Pd = c*Q + d, and let the line Pd be orthogonal to Ps. A vector directed along
a line y = a*x + b is represented by the "x component" yx and the "y component" yy as 
a list [yx, yy] = [x2 - x1, y2 - y1] = [x2 - x1, a*(x2 - x1)] = (x2 - x1)*[1,a]. A vector directed along
a second line y = c*x + d is the list [x2 - x2, c*(x2 - x1)] = (x2 - x1)*[1, c]. For these two
direction vectors to be orthogonal, the scalar dot product must be zero.  

For example, represent the direction vector Avec by the list of x and y components Ax and
Ay, and represent the direction vector Bvec by the list of x and y components Bx and By.

kill(Ax,Ay,Bx,By)$
Avec : [Ax,Ay];
Bvec : [Bx,By];
AdotB : Avec . Bvec;

If the scalar dot product is equal to zero, these two direction vectors are orthogonal (directed
at right angles to each other). Note the use of a single dot (period) between the lists to find
the "dot product" of these two lists representing two direction vectors. In the Maxima manual,
under Dot, you find:
 "The operator . represents noncommutative multiplication and scalar product."
The use of lists of components to represent vectors in calculations is not documented in
the Manual.

solve( [1,a] . [1,c] = 0, c);

So we choose the slope of the demand curve line to be -1/2, and use Pd = -Q/2 + d, with
d to be chosen to have a particular (Q,P) point be on the line. If we want the point
(Q,P) = (5,100) to be on the line, then we need d to be

solve (100 = -5/2 + d, d), numer;

We do that calculation for all the points in the list SDpts.

SDpts : [[5,100],[6,102],[7,104],[9,108] ];

wxdraw2d(xlabel = "Real GDP", ylabel = "Price Level", yrange = [98,110],
    proportional_axes = 'xy,key_pos = top_left, line_width = 3,
    xrange = [4,10], color = red, key = "Supply", explicit (2*x+90,x,4.5, 9.5), 
    color = blue, key = "demand D1", explicit (-x/2 + 102.5,x,4.3,7),color = brown, 
    key = " demand D2", explicit(-x/2 + 105,x,4.5,8),color = magenta,  key = "demand D3",
     explicit (-x/2 + 107.5, x, 5,9),color = dark_turquoise,key = "demand D4",  
    explicit (-x/2 + 112.5,x, 7.6, 9.8),color = black,key = "", points (SDpts),
    label (["A", 4.9, 100.3], ["B",5.9,102.3], ["C",7.1,103.5], ["D",9.1,107.5],
    ["S",9.6,109.3], ["S",4.4,98.7])), wxplot_size=[1024,768];

"Phillips Curve and Aggregate Demand: As aggregate demand increases from AD1 to AD4, 
[curve A to curve D on this plot] the price level and real GDP increases. This translates 
to corresponding movements along the Phillips curve as inflation increases and 
unemployment decreases."

### Plot of Phillips Curve with Shifting Demand Points

We use Maxima's interpolation functions to draw an approximate Phillips Curve from sample
data. We define pts to be a list of elements [Unemployment-rate, Price Level], each of which
respectively correspond to the points (A,B,C,D) on the supply-demand curve drawn above.

pts : [[2,8], [3,4], [4,2], [6,0.6],[8,0]];

The contributed Maxima package interpol.mac has a function lagrange(points) defined which
creates a "Lagrange interpolation" on the points provided.

load ("interpol");

lagrange (pts);

We can then define a Maxima function f1(x) by placing two single quotes (') in front of %.

f1(x) := ''%;

We now show the points on this theoretical Phillips curve which correspond to the
points labeled (A,B,C,D) on the AD-AS curve drawn above.

wxdraw2d (xlabel  = "Unemployment Rate (%)", xrange = [1, 9],yrange = [-1,9],
   ylabel = "Inflation Rate (%)", title = "Theoretical Phillips Curve", 
    explicit (f1(x), x, 2,8), color = red, points (pts), 
    color = black, line_width = 1,explicit(0,x,1,9),
    label(["A", 8, 0.5], ["B", 4.2,2.5],["C", 3.2, 4.5], ["D", 2.3,8])), wxplot_size = [1024,768];

"Phillips Curve and Aggregate Demand: As aggregate demand increases from AD1(A) 
to AD4 (D), the price level increases and the real GDP increases. This translates to 
corresponding movements along the Phillips curve as inflation increases and unemployment
decreases."

### Unstable Inflation-Unemployment Example

From the webpage quoted above:

"The theory of the Phillips curve seemed stable and predictable. Data from the 1960’s
 modeled the trade-off between unemployment and inflation fairly well. The Phillips curve 
offered potential economic policy outcomes: fiscal and monetary policy could be used to 
achieve full employment at the cost of higher price levels, or to lower inflation at the cost 
of lowered employment. However, when governments attempted to use the Phillips curve 
to control unemployment and inflation, the relationship fell apart. Data from the 1970’s and 
onward did not follow the trend of the classic Phillips curve. For many years, both the 
rate of inflation and the rate of unemployment were higher than the Phillips curve would
have predicted, a phenomenon known as “stagflation. ” Ultimately, the Phillips curve was
proved to be unstable, and therefore, not usable for policy purposes."

The quoted webpage has a plot in which the relation between inflation and unemployment
looks random, and they remark:

"US Phillips Curve (2000 – 2013): The data points in this graph span every month from 
January 2000 until April 2013. They do not form the classic L-shape the short-run Phillips 
curve would predict. Although it was shown to be stable from the 1860’s until the 1960’s, 
the Phillips curve relationship became unstable – and unusable for policy-making – in the 
1970’s."

We will use Maxima tools to produce a similar plot, using random numbers to add white 
noise to the mean value.

We first make of plot of the mean signal only, using the function shown on the plot, which
we define as the Maxima function f(x) here. To match the webpage plot, we choose to draw
the line over the interval U = 3.8 to U = 10.

f(x) := -0.3556*x + 4.7393$
wxdraw2d(xlabel = "Unemployment (%)", ylabel = "Inflation Rate (%)",
    xrange = [3,11], yrange = [-3, 6],
    title = "Mean Signal: Unstable Inflation-Unemployment Period",
    explicit (f(x), x, 3.8,10), color = black,line_width = 1,
    explicit (0,x,3,11))$

The contributed Maxima package distrib contains a function we can use here.
random_normal(m, s) returns a random number selected from a normal distribution with
mean value m and standard deviation s. random_normal(m, s, n) returns a list of n such
random numbers. Use load("distrib") first to be able to use this function.

load("distrib");

For example, here is a list of ten random numbers taken from a normal distribution with
a mean value equal to zero (0) and with a standard deviation equal to unity (1).

makelist(random_normal(0,1), j, 1,10);

We can get ten numbers as a list by adding the option n = 10 to random_normal:

random_normal(0,1,10);

We now create a noisy signal by adding random numbers taken from a normal distribution
with a mean value of 0 and a standard deviation 1.7 to the mean signal f(x) we plotted above.
Let dx be the small change in U between white noise points; this will depend on the 
parameter num (the number of points).

num : 168$
dx : (10 - 3.8)/num$
ss : 1.7$
whiteNoise : makelist([3.8 + j*dx, f(3.8 + j*dx) + random_normal(0,ss)],j,0,num)$

Here we show this white noise as small points.

wxdraw2d(xlabel = "Unemployment (%)", ylabel = "Inflation Rate (%)",
    xrange = [3,11], yrange = [-3, 8],
    title = "Unstable Inflation-Unemployment Period, 2000-2013",
    explicit (f(x), x, 3.8,10), point_size = 0.5,
    points (whiteNoise),  color = black,line_width = 1,
    explicit (0,x,3,11))$

Repeating the conclusions:

"US Phillips Curve (2000 – 2013): The data points in this graph span every month from 
January 2000 until April 2013. They do not form the classic L-shape the short-run Phillips 
curve would predict. Although it was shown to be stable from the 1860’s until the 1960’s, 
the Phillips curve relationship became unstable – and unusable for policy-making – in the 
1970’s."

## Differential Equations with a Variable Term

We follow some of the discussion in Chiang and Wainwright, Sec. 16.6.

### Example 1: y'' + 5y' + 3y = 6t^2 - t - 1

Find a solution of the second order ODE
  y''(t) + 5 y'(t) + 3 y(t) = 6 t^2 - t - 1

This is not a linear ordinary differential equation (ode) with constant coefficients of the
form y'' + A y' + B y = C, because C is not a constant, so we can't use Lode2, but we 
can use ode2.

croots(5, 3), numer;

de : 'diff (y, t, 2) + 5*'diff (y,t) + 3*y = 6*t^2 - t -1;
gsoln : ode2 (de, y, t);

The "particular solution" yp = 2 t^2 - 7 t + 10.

We now assume y(0) = 3, y'(0) = 0.5, to find a definite solution.

dsoln : ic2 (gsoln, t = 0,y = 3, 'diff (y, t) = 0.5), numer;

yex1 : rhs (dsoln);

Ex1(t) := ''yex1;

wxdraw2d (xlabel = "t", ylabel = "y(t)", yrange = [0,14],
    explicit (Ex1(t), t, 0, 4))$

We can look for the value of t for which y'(t) = 0 and obtain the minimum value of y.

Ex1(0);

map ('Ex1, [1.1, 1.2, 1.3,1.4]);

So roughly, tmin = 1.3. We can get a more accurate value of tmin by looking for the value
of t for which dEx1(t)/dt = 0.

dy : diff(Ex1(t), t);

solve can't handle this type of expression, so use purely numerical method using
find_root(expr, var, v1, v2).

tmin : find_root (dy, t, 1.2, 1.4);

ymin : Ex1(tmin);

wxdraw2d (xlabel = "t", ylabel = "y(t)", yrange = [0,14], key_pos = top_left,
    title = "y'' + 5*y + 3*y = 6*t^2 - t - 1", key = "y(t)",
 explicit (Ex1(t), t, 0, 4), color = red, key = sconcat("ymin = ", ymin),explicit (ymin, t, 0, 4))$

### Example 3: y'' + 3y' - 4y = 2 exp(-4t)

Find the solution to the second order ODE
 y''(t) + 3 y'(t) - 4 y = 2 exp(-4 t).

croots (3,-4);

One of the real distinct roots of the characteristic equation is positive, and will dominate the
solution for long enough times, leading to an unstable solution.

de : 'diff (y,t,2) + 3*'diff (y,t) - 4*y = 2*exp (-4*t);
gsoln : ode2 (de, y, t);

Look for a definite solution for which y(0) = 5, y'(0) = -5.

dsoln : ic2 (gsoln, t = 0, y = 5, 'diff (y,t) = -5);

yex3 : rhs (dsoln);

Ex3(t) := ''yex3;

wxdraw2d (xlabel = "t", ylabel = "y(t)", yrange = [4,9],
    explicit (Ex3(t), t, 0, 1))$

dy : diff (Ex3(t),t);

tmin : find_root (dy, t, 0.1, 0.3);

ymin : Ex3(tmin);

wxdraw2d (xlabel = "t", ylabel = "y(t)", yrange = [4,9], key_pos = top_left,
    key = "y(t)", explicit (Ex3(t), t, 0, 1),color = red, key = sconcat (" ymin = ", ymin),
    explicit (ymin, t, 0, 1))$

## Higher Order ODE's With desolve (eqn, [y(t)])

We follow Sec. 16.7 in Chiang and Wainwright.

ode2 cannot handle third order or higher ODE's, so we must use desolve, which
has a slightly different syntax. We presented an example of using desolve at the end of
Dowling18A.wxmx.

### Example 4: y''''(t) + 6 y'''(t) + 14 y''(t) + 16 y'(t) + 8 y = 24

Notice that we must write 'diff(y(t),t,4), and not 'diff (y,t,4) when defining the differential
equation de. If we don't use atvalue to define initial values, we get a formal solution.

de :  'diff(y(t),t,4) + 6*'diff (y(t),t,3) + 14*'diff (y(t),t,2) + 16*'diff (y(t),t) + 8*y(t) = 24;
gsoln : desolve (de, [y(t)]);

To get a solution corresponding to initial values, use atvalue.

atvalue ('diff (y(t),t,3), t = 0, 0.5);
atvalue ('diff (y(t),t,2), t = 0, 1);
atvalue ('diff (y(t), t), t = 0, 2);
atvalue (y(t), t = 0, 4);
dsoln : desolve (de, [y(t)]);

yex4 : rhs (dsoln);

Ex4(t) := ''yex4;

wxdraw2d (xlabel = "t", ylabel = "y(t)", yrange = [2,6], explicit (Ex4(t), t, 0, 12))$

dy : diff(yex4, t);

tmax : find_root (dy, t, 0.5, 1.5);

ymax : Ex4(tmax);

tmin : find_root (dy, t, 4.5, 5.5);

ymin : Ex4(tmin);

wxdraw2d (xlabel = "t", ylabel = "y(t)", yrange = [0,6], key_pos = bottom_right,
    key = "y(t)", explicit (Ex4(t), t, 0, 12),color = red, line_width = 1, key = "ye = 3",
    explicit (3, t, 0, 12), color = dark_turquoise, key = sconcat("ymin = ",ymin),
    explicit (ymin, t, 0, 12), color = brown, key = sconcat ("ymax = ", ymax),
    explicit (ymax, t, 0, 12))$

## 2nd Order Difference Eqn., Char. Eqn. Roots r1, r2

We follow Chiang and Wainwright, Sec. 16.6. except for notation!
Given the 2nd order difference equation with constant coefficients
   y[t]  + b1 y[t-1] + b2 y[t-2] = a,  equivalent to y[t+2]  + b1 y[t+1] + b2 y[t] = a,
the characteristic roots are the solutions of the characteristic equation r^2 + b1 r + b2 = 0.
We call these characteristic roots r1 and r2.

croots (b1, b2) returns the roots of the characteristic equation and is defined in Econ2.mac.

[r1,r2] : croots (b1,b2);

grind(%)$

r1 and r2 then are related to b1 and b2 by the following:

1.   r1+r2 = -b1
2.   r1*r2 = b2
3.   (1 - r1)*(1 - r2) = 1 + b1 + b2.

The third relation is easy to derive by hand, using the first two relations. Getting Maxima to
do the simplifications requires the use of ratsimp, or first expand, followed by ratsimp.

r1 + r2, ratsimp;

r1*r2, ratsimp;

(1-r1)*(1-r2), expand, ratsimp;

## Inflation and Unemployment in Discrete Time

We follow Chiang and Wainwright, Sec. 18.3.
With p the actual inflation rate (%), π the expected inflation rate (%), and U the 
unemployment rate (%), m the fixed rate of increase of money [m = (1/M)*dM/dt]
set as a monetary policy decision by the FED, we translate the continuous time
equations used above into a discrete time version. Note that with Δt = 1 unit, 
Δx/Δt = (x[t+Δt] - x[t])/Δt =  x[t+1] - x[t], and we are replacing dx/dt ~ Δx/Δt = x[t+1] - x[t].


The continuous time equation
     p(t) = α - T - β*U(t) + g*π(t),      ( 0 < g <= 1),       (α, β > 0),                                (1)
is written as (t = 0, 1, 2, ...):
     p[t] = α - T - β*U[t] + g*π[t],      ( 0 < g <= 1),       (α, β > 0).                               (10) 

The continuous time equation
      dπ/dt = j*(p -  π),    (0 < j <= 1),                                                                          (2)
is written as
      π[t+1] - π[t] = j*(p[t] -  π[t]),    (0 < j <= 1).                                                         (11)

The continuous time equation
      dU/dt = -k*(m - p(t)),   (k > 0),                                                                             (3)
is written as
      U[t+1] - U[t] = -k*(m - p[t+1]),   (k > 0).                                                               (12)

Eqn (12) assumes the change in U[t] over the time interval [t, t + Δt] depends on 
(m - p[t+1]), the rate of growth of "real money" in period (t+1). Chiang and Wainwright 
leave as a problem the consequences of assuming the change in U[t] depends on 
(m - p[t]) instead.

### Time Path Difference Equation for Inflation p[t], IUparam

Rewrite (10) as:
   p[t+1] = α - T - β*U[t+1] + g*π[t+1],                                                                       (10')
then subtract (10) from (10') to get the first difference of p[t]:
   p[t+1] - p[t] = - β*(U[t+1] - U[t]) + g*(π[t+1] - π[t]).
Use (12) for the first difference of U on the rhs and (11) for the first difference of π on the
rhs to get
   p[t+1] - p[t] = β k (m - p[t+1])  + g j (p[t] -  π[t]).
Solve (10) for g*π[t]:
  g*π[t] = p[t] - α + T + β*U[t], 
and use this to get, after collecting terms,
     A*p[t+1] - B*p[t] + j*β*U[t] = C, 
     with A = (1 + β*k), B = [1 - j*(1 - g)], C = j*(α - T) + β*k*m.
which implies
     A*p[t+2] - B*p[t+1] + j*β*U[t+1] = C, and after subtraction, we get
         
A*(p[t+2] - p[t+1]) - B*(p[t+1] - p[t]) + j*β*(U[t+1] - U[t]) = 0

Now replace the difference in U[t] using (3) to get, after collecting terms,

(1 + β*k)*p[t+2] - [1 + g*j + (1 - j)*(1+β*k)]*p[t+1] + [1 - j*(1 - g)]*p[t] = j*β*k*m.

We write this in the standard form we used in Dowling18B.wxmx:
      
         p[t] + b1*p[t-1] + b2*p[t-2] = a,                                                                           (13)

in which  b1 = - [1 + g*j + (1 - j)*(1+β*k)]/(1 + β*k),
               b2 = [1 - j*(1 - g)]/(1 + β*k),                                                                          (14)
               a = j*β*k*m/(1 + β*k).
  
Note that 0 < b2 < 1/(1+β*k) < 1, given, 0 < g <= 1, 0 < j <= 1, k > 0, and β > 0. 
This implies that  sqrt(b2) < 1.  We discuss 0 < b2 < 1 in more detail below.

Define the list IUparam, taking care to replace in Eq. (14) brackets [, ], by parentheses (,  ).

IUparam : [b1 = - (1 + g*j + (1 - j)*(1+β*k))/(1 + β*k), b2 = (1 - j*(1 - g))/(1 + β*k),
                        a = j*β*k*m/(1 + β*k)];

grind(%)$

The intertemporal equilibrium value of p is the particular integral of (13), which can be
found using our Maxima function ypart (b1, b2, a, t) provided a is independent of t.

yp : ypart(b1,b2,a,t);

and we can use the list IUparam to express this, after simplification, in terms of our
original parameters.

at(yp, IUparam),ratsimp;

As in the continuous time analysis, the equilibrium rate of inflation p = (1/P)*dP/dt is
exactly equal to the rate of monetary expansion m = (1/M)*dM/dt. If m is 5%, then,
after waiting for the smoke to clear, p = 5%.

The complete solution, as we saw in Dowling18B.wxmx, has a form which depends on the
relative values of b1^2 and 4*b2. If b1^2 > 4*b2 then we have distinct real roots of the
characteristic equation r^2 + b1*r + b2 = 0. For this real root case, if the magnitude of the
dominant root is greater than 1, the time path is divergent, since yc = %A1*r1^t + %A2*r2^t.
 
If 4*b2 > b1^2 we have complex characteristic roots, and k = sqrt(b2) must be less than 1 
for convergence, since yc = k^t*(oscillating fac). We have already seen above that
0 < b2 < 1, so k = sqrt(b2) < 1, and complex solutions are convergent.

### Review of croots, ypart, ycompl, ysoln, ytlist(b1,b2,a,y0,y1,tmax)

As a mini-review of second order difference equation solutions, consider the case
of Dowling Prob. 18.18, which has b1 = 7, b2 = 6, a = 42, b1^2 = 49 > 4*b2 = 24, so there
are two real (distinct) roots. Below we use Maxima functions defined in Econ2.mac.

croots (b1, b2) returns the roots of the characteristic equation r^2 + b1*r + b2 = 0.

croots(7,6);

The particular solution ypart(b1,b2, a, t), even if a is constant, sometimes depends on t.
If a does depend on t, ypart will not find the correct solution. (See our later section on
2nd order linear difference equations with a variable term.)

ypart (7, 6, 42, t);

The complementary equation solution ycompl (b1,b2,t) contains two arbitrary constants 
%A1 and %A2, whose value depends on initial conditions (for example).

ycompl(7,6, t);

The indefinite solution ysoln (b1,b2, a, t) is the sum of the particular solution and the 
solution of the complementary equation. Our function ysoln calls ypart and ycompl and
will not return the correct solution if the right hand side 'a' depends on time t.

ysoln(7, 6, 42, t);

The dominant root is -6, whose magnitude is greater than 1, hence any definite solution
will diverge from ye = 3.

In Dowling18B.wxmx we show you how to evaluate the constants for various initial
conditions.

A shorter route to viewing the solution behavior is to use ytlist(b1,b2,a,y0,y1,tmax), which
will produce the evolution of y[t] whether or not 'a' depends on time.

Suppose we assume y[0] = 5, y[1] = 4. ytlist(b1,b2,a,y0,y1,tmax) produces a list of values
of y[t] starting with y[0].

ytL : ytlist(7,6,42, 5, 4, 4);

Obviously divergent.

tL : makelist (t, t, 0, 4);

wxdraw2d(xlabel = "t", ylabel = "y[t]", yrange = [-1e3, 200],
    points (tL, ytL), color = black, line_width = 1, explicit(0, t, 0, 4))$

### b1^2 - 4*b2  <---    R(g, j, kbeta) 

We define R (g, j, kbeta) to be the difference b1^2 - 4*b2, expressed in terms of g, j, and
kbeta = k*β = β*k.

Rexpr : at (b1^2 - 4*b2, IUparam);

Rexpr : ratsubst (kbeta, β*k, Rexpr);

R(g, j, kbeta) := ''Rexpr;

### Real Roots Example for p[t]

We consider the case g = 1/2, j = 1/3, β*k = 5, m = 2. We can find the b1, b2, and 'a' values
associated with the 2nd order difference equation for p[t]: p[t] + b1*p[t-1] + b2*p[t-2] = a.

R (1/2, 1/3, 5),numer;

Since R > 0, b1^2 > 4*b2, and we have real roots.

IUparam;

We use ratsubst to replace β*k, and use a list Case_p1 to replace the other three 
parameters.

Case_p1 : [g = 1/2, j = 1/3, m = 2];

at (IUparam, Case_p1);

ratsubst (5, β*k, %);

Check that the long term value of p[t] is m = 2.

at (a/(1 + b1 + b2), %);

Check that the roots are real.

croots (-31/36, 5/36), numer;

The roots of the characteristic equation are real, positive, and have magnitudes < 1, 
implying a convergent solution case.

tL : makelist (t, t, 0, 5);

We can now use ytlist (b1, b2, a, y0, y1, tmax), letting p[0] = 4, p[1] = 3.5, tmax = 5.

ptL1 : ytlist (-31/36, 5/36, 5/9, 4, 3.5, 5);

A second list ptL2 with p[0] = 0.5, p[1] = 1.

ptL2 : ytlist (-31/36, 5/36, 5/9, 0.5, 1, 5);

wxdraw2d  (xlabel = "t", ylabel = "p[t]", xrange = [0,6], yrange = [0, 5], 
    title = " Inflation Rate p[t] ", points_joined = true,
    points (tL, ptL1), points (tL, ptL2), color = red, line_width = 1, 
    key = "m = 2",explicit (2, t, 0, 6))$

which shows convergence to m = 2.

### Complex Roots Example for p[t]

Now consider an example with g = 0.9, j = 0.9, β*k = 5, and m = 2, using two different 
starting conditions (y0, y1), and going out to t  = 5. This case produces complex 
characteristic roots, so the convergence includes some oscillation.
Check with R (g, j, kbeta).

R (0.9, 0.9, 5),numer;

This shows that b1^2 < 4*b2, so we have complex roots.

IUparam;

Case_p2 : [g = 0.9, j = 0.9, m = 2];

at (IUparam, Case_p2);

ratsubst (5, β*k, %);

Check that the long term value of p[t] is m = 2.

at (a/(1 + b1 + b2), %);

croots (-241/600, 91/600);

The magnitude of K = sqrt(b2) must be less than 1 for convergence in a complex root case.

sqrt(91/600), numer;

A complex root solution has the form K^t * (oscillating factor), where K = sqrt(b2).
If the magnitude of K is less than 1 the solution converges in the intertemporal value
pe = m = 2. The above cell shows we have convergence. Since the roots are complex,
the convergence will show some oscillation around pe = 2.

Using ytlist (b1, b2, a, p0, p1, tmax) with p[0] = 4, p[1] = 3.5, tmax = 5, p decreasing initially.

ptL1_complex : ytlist (-241/600, 91/600, 3/2, 4, 3.5, 5);

A second solution has p[0] = 0.5, p[1] = 1, tmax = 5, p increasing initially.

ptL2_complex : ytlist (-241/600, 91/600, 3/2, 0.5, 1, 5);

Here we show a plot of both initially decreasing and initially increasing.

wxdraw2d  (xlabel = "t", ylabel = "p[t]", xrange = [0,6], yrange = [0, 5], 
    title = " Inflation Rate p[t] ", key = "p[0] = 4", points (tL, ptL1_complex), 
    color = magenta, key = "p[0] = 0.5", points (tL, ptL2_complex), 
    color = black, line_width = 1,  key = "m = 2",explicit (2, t, 0, 6))$

With complex characteristic roots, the convergence includes some oscillation.

### General Characteristic Root Properties for Time Path of p[t]

We start with the general case.

[r1,r2] : croots (b1,b2);

First look at the sum of the roots.

r1 + r2, ratsimp;

Now we look at our case.

r1pr2 : at(ratsimp(r1+r2), IUparam);

We can absorb the overall negative sign out front into a multiplication of the numerator by
that minus sign using num(expr)/denom(expr) (this doesn't always work!).

num(r1pr2)/denom(r1pr2);

With β > 0, k > 0, 0 < g <= 1, 0 < j <= 1, both the denominator and numerator are positive
and (r1 + r2) > 0.

Now look at the product of the roots.

r1*r2, ratsimp;

at(b2, IUparam);

Again, both the numerator and denominator are positive for all allowed values of β, g, k, j.
Hence r1*r2 > 0. We can also conclude that r1*r2 < 1, as follows:
First we have a positive denominator which reduces the magnitude of the positive numerator,
and hence b2 > 0 and r1*r2 > 0.

a) If g = 1, the numerator = 1, r1*r2 = 1/(+) < 1.
b) If g < 1, consider first j = 1 also, then the numerator = 1-(1-g)  which is positive but less
than one,  so r1*r2 = (<1)/(+) < 1.
Consider next both g < 1 and j < 1, j(1-g) < j < 1, the numerator is a positive number less 
 one and the fraction is even smaller. So we have 0 < b2 < 1 and 0 < r1*r2 < 1.

What generally can we say about (1 - r1)*(1 - r2)?

(1-r1)*(1-r2), expand, ratsimp;

at(1+b1+b2,IUparam),expand,ratsimp;

Both the numerator and denominator are positive, so we have the third general property
        (1-r1)*(1-r2) > 0

To summarize, for our problem, the roots r1 and r2 of the characteristic equation have
the three properties:
     P1:  r1 + r2 > 0,
     P2:   0 < r1*r2 < 1,
     P3:  (1 - r1)*(1 - r2) > 0.

### Convergence Proved for Distinct Real Roots Case

For case 1, b1^2 > 4*b2, we have two distinct real roots.
For this case if the magnitude of the dominant root is > 1, the time path is divergent.

It turns out that for this inflation-unemployment model, we can prove that both roots
are positive and less than 1, so that the time path is convergent for any choice of 
parameters (g, j, k, β). (We have already concluded the complex root case is always
convergent.)

Since property P2 includes the requirement that r1*r2 > 0, both roots must have the 
same sign (since if the signs of the roots were opposite, the product would be negative).

If both roots were negative, their sum would be negative, violating P1. Hence both roots must
be positive. This means that no oscillation of the solution around ye occurs.

If either root were equal to one, a violation of P3 would occur, so neither r1 or r2 can equal
one.

It is also impossible that one root is greater than one and the other root is less than one, for
that would result in (1-r1)*(1-r2) < 0 in violation of property P3.

It follows that both roots (which must be positive) are either less than one, or both roots
are greater than one.

However we can rule out the latter case, since if r1 > 1 and r2 > 1 then r1*r2 > 1 which 
violates property P2.

Hence each root must be positive and less than one. This means than we have smooth 
convergence to ye.

We conclude that 0 < r1 < 1 and 0 < r2 < 1 and hence the time path of p[t] is convergent
for case 1, b1^2 > 4*b2.

### Time Path Difference Equation for Unemployment U[t]

This will take some work.

We repeat our starting set of discrete time relations
     p[t] = α - T - β*U[t] + g*π[t],      ( 0 < g <= 1),       (α, β > 0).                               (10) 

      π[t+1] - π[t] = j*(p[t] -  π[t]),    (0 < j <= 1).                                                         (11)

      U[t+1] - U[t] = -k*(m - p[t+1]),   (k > 0).                                                               (12)

We replace p[t+1] in (12) by solving for p[t] in (11) and translating in time to get
 (1 + β*k) U[t+1] - U[t] = k*(α - T - m) + k*g*π[t+1].
We difference this equation by replacing t by t+1 everywhere,
 (1 + β*k) U[t+2] - U[t+1] = k*(α - T - m) + k*g*π[t+2], and subtracting to get

   (1+β*k) (U[t+2] - U[t+1]) - U[t+1] + U[t] = k*g*(π[t+2] - π[t+1]).                              (12a)

On the left hand side the terms proportional to U[t+1] can be combined as
   -(2 + β*k) U[t+1], leaving
(1+β*k) U[t+2]  -(2 + β*k) U[t+1] + U[t] = k*g*(π[t+2] - π[t+1]).                                  (12b)

Replace t by t+1 throughout in (11) to get
 π[t+2] - π[t+1] = j*(p[t+1] -  π[t+1]).                                                                         (11')
Use this in the rhs of (12b) to get
(1+β*k) U[t+2]  -(2 + β*k) U[t+1] + U[t] = k*g*j*(p[t+1] -  π[t+1])                            
                                                          =  k*g*j*p[t+1]  -  k*g*j*π[t+1]).                      (12c) 
Solve (12) for k*p[t+1]:
k*p[t+1] = U[t+1] - U[t] + k*m.                                                                                  (12p)

Solve (10) for g*π[t+1],
g*π[t+1] = p[t+1] - α + T + β*U[t+1],
so
-k*g*j*π[t+1] = -k*j*(p[t+1] - α + T + β*U[t+1]) 
                    = -k*j*p[t+1]  -k*j*(- α + T + β*U[t+1])
                    = -k*j*p[t+1]  + k*j*(α - T) - k*j*β*U[t+1]
and the rhs of (12c) becomes
rhs(12c) = k*g*j*p[t+1]  -k*j*p[t+1]  + k*j*(α - T) - k*j*β*U[t+1]
             = (g - 1)*j*k*p[t+1] + k*j*(α - T) - k*j*β*U[t+1], now use (12p):
             = (g - 1)*j*(U[t+1] - U[t] + k*m) + k*j*(α - T) - k*j*β*U[t+1]
             =  j*(g -1- k*β)*U[t+1] - j*(g - 1)*U[t]  + j*k*m*(g - 1) + k*j*(α - T).

Using this rhs(12c) in (12c) we get
(1+β*k) U[t+2]  -(2 + β*k) U[t+1] + U[t] = j*(g -1- k*β)*U[t+1] - j*(g - 1)*U[t] 
                                                                             + j*k*m*(g - 1) + k*j*(α - T)        (12d)

Setting all U[t+1] terms on the lhs, those terms are
( - 2 - β*k - j*g + j + j*k*β)*U[t+1].
Bringing U[t] also to the left side we get for (12d)

 (1+β*k) U[t+2]  +  ( - 2 - β*k - j*g + j + j*k*β)*U[t+1]  + 
                                                                (1 + j*(g - 1))*U[t]  = j*k*m*(g - 1) + k*j*(α - T)
which can be rewritten as
(1+β*k) U[t+2]  -  [1 + g*j + (1 - j)*(1 + β*k)]*U[t+1] + [1 - j*(1 - g)]*U[t]  = 
                                                                                  k*j*[α - T  - (1 - g)*m].          (12e)

Dividing by (1+β*k) and replacing t by (t-2) everywhere to get our standard form,
 
     U[t] + b1*U[t-1] + b2*U[t-2] = a,                                                                         (15)
   
in which the constants (b1, b2, a)  are:

     b1 = -  [1 + g*j + (1 - j)*(1 + β*k)]/(1+β*k),

     b2 =  [1 - j*(1 - g)]/(1+β*k),                                                                                (16)      

     a =   k*j*[α - T  - (1 - g)*m]/(1+β*k).

The expressions for b1 and b2 found for the difference equation governing U[t] are exactly
the same as those found for the difference equation governing the difference equation
for p[t].  Our convergence arguments then apply to U[t] solutions also. 

However, the value of a on the right hand side of the U[t] difference equation is not the
same as that on the right hand side of the p[t] difference equation. There is no particular
reason why the intertemporal value of inflation rate should be the same as the intertemporal
value of the unemployment rate.

We define uIUparam as a list containing the same b1 and b2 used with IUparam above,
but with a different definition of a. To distinguish the value of the a parameter for U[t] from
the a parameter for p[t], we us the symbol aU.

uIUparam : [b1 = - (1 + g*j + (1 - j)*(1+β*k))/(1 + β*k), b2 = (1 - j*(1 - g))/(1 + β*k),
                        aU = k*j*(α - T  - (1 - g)*m)/(1+β*k)];

To identify the long term intertemporal value of U we look for the particular solution of the
difference equation for U[t], (15), by assuming that U[t] is a constant C.
Then C + b1*C + b2*C = aU, and C = aU/(1 + b1 + b2).

We can also use our function ypart (b1,b2,aU,t), since aU is a constant, independent of t.

yp : ypart (b1,b2, aU, t);

Ue : at (yp, uIUparam), ratsimp;

The long run unemployment Ue thus can be written as
   Ue = [α - T - (1 - g)*m]/β  = [α - T - (1 - g)*pe]/β,                                                    (17)
since the long run rate of inflation is pe = m.

This equation relates the long run unemployment with the long run inflation rate, and is
called the long run Phillips relation.

In the special case g = 1, the dependence of Ue on pe disappears, and Ue will be the same
irregardless of the long-run rate of inflation. The resulting expression then predicts the 
"natural rate of inflation":
     Unat = (α - T)/β.                   

If g = 1, the long run Phillips Curve is a vertical line parallel to the long-run inflation rate
 (pe) axis, showing the same long-run unemployment no matter the long-run inflation rate pe.

If g < 1, the coefficient of pe in (17) is negative and d(Ue)/d(pe) = -(1-g)/β < 0, and an
increase in the long run inflation rate will lead to a decrease in the long term unemployment.

Quoting Chiang and Wainwright p. 585:
"Then the long run Phillips Curve will turn out to be downward sloping, thereby still providing
a trade-off relation between inflation and unemployment. Whether the long-run Phillips
curve is vertical or negatively sloped is, therefore, critically dependent on the value of the
g parameter, which, according to the expectations-augmented Phillips relation":
    p(t) = α - T - β*U(t) + g*π(t),      ( 0 < g <= 1),       (α, β > 0)
first used in the continuous time treatment of our inflation-unemployment model.
The value of the g parameter "measures the extent to which the expected rate of inflation π
can work its way into the wage structure and the actual rate of inflation p."

### Case 1: Real Roots Example for U[t]

Consider the case with parameter choices: β*k = 5, β = 2, k = 5/2, g = 1/2, j = 1/3, 
α - T = 11, m = 2. We use ratsubst in a separate step to replace α - T = 11.

The choices β*k = 5, g = 1/2, and j = 1/3 agree with our real roots example for p[t] above, 
and we can use the  lists ptL1 (p[0] = 4, p[1] = 3.5, with p[t] decreasing to m = 2) and 
ptL2 (p[0]=0.5,p[1]=1, with p[t] increasing) here to plot both p[t] converging to m = 2,
and U[t] converging to Ue = 5.

Note that in the short run, unemployment U[t] can be increasing (if U[0] < 5) or decreasing
(if U[0] > 5), whether or not p[t] is decreasing (p[0] > 2) or increasing (p[0] < 2).

Ucase1 : [β = 2, k = 5/2, g = 1/2, j = 1/3,  m = 2];

at (uIUparam, Ucase1);

ratsubst(11, α - T, %);

croots(-31/36, 5/36);

float(%);

at (Ue, Ucase1);

ratsubst (11, α - T, %);

tL;

Let UtL1 correspond to U[t] increasing, using U[0] = 1, U[1] = 1.5, tmax = 5,
calling ytlist (b1, b2, aU, y0, y1, tmax).

UtL1 : ytlist (-31/36, 5/36, 25/18, 1,1.5, 5);

Likewise let UtL2 correspond to U[t] decreasing to Ue = 5, using U[0] = 7, U[1] = 6.5, tmax=5.

UtL2 : ytlist (-31/36, 5/36, 25/18, 7, 6.5, 5);

Two plots using solution ptL1 (p decreasing). First using solution UtL1 (U increasing)

ptL1;

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
    title = "short-run inflation decreasing, unemployment increasing",
    key = "U[t]", points (tL, UtL1), color = magenta, key = "p[t]",
    points (tL, ptL1), color = black, line_width = 1, key = "Ue = 5",
    explicit (5, t, 0, 6), color = brown, key = "m = 2", explicit (2, t, 0, 6))$

Next ptL1 (p decreasing) with UtL2 (U decreasing).

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
    title = "short-run inflation decreasing, unemployment decreasing",
    key = "U[t]", points (tL, UtL2), color = magenta, key = "p[t]",
    points (tL, ptL1), color = black, line_width = 1, key = "Ue = 5",
    explicit (5, t, 0, 6), color = brown, key = "m = 2", explicit (2, t, 0, 6))$

Next two plots using list ptL2 (p[t] increasing).

ptL2;

First ptL2 (p increasing) with UtL2 (U decreasing).

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
    title = "short-run inflation increasing, unemployment decreasing",
    key = "U[t]", points (tL, UtL2), color = magenta, key = "p[t]",
    points (tL, ptL2), color = black, line_width = 1, key = "Ue = 5",
    explicit (5, t, 0, 6), color = brown, key = "m = 2", explicit (2, t, 0, 6))$

Next ptL2 (p increasing) with UtL1 (U increasing).

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
    title = "short-run inflation increasing, unemployment increasing",
    key = "U[t]", points (tL, UtL1), color = magenta, key = "p[t]",
    points (tL, ptL2), color = black, line_width = 1, key = "Ue = 5",
    explicit (5, t, 0, 6), color = brown, key = "m = 2", explicit (2, t, 0, 6))$

### Case2: Complex Roots Example for U[t]

We consider the case g = 0.9, j = 0.9, β*k = 5, 

Consider the case with parameter choices: β*k = 5, β = 2, k = 5/2, g = 0.9, j = 0.9, 
α - T = 11, m = 2. We use ratsubst in a separate step to replace α - T = 11.

The choices β*k = 5, g = 0.9, and j = 0.9 agree with our complex roots example for p[t] above, 
and we can use the  lists ptL1_complex (p[0] = 4, p[1] = 3.5, with p[t] decreasing to m = 2)
and ptL2_complex (p[0]=0.5,p[1]=1, with p[t] increasing) here to plot both p[t] converging 
to m = 2, and U[t] converging to Ue = 5.

Note that in the short run, unemployment U[t] can be initially increasing (if U[0] < 5) or 
initially decreasing (if U[0] > 5), whether or not p[t] is initially decreasing (p[0] > 2) or 
initially increasing (p[0] < 2).

Check that b1^2 < 4*b2 (with these parameter choices) using R(g, j, kbeta).

R(0.9, 0.9, 5), numer;

The negative sign says b1^2 < 4*b2 and we have complex roots (as we saw in the complex
roots example for p[t]).

IUparam;

uIUparam;

Ucase2 : [β = 2, k = 5/2, g = 9/10, j = 9/10,  m = 2];

at (uIUparam, Ucase2);

ratsubst(11, α - T, %);

croots(-241/600, 91/600);

at (Ue, Ucase2);

Ue2 : ratsubst (11, α - T, %), numer;

UtL1_complex : ytlist (-241/600, 91/600,81/20,1, 1.5, 5);

UtL2_complex : ytlist (-241/600, 91/600, 81/20,7.5, 7, 5);

Plot of UtL1_complex (increases) and ptL1_complex (decreases)

wxdraw2d (xlabel = "t",  yrange = [0,10], xrange = [0, 6],
    title = "short-run inflation decreasing, unemployment increasing",
    key = "U[t]", points (tL, UtL1_complex), color = magenta, key = "p[t]",
    points (tL, ptL1_complex), color = black, line_width = 1, key = "Ue = 5.4",
    explicit (Ue2, t, 0, 6), color = brown, key = "m = 2", explicit (2, t, 0, 6))$

Plot of UtL2_complex (decreases) and ptL1_complex (decreases)

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
    title = "short-run inflation decreasing, unemployment decreasing",
    key = "U[t]", points (tL, UtL2_complex), color = magenta, key = "p[t]",
    points (tL, ptL1_complex), color = black, line_width = 1, key = "Ue = 5.4",
    explicit (Ue2, t, 0, 6), color = brown, key = "m = 2", explicit (2, t, 0, 6))$

Plot of UtL1_complex (increases) and ptL2_complex (increases)

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
    title = "short-run inflation increasing, unemployment increasing",
    key = "U[t]", points (tL, UtL1_complex), color = magenta, key = "p[t]",
    points (tL, ptL2_complex), color = black, line_width = 1, key = "Ue = 5.4",
    explicit (Ue2, t, 0, 6), color = brown, key = "m = 2", explicit (2, t, 0, 6))$

Plot of UtL2_complex (decreases) and ptL2_complex (increases)

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
    title = "short-run inflation increasing, unemployment decreasing",
    key = "U[t]", points (tL, UtL2_complex), color = magenta, key = "p[t]",
    points (tL, ptL2_complex), color = black, line_width = 1, key = "Ue = 5.4",
    explicit (Ue2, t, 0, 6), color = brown, key = "m = 2", explicit (2, t, 0, 6))$

## Variable Term and Higher Order Difference Equations

In this section we follow Chiang and Wainwright, Ch. 18, Sec. 4.

### 2nd Order Linear Difference Equations with a Variable Term

"A variable term" means that 'a' depends on the time t, in the equations
y[t] + b1*y[t-1] + b2*y[t-2]  = a, or y[t] + b1*y[t-1] + b2*y[t-2] + b3*y[t-3]  = a, etc.

### y[t] + y[t-1] -3*y[t-2] = 7^t or y[t+2] + y[t+1] - 3*y[t] = 7^(t+2)

The quickest route to a numerical solution is to use our Maxima function 
    ytlist (b1, b2, a, y0, y1, tmax). We assume y[0] = 2, y[1] = 1 here.
The use of ytlist with this syntax depends on starting with a difference equation
first being written in the form:  y[t] + b1*y[t-1] + b2*y[t-2] = a.
Starting with y[t] + y[t-1] - 3*y[t-1] = 7^t, if we let t --> t + 2 in every term, we get
 y[t+2] + y[t+1] - 3*y[t] = 7^(t+2). However, using ytlist requires we use the 'a' 
which matches the form y[t] + y[t-1] - 3*y[t-1] = 7^t.

tL : makelist (tt, tt, 0, 6);
ytL : ytlist (1,-3, 7^t, 2, 1, 6);

print ("y[t] + y[t-1] - 3*y[t-1] = 7^t  or  y[t+2] + y[t+1] - 3*y[t] = 7^(t+2)")$
wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0, 3],  yrange = [-10, 300],
    title = "y0 = 2, y1 = 1 solution",
    points (tL, ytL), color = black, line_width = 1, explicit (0, t, 0, 3))$

A longer route to a solution, which has the advantage of revealing the particular solution,
is to use the Maxima function solve_rec (difference-eqn, y[t]) defined in the contributed
code file solve_rec.mac, which we must first load. When using solve_rec we can use
either form of the difference equation.

load(solve_rec);

gsoln : solve_rec (y[t] + y[t-1] - 3*y[t-2] = 7^t, y[t]);

We see that the particular solution yp[t] = 7^(t+2)/53.

The right hand side of gsoln is our indefinite solution, which contains two constants, 
%k[1] and %k[2].

yindef : rhs(gsoln);

grind (yindef);

Define a Maxima function ys(t) by placing two single quotes (') in front of yindef.

ys(t) := ''yindef;

We look for a solution for which y[0] = 2 and y[1] = 1.

solns : solve ([ys(0) = 2, ys(1) = 1 ]);

ydef : at (yindef, solns[1]);

yd (t) := ''ydef;

tL;

ytL : map ('yd, tL), numer;

Note that this list is arrived at much more quickly using ytlist(b1,b2, a, y0, y1, tmax),
as we showed above.

print ("y[t] + y[t-1] - 3*y[t-1] = 7^t  or  y[t+2] + y[t+1] - 3*y[t] = 7^(t+2)")$
wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0, 3],  yrange = [-10, 300],
    title = "y0 = 2, y1 = 1 solution",
    points (tL, ytL), color = black, line_width = 1, explicit (0, t, 0, 3))$

### y[t+2] + y[t+1] - 3*y[t] = 7^t, y[t] + y[t-1] - 3*y[t-2] = 7^(t-2)

Chiang and Wainwright, Sec. 18.4, start with a difference equation expressed in the form
  y[t+2] + y[t+1] - 3*y[t] = 7^t. If we use solve_rec (difference-eqn, y[t]) with this form, we
get:

gsoln : solve_rec (y[t + 2] + y[t + 1] - 3*y[t] = 7^t, y[t]);

yindef : rhs (gsoln);

grind(yindef)$

ys(t) := ''yindef;

solns : solve ([ys(0) = 2, ys(1) = 1]);

ydef : at (yindef, solns[1]);

yd(t) := ''ydef;

tL;

ytL : map ('yd, tL), numer;

To instead use ytlist (b1, b2, a, y0,y1, tmax), we need to translate the difference equation
into the form: y[t] + y[t-1] - 3*y[t] = 7^(t-2), in which b1 = 1, b2 = -3, 
a = 7^(t-2). 

ytL : ytlist (1, -3, 7^(t-2), 2, 1, 6);

We see we get the same solution as gotten the long way using solve_rec.
We make a plot of our solution.

tL;

print ("y[t+2] + y[t+1] - 3*y[t] = 7^t or y[t] + y[t-1] -3*y[t] = 7^(t-2)")$
wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0, 5],  yrange = [-10, 300],
    title = " y[0] = 2, y[1] = 1",
    points (tL, ytL), color = black, line_width = 1, explicit (0, t, 0, 5))$

### y[t+2] + 5*y[t+1] + 2*y[t] = t^2 or y[t] + 5*y[t-1] + 2*y[t-2] = (t-2)^2

tL;

In order to use ytlist, we use the form y[t] + 5*y[t-1] + 2*y[t-2] = (t-2)^2 to identify b1, b2, 
and 'a'.

ytL : ytlist (5, 2, (t-2)^2, 2, 1, 6);

print ("y[t+2] + 5*y[t+1] + 2*y[t] = t^2 or y[t] + 5*y[t-1] + 2*y[t-2] = (t-2)^2")$
wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0, 4],  yrange = [-200, 50],
    title = " y[0] = 2, y[1] = 1",
    points (tL, ytL), color = black, line_width = 1, explicit (0, t, 0, 4))$

To check our methods, we use the the translated form: y[t+2] + 5*y[t+1] + 2*y[t] = t^2 
with solve_rec.

gsoln : solve_rec (y[t + 2] + 5*y[t + 1] + 2*y[t] = t^2, y[t]);

yindef : rhs (gsoln);

ys(t) := ''yindef;

solns : solve ([ys(0) = 2, ys(1) = 1]);

ydef : at (yindef, solns[1]);

yd(t) := ''ydef;

ytL2 : map ('yd, tL), numer;

This is the same list of values as we got by using ytlist above, so the plot
will look the same.

### Higher-Order Linear Difference Equations

### y[t+3] - (7/8)*y[t+2] + (1/8)*y[t+1]  + y[t]/32 = 9

gsoln : solve_rec (y[t+3] - 7*y[t+2]/8 + y[t+1]/8 + y[t]/32 = 9, y[t]);

yindef : rhs (gsoln);

ys(t) := ''yindef;

solns : solve ([ys(0) = 2, ys(1) = 3, ys(2) = 4]);

ydef : at (yindef, solns[1]);

yd(t) := ''ydef;

tL;

ytL : map ('yd, tL), numer;

print ("y[t+3] - (7/8)*y[t+2] + (1/8)*y[t+1] + y[t]/32 = 9 ")$
wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0, 6],  yrange = [0, 30],
    title = " y[0] = 2, y[1] = 3, y[2] = 4", points (tL, ytL))$