# Solving Simultaneous Difference Eqns.

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

load ("Econ2.mac");

## Preface

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

In Dowling19B.wxmx, we use Maxima to discuss the matrix solution of a set of n linear
first order difference equations. We end with a reconsideration of the inflation and 
unemployment model discussed in Dowling18C.wxmx, following Chiang and 
Wainwright's Ch. 19, Sec. 4.

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

Ted Woollett
https://home.csulb.edu/~woollett/
woollett@charter.net
Mar. 9, 2022

## References

Introduction to Mathematical Economics, 3rd ed., Edward T. Dowling, 2012, Schaum's 
Outline Series, McGraw-Hill.

Fundamental Methods of Mathematical Economics, Alpha C. Chiang and Kevin Wainwright, 
4th ed., 2005, McGraw-Hill

## Matrix Solution of Y[t] = A . Y[t-1] + B

We have provided a brief review of Maxima matrix functions in Dowling19A.wxmx, 
which you should consult first.

If we let I represent the n x n unit matrix, each component of Ys[t] will be the matrix sum
of  the particular solution, which applies to a hypothetical time in the future when
 Y[t] = Y[t - 1] and the matrix equation becomes Yp = A . Yp + B, or (I - A) . Yp = B, 
whence
      Yp = invert(I - A) . B,
and Yc, the solution of the complementary equation Y[t] = A . Y[t - 1], which solution will
have  the form (if we assume distinct real eigenvalues)
      Yc = sum (k[j]*V[j]*r[j]^t, j, 1, n), 
where k[j] are adjustable constants, and r[j] and V[j] are respectively eigenvalues and 
corresponding eigenvectors must satisfy the equations (j = 1 thru n)
      1.   determinant (A - r[j]*I)  = 0
      2.   A . V[j] = r[j]*V[j]
The constants k[j] in Yc are determined by using prescribed initial values Y[0].
Again assuming distinct real eigenvalues, the dominant root is the root with the largest
absolute value. For convergence we need the magnitude of the dominant root to be 
less than 1.

### Proof of the Yc Solution Form for n = 2

For n = 2 we assert Yc[t] has the form
   Yc[t] = k1*r1^t*V1 + k2*r2^t*V2.
We need to show that Yc[t] = A . Yc[t-1].
Now A . Yc[t-1] = k1*r1^(t-1)* (A . V1) + k2*r2^(t-1)* (A . V2)
                        = k1*r1^(t-1)*r1*V1 + k2*r2^(t-1)*r2*V2
                        = k1*r1^t*V1 + k2*r2^t*V2 = Yc[t]. 

### Example 5:  n = 2 Distinct Real Roots

Solve the following system of first-order linear difference equations:
   x[t] = - 4*x[t-1] + y[t-1] + 12, x[0] = 16,
   y[t] =   2*x[t-1] - 3*y[t-1] + 6,  y[0] = 8.

### Using Ytlist (A, B, Y0, tmax)

The quickest way to produce a visual look at the solution of a set of first order 
difference equations written in the matrix form
   Y[t] = A . Y[t-1] + B
is to use the Maxima function (defined in Econ2.mac) Ytlist (A, B, Y0, tmax).

Ytlist (A, B, Y0, tmax) produces a column vector whose jth component is the list of values
         [yj[0], yj[1], yj[2],...,yj[tmax] ],
using directly the first order difference equation Y[t] = A . Y[t-1] + B repeatedly, 
starting with Y[1] = A . Y0 + B, without using the eigenvalues or eigenvectors of A.

We revert to the notation Ye (instead of Yp) for the particular solution, since its physical
meaning is the (hypothetical) long time equilibrium value of the system vector Y[t].

A : matrix ( [-4, 1], [2, -3] );
eigenvalues (A);
B : cvec ( [12, 6]);
Y0 : cvec ([16,8] );
Ye : invert (ident(2) - A) . B;
Yt : Ytlist (A, B, Y0, 5);

grind (Yt)$

The dominant root is -5 whose magnitude is greater than 1, hence a divergent system.

