In [ ]:
#import plt
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(xs[:,0], label='Employed')
ax.plot(xs[:,1], label='Unemployed')
ax.legend()
Out[ ]:
<matplotlib.legend.Legend at 0x7feb41e5d220>
Using only the qr function and matrix multiplication, use the QR algorithm to find the eigenvalue of the following matrix:
$$\begin{bmatrix} 1,3,4\\3,1,0\\4,0,1 \end{bmatrix}$$Use a while loop to automatically loop through the iterations. What would be an appropriate stopping/exit decision for the loop?
import numpy as np
from numpy.linalg import qr
A = np.array([[1,3,4],[3,1,0],[4,0,1]])
q, r = qr(A)
#iterate this
A_new = np.dot(r, q)
q,r = qr(A_new)
print('QR:', A_new)
QR: [[ 2.92307692 -4.44064833 -1.2579418 ] [-4.44064833 -0.78021978 -0.50429863] [-1.2579418 -0.50429863 0.85714286]]
Until convergence towards:
$$\begin{bmatrix} 6, 0, 0 \\0, -4, 0 \\0, 0, 1 \end{bmatrix}$$A_new[1,0]
-4.440648327160281
#or instead
e = np.abs(A_new[1,0])
while e>.0001:
A_new = np.dot(r, q)
q,r = qr(A_new)
e = np.abs(A_new[1,0])
print(A_new)
[[ 6.00000000e+00 -7.82264258e-05 -7.34780335e-17] [-7.82264258e-05 -4.00000000e+00 -7.51391859e-17] [-2.55880669e-22 -1.63551298e-17 1.00000000e+00]]
This is based on a lecture from the the datascience quantecon course
We consider an economy where a certain percentage, $\alpha$, of people in any given year lose their jobs and a certain percentage of people, $\theta$ find a job.
We define a transition matrix
$$A = \begin{bmatrix} 1-\alpha & \alpha \\ \theta & 1-\theta \end{bmatrix}$$Let's say that initially (period 0) there are 1,000,000 employed people and 150,000 unemployed people: we can put this in a vector:
$$x0 = \begin{bmatrix} 1.000.000 & 150.000 \end{bmatrix}$$a.) Assuming that $\alpha=0.04$ and $\theta=0.10$ Use matrix multipication to find how many employed and unemployed people there will be in period 1, $x1$. Explain the calculation.
b.) What about in periods $x2$, $x3$, ... and so on? Do the values converge to some stable numbers?
c.) What are the eigenvectors of A? How does this relate to b?
# a.) we first define our transition matrix and our vector of initial values:
import numpy as np
alpha = .04
theta = .10
A = np.array([[1-alpha, alpha],[theta, 1-theta]])
x0 = np.array([1000000,150000]).T
x0
array([1000000, 150000])
np.dot(A.T, x0)
array([975000., 175000.])
When we do this calculation we are creating the vector:
$$ x_1 = \begin{bmatrix} 1.000,000 * 1-\alpha + 150.000* \theta \\1.000,000 * 1-\alpha + 150.000* 1-\theta \end{bmatrix} = \begin{bmatrix} 975.000 \\ 175.000 \end{bmatrix}$$The first row of the new vector shows that $1-\alpha = 96\%$ of the 1 million that initially had a job still have a job after the first period, plus we add the $theta = 10\%$ of the unemployed who find a job. In the second row, we add the 4% of people who lost their job to the 90% of those who were unemployed who remain unemployed.
Now we write a little routine to do this over many periods:
T=50
#create empty array Tx2
xs = np.zeros((T, 2))
#set initial value
xs[0, :] = x0
#iterate over time
for t in range(1, T):
xs[t, :] = np.dot(A.T, xs[t-1, :])
xs
array([[1000000. , 150000. ],
[ 975000. , 175000. ],
[ 953500. , 196500. ],
[ 935010. , 214990. ],
[ 919108.6 , 230891.4 ],
[ 905433.396 , 244566.604 ],
[ 893672.72056 , 256327.27944 ],
[ 883558.5396816 , 266441.4603184 ],
[ 874860.34412618, 275139.65587382],
[ 867379.89594851, 282620.10405149],
[ 860946.71051572, 289053.28948428],
[ 855414.17104352, 294585.82895648],
[ 850656.18709743, 299343.81290257],
[ 846564.32090379, 303435.67909621],
[ 843045.31597726, 306954.68402274],
[ 840018.97174044, 309981.02825956],
[ 837416.31569678, 312583.68430322],
[ 835178.03149923, 314821.96850077],
[ 833253.10708934, 316746.89291066],
[ 831597.67209683, 318402.32790317],
[ 830173.99800327, 319826.00199673],
[ 828949.63828282, 321050.36171718],
[ 827896.68892322, 322103.31107678],
[ 826991.15247397, 323008.84752603],
[ 826212.39112761, 323787.60887239],
[ 825542.65636975, 324457.34363025],
[ 824966.68447798, 325033.31552202],
[ 824471.34865107, 325528.65134893],
[ 824045.35983992, 325954.64016008],
[ 823679.00946233, 326320.99053767],
[ 823363.9481376 , 326636.0518624 ],
[ 823092.99539834, 326907.00460166],
[ 822859.97604257, 327140.02395743],
[ 822659.57939661, 327340.42060339],
[ 822487.23828109, 327512.76171891],
[ 822339.02492173, 327660.97507827],
[ 822211.56143269, 327788.43856731],
[ 822101.94283211, 327898.05716789],
[ 822007.67083562, 327992.32916438],
[ 821926.59691863, 328073.40308137],
[ 821856.87335002, 328143.12664998],
[ 821796.91108102, 328203.08891898],
[ 821745.34352968, 328254.65647032],
[ 821700.99543552, 328299.00456448],
[ 821662.85607455, 328337.14392545],
[ 821630.05622411, 328369.94377589],
[ 821601.84835274, 328398.15164726],
[ 821577.58958335, 328422.41041665],
[ 821556.72704168, 328443.27295832],
[ 821538.78525585, 328461.21474415]])
#import plt
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(xs[:,0], label='Employed')
ax.plot(xs[:,1], label='Unemployed')
ax.legend()
<matplotlib.legend.Legend at 0x7feb41e5d220>
We see that we are converging on constant values for employment and unemployment. As we continue to iterate, we get to a situation where we get:
$$\mathbf{A^T}x = x$$That should look familier. x is then a eigenvector of $\mathbf{A}$
Let us calculate the eigenvectors and eigenvalues of A:
from numpy.linalg import eig
evalues,evectors=eig(A.T)
print('Eigen value:', evalues)
print('Eigen vector', evectors)
Eigen value: [1. 0.86] Eigen vector [[ 0.92847669 -0.70710678] [ 0.37139068 0.70710678]]
# taking the eigen vector in the 0-column
evector = evectors[:,0]
#normalize so that the two values sum to 1
evector = evector/evector.sum()
evector
array([0.71428571, 0.28571429])
# we look at the last iteration of our calculation, which we also normalize
xs[-1,:]/np.sum(xs[-1,:])
array([0.71438155, 0.28561845])
We can confirm, we managed to iterate to eigenvector, representing the steady state of our system.