Let xs be the list of values of x[t] from t = 0 to t = 5. Likewise, ys is the list of
values of y[t].

lme is our alias (defined in Econ2.mac) for list_matrix_entries.

[xs, ys] : lme (Yt);

xs;

The list of times tL should be the same length as xs or ys.

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

We are only going to display the values of x[t] and y[t] out to the time  t = 3 on this plot.

wxdraw2d (xlabel = "t", xrange = [0, 3.5], yrange = [-1e3, 1e3],
    key_pos = top_left,  key = "x[t]", points (tL, xs), color = red, key = "y[t]",
    points (tL, ys), color = black, line_width = 1, key = "",
    explicit (0, t, 0,3.5))$

We see the expected divergent behavior of this system.

### Using DE1S (A, B, t, Y0)

The Maxima function DE1S (A, B, t, Y0), defined in Econ2.mac, uses the eigenvalues 
and corresponding eigenvectors of the square matrix A to find a solution to the 
complementary equation Yc[t] = A . Yc[t-1], and that solution has the form
  Yc = k1*r1^t*V1 + k2*r2^t*V2 (k1 and k2 being arbitrary constants)
(if the eigenvalues are real and distinct) for n = 2, where V1 and V2 are the corresponding
eigenvectors: A . V1 = r1*V1, A . V2 = r2*V2.

The indefinite solution is then Yc + Ye, where Ye is the hypothetical long time equilibrium
state in which Y[t] = Y[t-1], so Ye = A . Ye + B, and hence
   Ye = invert (I - A) . B,
where I = ident(n) is the n x n unit matrix.

DE1S (A, B, t, Y0) or DE1S (A, B, t)  returns Y[t], the matrix column vector solution
      of a set of n first order difference equations expressed in matrix form as
        Y[t] = A . Y[t-1] + B.
     Y0 is the optional matrix column vector of initial values.   
     Y[t] has with n components if A is a square n x n matrix, and is returned in 
         definite form if the fourth arg Y0 is included, else in terms of n constants _k[j].

DE1S (A, B, t);

Ys : DE1S (A, B, t, Y0);

Because we have included Y0, the definite form is returned which is written using the
eigenvalues -2 and -5 (call them r1 and r2 for this example)

[xs, ys] : lme (Ys);

xs;

Here, xs and ys are not lists, rather they are each a Maxima expression depending on 
the parameter t.

at ( [xs, ys], t = 0 );

txL : makelist ([tt, at (xs, t = tt)], tt, 0, 6);

tyL : makelist ([tt, at (ys, t = tt)], tt, 0, 6);

The divergence is clear from these values of x[t] and y[t].

We again only look out to t = 3.

wxdraw2d (xlabel = "t", xrange = [0, 3.5], yrange = [-1e3, 1e3],
    key_pos = top_left,  key = "x[t]", points (txL), color = red, key = "y[t]",
    points (tyL), color = black, line_width = 1, key = "",
    explicit (0, t, 0,3.5))$

This plot reproduces the plot made above using Ytlist.

### Matrix Details of DE1S Code

The solution Yc of the complementary equation
  Yc[t] = A . Yc[t-1]
depends on the eigenvalues and corresponding eigenvectors of A.

[eval, evec] : eigenvectors (A);

evec;

evec[1][1];

evec[2][1];

eval;

eval[1];

eval[1][1];

The eigenvalues of A are real and distinct and each have magnitudes greater than 1,
which leads to a divergent system time path.

r1 : -5;
V1 : cvec ([1, -1]);
A . V1;
r2 : -2;
V2 : cvec ([1,2]);
A . V2;

Here is the complementary equation solution Yc.

Yc : k1*r1^t*V1 + k2*r2^t*V2;

Adding up the matrix column vectors Yc and Ye (defined way up above) defines the 
"indefinite solution" Yindef:

Yindef : Yc + Ye;

The conponents of Y0 are x[0] and y[0]. (We defined this way above).

Y0 ;

We use our Maxima function colVecSolve (col1, col2), defined in Econ2.mac, to require
Yindef = Y0 when t = 0.

solns : colVecSolve (at (Yindef, t = 0), Y0);

Let Ys be our solution Y[t] taking into account our initial conditions.

Ys : at (Yindef, solns);

grind (Ys)$

Let xs and ys be the components of Ys. lme is our alias for list_matrix_entries and is
defined in Econ2.mac.

[xs, ys] : lme (Ys);

at ( [xs, ys], t = 0 );

txL : makelist ([tt, at (xs, t = tt)],tt, 0, 6);

tyL : makelist ([tt, at (ys, t = tt)],tt, 0, 6);

The divergence is clear from these values of x[t] and y[t].

We can use points (txL) , for example, to show the values of x[t] on a plot
of x[t] versus t, for discrete values of t = 0, 1, 2, ...

wxdraw2d (xlabel = "t", xrange = [0, 3.5], yrange = [-1e3, 1e3],
    key_pos = top_left,  key = "x[t]", points (txL), color = red, key = "y[t]",
    points (tyL), color = black, line_width = 1, key = "",
    explicit (0, t, 0,3.5))$

## Example with n = 2 Complex Roots

kill (a,b,c,d)$
A : matrix ([a,b],[c,d]);
eigenvalues (A);

To have complex eigenvalues we need the arg of the square root equal to a negative
number, or  (d - a)^2 + 4*b*c < 0, or (d - a)^2 < - 4*b*c, so we need b # 0 and c # 0 and
b and c to have opposite signs. If we take c = - e*b, then we need e > 0 and b # 0.
Example: take c = - b, then we need (d - a)^2 < 4*b^2. If we take b = 2, we need
(d - a)^2 <  16, which is satisfied if we take a = 2, d = 4, d - a = 2, (d - a)^2 = 4 < 16.

A : matrix ([2, 2], [-2, 4]);
[r1, r2] : eigenvalues (A)[1];

map ('abs, [r1, r2]);

float(%);

Since we can write, using |r1| = sqrt (x1^2 + y1^2) if r1 = x1 + %i*y1,
  r1 = |r1|*exp (%i*th1) =  |r1|* (sin (th1) + %i*cos (th1)), and when we raise r1 to the
integral power t, we get 
    r1^t =  |r1|^t*exp (%i*t*th1)
and since the part of Yc proportional to r1^t is
        k1*r1^t*V1 = k1*|r1|^t*exp (%i*t*th1)*V1,
we need |r1| < 1 for this part of Yc to converge. Since r1 and r2 will be complex conjugates
of each other, we need the absolute value of either root to be less than 1 for convergence.

Maxima's abs or cabs functions will return the absolute value of a number,
  where if r = x + %i*y, |r| = sqrt (x^2 + y^2).

abs (3 - %i*sqrt(3));

float(%);

Since this is greater than 1, the solution is divergent.

We take the same B and Y0.

B;
Y0;

Ye : invert (ident(2) - A) . B, numer;

### Using Ytlist (A, B, Y0, tmax)

Yt : Ytlist (A, B, Y0, 6);

[xs, ys] : lme (Yt);

xs;

txL : makelist ([t, xs[t+1]], t, 0, 4);

tyL : makelist ([t, ys[t+1]], t, 0, 4);

wxdraw2d (xlabel = "t", xrange = [0, 4.5], yrange = [-3e3, 300],
    key_pos = bottom_left,  key = "x[t]", points (txL),  color = red, key = "y[t]",
    points (tyL), color = black, line_width = 1, key = "",
    explicit (0, t, 0, 4.5))$

### Using DE1S (A, B, t, Y0)

Ys : DE1S (A, B, t, Y0);

Because we have included Y0, the definite form is returned which is written using the
eigenvalues (call them r1 and r2 for this example)

[xs, ys] : lme (Ys);

xs;

at ( [xs, ys], t = 0 ), ratsimp;

txL : makelist ( [tt, ratsimp ( at (xs, t = tt))], tt, 0, 6);

tyL : makelist ( [tt, ratsimp ( at (ys, t = tt))], tt, 0, 6);

wxdraw2d (xlabel = "t", xrange = [0, 4.5], yrange = [-3e3, 300],
    key_pos = bottom_left,  key = "x[t]", points (txL),  color = red, key = "y[t]",
    points (tyL), color = black, line_width = 1, key = "",
    explicit (0, t, 0, 4.5))$

### n = 2 Complex Root Case Convergence

Begin with the form of the matrix A developed above to produce complex roots.
     A : matrix ([a, b], [-e*b, d] )
with e > 0 and b # 0. What must be true to now have convergence?

A : matrix ( [a, b], [- e*b, d]);
[r1, r2 ] : eigenvalues (A)[1];

Assume 4*e*b^2 > (d - a)^2 to have complex roots, and rewrite r2 as
  r2 = (1/2)*(d + a) + %i*(1/2)*sqrt(4*e*b^2  - (d-a)^2) = x2 + %i*y2.
Then |r2|^2 = x2^2 + y2^2:

expand ((d + a)^2/4 + (4*e*b^2 - (d - a)^2)/4);

Thus |r2| = sqrt (a*d + e*b^2), and e > 0, b^2 > 0.
Thus for convergence we need (a*d + e*b^2) < 1.

With a = b = d = 1/2 and e = 1 we then expect convergence.

A : matrix ( [1/2, 1/2], [- 1/2, 1/2]);
[r1, r2] : eigenvalues (A)[1];

map ('abs, [r1, r2]);

float (%);

Since the magnitude of the complex roots is ~ 0.71 < 1, we expect convergence.

Again we assume the same B and Y0.

B;
Y0;

Ye : invert (ident(2) - A) . B, numer;

Yt : Ytlist (A, B, Y0, 8);

[xs, ys] : lme (Yt);

xs,numer;

ys,numer;

txL : makelist ([t, xs[t+1]], t, 0, 8);

tyL : makelist ([t, ys[t+1]], t, 0, 8);

tmax : 8.5$
wxdraw2d (xlabel = "t", xrange = [0, tmax], yrange = [-30, 30],
    key_pos = bottom_left,  key = "x[t]", points (txL),  color = red, key = "y[t]",
    points (tyL), line_width = 1,  color = blue, key = "xe = 18", explicit (18, t, 0, tmax),
    color = red, key = "ye = - 6", explicit (-6, t, 0, tmax), color = black,  key = "",
    explicit (0, t, 0, tmax))$

We see convergence (with oscillation) in this complex root example.

## n = 2 Repeated Real Roots Case

kill (a, b, c, d)$
A : matrix ([ a, b], [c, d] );
[r1, r2 ] : eigenvalues (A)[1];

What do we need to get r1 = r2?

r1mr2 : expand (2*r1 - 2*r2);

r1mr2 : r1mr2/(-2);

r1mr2 : r1mr2^2;

So for two equal roots we need (d-a)^2 + 4*b*c = 0, so the sign of c must be opposite
to the sign of b, so replace c by -e  so A = matrix ([a,b],[-e,d]), and then we need
4*b*e = (d-a)^2, or A = matrix ([a, b], [-(d-a)^2/(4*b), d]).

A : matrix ([a, b], [- (d-a)^2/(4*b), d]);
eigenvalues (A);

So this is the case we want, one root  r = (d+a)/2 with multiplicity 2, and we can 
assume a solution of the form
    Yc[t] =  (k1*r^t + k2*r^(t + 1) )*V, where A . V = r*V, 
since Yc[t-1] = (k1*r^(t-1) + k2*r^t) * V, and then  A . Yc[t-1] = Yc[t].

Here is an example which results in a repeated real root r = 1/2, with |r| < 1, and
a convergent system. We continue to use the same matrix column vectors for B and
Y0 as above. We find a brute force solution by using Ytlist(A,B,Y0, tmax).

A : matrix ([1/2, 1/2], [ 0, 1/2]);
eigenvalues (A);

B;
Y0;

Ye : invert (ident(2) - A) . B, numer;

Yt : Ytlist (A, B, Y0, 8), numer;

[xs, ys] : lme (Yt);

xs;

txL : makelist ([t, xs[t+1]], t, 0, 8);

tyL : makelist ([t, ys[t+1]], t, 0, 8);

tmax : 8.5$
wxdraw2d (xlabel = "t", xrange = [0, tmax], yrange = [0, 70],
    key_pos = top_left,  key = "x[t]", points (txL),  color = red, key = "y[t]",
    points (tyL), color = blue, line_width = 1, key = "xe = 36",
    explicit (36, t, 0, tmax), color = red, key =  "ye = 12",
    explicit (12, t, 0, tmax))$

DE1S (A, B, t, Y0);

## Matrix Solution of  A1 . Y[t] = A2 . Y[t-1] + B

DE2S (_A1, _A2, _B, _t, _Y0) :=
block (
   if det (_A1) = 0 then return (" DE2S method depends on det (A1) not being zero"),
   DE1S (invert (_A1) . _A2, invert (_A1) . _B, _t, _Y0))$

### Dowling Ex. 7 & 8

Solve the set of first order difference equations
  x[t] = 4*x[t-1] - 2*y[t-1] + y[t] - 10, x[0] = 20,
  y[t] = 3*x[t-1] + 6*y[t-1] - 4, y[0] = 3.

Writing this set in matrix form
  A1 . Y[t] = A2 . Y[t-1] + B,
with the initial values in Y0,

A1 : matrix ([1,-1], [0, 1]);
A2 : matrix ([4, -2], [3, 6]);
B : cvec ([-10, -4]);
Y0 : cvec ([20, 3]);

### Using Ytlist (D, E, Y0, tmax)

We can convert the given matrix equation to the form
Y[t] = D . Y[t-1] + E, with Y0 unchanged, by multiplying the given matrix equation
from the left by invert(A1), and using invert (A1) . A1 = I = unit matrix and then using
I . Y[t] = Y[t],  and defining D and E matrices:

D : invert (A1) . A2;
E : invert (A1) . B;

eigenvalues (D);

Real eigenvalues with magnitudes much greater than 1 implies a very divergent system.

Ye : invert (ident(2) - D) . E;
Ys : Ytlist (D, E, Y0, 6);

[xs, ys] : lme (Ys);

tL : makelist (tt, tt, 0, 6);

tmax : 2.5$
wxdraw2d (xlabel = "t", xrange = [0, tmax], yrange = [0, 1.5e3],
    key_pos = top_left,  key = "x[t]", points (tL, xs), color = red, key = "y[t]",
    points (tL, ys) )$

### Using DE2S (A1, A2, B, t, Y0)

DE2S (A1,A2,B, t,Y0);

## Inflation-Unemployment Set of 1st Order Difference Eqns

We follow Chaing and Wainwright's Ch. 19 treatment of this topic. In Dowling18C.wxmx
we arrived at a set of discrete time equations, which we conbined into one second order
difference equation. Here we want to work with a set of two first order difference
equations.

With p the actual inflation rate (%), π the expected inflation rate (%), and U the 
unemployment rate (%), μ the fixed rate of increase of money [μ = (1/M)*dM/dt]
set as a monetary policy decision by the FED, in 18C  we translated the continuous time
equations into a discrete time version (here we have replaced  m --> μ):

p[t] = α - T - β*U[t] + g*π[t], ( 0 < g <= 1), (α, β > 0),

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

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

We eliminate p[t] and p[t+1] in the second and third equations:
π[t+1] - π[t] = j*(α - T - β*U[t] + g*π[t])  - j*π[t],                                               (1)

U[t+1] - U[t] = -k*μ + k* (α - T - β*U[t+1] + g*π[t+1])                                        (2)

Our methods in this chapter 19B are based on the two matrix forms
  Y[t] = A . Y[t-1] + B, or
  A1. Y[t] = A2 . Y[t-1] + B.

So we replace t --> t - 1 in both (1) and (2) to get

π[t] = (1 - j*(1-g))*π[t-1] - j*β*U[t-1] + j*(α - T),                                               (3)

- k*g*π[t] + (1 + β*k)*U[t] = U[t-1] + k*(α - T - μ).                                            (4)

We can then write the pair of equations (3) and (4) in matrix form as

  A1 . Y[t] = A2 . Y[t-1] + B                                                                            (5)

with the following assignments:

A1 : matrix ([1, 0], [-k*g, 1 + β*k]);
A2 : matrix ([1 - j*(1 - g), - j*β], [0, 1]);
B : cvec ([j*(α  - T), k*(α - T - μ)]);

In Dowling18C.wxmx we considered two cases, which we will call case1 and case2 here.
We defined a list of parameter replacements, but replace α - T = 11 separately using
ratsubst.

### Case 1: Real Eigenvalues

case1 : [β = 2, k = 5/2, g = 1/2, j = 1/3, μ = 2];

A1_case1 : at (A1, case1);
A2_case1 : at (A2, case1);
B_case1 : at (B, case1);
B_case1 : ratsubst (11, α - T, B_case1);

### Using Ytlist (D, E, Y0, tmax)

One approach is to multiply the matrix equation from the left by invert (A1), arriving at the
equation Y[t] = D . Y[t-1] + E, (using invert(A1) . A1 = ident(2) = 2 x 2 unit matrix, and
if I is the appropriate unit matix, I . A = A).

D : invert (A1_case1) . A2_case1;
E : invert (A1_case1) . B_case1;

eigenvalues (D), numer;

We have real eigenvalues with magnitudes less than 1, so we expect eventual convergence.

Ye : invert (ident(2) - D) . E;

So we expect long time convergence with πe = 2 and Ue = 5.

For our first plot, we assume π[0] = 4 and U[0] = 1.

Y0 : cvec ([4, 1]);

Ys : Ytlist (D, E, Y0, 6), numer;

[πs, Us] : lme (Ys);

We need the list of times (tL) to be the same length as πs and Us.

tL : makelist (tt, tt, 0, 6);

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
 title = "g = 1/2, j = 1/3, β = 2, k = 5/2, μ=2, α - T = 11",
 key = "U[t]", points (tL, Us), color = magenta, key = "pi[t]",
 points (tL, πs), color = black, line_width = 1, key = "Ue = 5",
 explicit (5, t, 0, 6), color = brown, key = "μ = 2", explicit (2, t, 0, 6))$

### Using DE1S (D, E, t, Y0)

Ys : DE1S (D, E, t, Y0);

[πs, Us] : lme (Ys),numer;

tπL : makelist ([tt, at (πs, t = tt)], tt, 0, 6);

tUL : makelist ([tt, at (Us, t = tt)], tt, 0, 6);

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
 title = "g = 1/2, j = 1/3, β = 2, k = 5/2, μ=2, α - T = 11",
 key = "U[t]", points (tUL), color = magenta, key = "pi[t]",
 points (tπL), color = black, line_width = 1, key = "Ue = 5",
 explicit (5, t, 0, 6), color = brown, key = "μ = 2", explicit (2, t, 0, 6))$

### Using DE2S (A1, A2, B, t, Y0)

Ys : DE2S (A1_case1, A2_case1, B_case1, t, Y0);

[πs, Us] : lme (Ys),numer;

tπL : makelist ([tt, at (πs, t = tt)], tt, 0, 6);

tUL : makelist ([tt, at (Us, t = tt)], tt, 0, 6);

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
 title = "g = 1/2, j = 1/3, β = 2, k = 5/2, μ=2, α - T = 11",
 key = "U[t]", points (tUL), color = magenta, key = "pi[t]",
 points (tπL), color = black, line_width = 1, key = "Ue = 5",
 explicit (5, t, 0, 6), color = brown, key = "μ = 2", explicit (2, t, 0, 6))$

### Changing Y0

We continue with the same case1 parameters but change the initial conditions, using
Ytlist for simplicity, first with π[0] = 4, U[0] = 7.

Y0 : cvec ([4, 7]);

Ys : Ytlist (D, E, Y0, 6), numer;

[πs, Us] : lme (Ys);

We need the list of times (tL) to be the same length as πs and Us.

tL : makelist (tt, tt, 0, 6);

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
 title = "g = 1/2, j = 1/3, β = 2, k = 5/2, μ=2, α - T = 11",
 key = "U[t]", points (tL, Us), color = magenta, key = "pi[t]",
 points (tL, πs), color = black, line_width = 1, key = "Ue = 5",
 explicit (5, t, 0, 6), color = brown, key = "μ = 2", explicit (2, t, 0, 6))$

We see both inflation and unemployment decreasing smoothly to their long term values.

We next start with U[0] = 7, π[0] = 0.5.

Y0 : cvec ([0.5, 7]);

Ys : Ytlist (D, E, Y0, 6), numer;

[πs, Us] : lme (Ys);

We need the list of times (tL) to be the same length as πs and Us.

tL : makelist (tt, tt, 0, 6);

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
 title = "g = 1/2, j = 1/3, β = 2, k = 5/2, μ=2, α - T = 11",
 key = "U[t]", points (tL, Us), color = magenta, key = "pi[t]",
 points (tL, πs), color = black, line_width = 1, key = "Ue = 5",
 explicit (5, t, 0, 6), color = brown, key = "μ = 2", explicit (2, t, 0, 6))$

We see some oscillation of unemployment around Ue = 5 but eventual convergence of
both U[t] and π[t].

### Case 2: Complex Eigenvalues

case2 : [β = 2, k = 5/2, g = 0.9, j = 0.9, μ = 2];

A1_case2 : at (A1, case2);
A2_case2 : at (A2, case2);
B_case2 : at (B, case2);
B_case2 : ratsubst (11, α - T, B_case2);

Redefine D and E for case2.

D : invert (A1_case2) . A2_case2;
E : invert (A1_case2) . B_case2;

[r1, r2] : eigenvalues (D)[1];

abs(r1), numer;

We have complex eigenvalues rather than real as in case1. r1 and r2 are complex
conjugates of each other. For convergence we need the absolute value of either root
to be less than 1. Since |r1| = 0.4 < 1, we have a convergent case.

Ye : invert (ident(2) - D) . E;

We expect convergence to the values πe = 2 = μ and Ue = 5.4.

For our first plot with case2, we assume π[0] = 4 and U[0] = 1.

Y0 : cvec ([4, 1]);

Ys : Ytlist (D, E, Y0, 6), numer;

[πs, Us] : lme (Ys);

We need the list of times (tL) to be the same length as πs and Us.

tL : makelist (tt, tt, 0, 6);

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
 title = "g = 1/2, j = 1/3, β = 2, k = 5/2, μ=2, α - T = 11",
 key = "U[t]", points (tL, Us), color = magenta, key = "pi[t]",
 points (tL, πs), color = black, line_width = 1, key = "Ue = 5.4",
 explicit (5.4, t, 0, 6), color = brown, key = "μ = 2", explicit (2, t, 0, 6))$

Both π[t] and U[t] undergo one large oscillation before settling down.

Changing Y0 in case 2, with π[0] = 4, U[0] = 7.

Y0 : cvec ([4, 7]);

Ys : Ytlist (D, E, Y0, 6), numer;

[πs, Us] : lme (Ys);

wxdraw2d (xlabel = "t", yrange = [0,10], xrange = [0, 6],
 title = "g = 1/2, j = 1/3, β = 2, k = 5/2, μ=2, α - T = 11",
 key = "U[t]", points (tL, Us), color = magenta, key = "pi[t]",
 points (tL, πs), color = black, line_width = 1, key = "Ue = 5.4",
 explicit (5.4, t, 0, 6), color = brown, key = "μ = 2", explicit (2, t, 0, 6))$

After one oscillation, both π[t] and U[t] converge to their respective long term values.