import numpy as np
import matplotlib.pyplot as plt
import scipy
from numpy.linalg import cond, matrix_rank, norm
from scipy.linalg import inv, solve, det, eig, lu, eigvals
from scipy.linalg import solve_triangular, eigvalsh, cholesky
# The rank in Question 4 depends on an almost-zero singular value, so
# other seeds may report rank 4 instead of 3
np.random.seed(4)ECON526: Problem Set 3
Student Name/Number:
Instructions
- Ensure you modify the field above with your name and student number above immediately
- Modify directly and do not rename the file as Canvas will automatically append your name to it.
- Submit both the
.ipynband an exportedhtmlversion of this file. - See here for instructions.
- Do not change the filenames. Canvas automatically appends your name/student number to submitted files.
Setup
Feel free to use the following packages (and we have added a few convenience imports)
Question 1
Consider the unemployment dynamics from the lecture, where each year a fraction \(\alpha\) of employed workers lose their job and a fraction \(\phi\) of unemployed workers find one, with no births or deaths. Take \(\alpha = 0.03\), \(\phi = 0.15\), and \(X_0 = \begin{bmatrix} E_0 \\ U_0 \end{bmatrix} = \begin{bmatrix} 900000 \\ 100000 \end{bmatrix}\).
\[ \underbrace{\begin{bmatrix} E_{t+1}\\U_{t+1}\end{bmatrix}}_{X_{t+1}} = \underbrace{\begin{bmatrix} 1-\alpha & \phi \\ \alpha & 1-\phi \end{bmatrix}}_{A} \underbrace{\begin{bmatrix} E_{t}\\U_{t}\end{bmatrix}}_{X_t} \]
Question 1.1
Using the simulate function from the lecture, simulate \(X_t\) for \(T = 50\) periods and report \(X_T\).
def simulate(A, X_0, T):
X = np.zeros((len(X_0), T+1))
X[:,0] = X_0
for t in range(T):
X[:,t+1] = A @ X[:,t]
return X
alpha, phi = 0.03, 0.15
X_0 = np.array([900000, 100000])
T = 50
# modify here
# A = ...
# X = simulate(A, X_0, T)Question 1.2
Calculate the eigenvalues and eigenvectors of \(A\) and check that one of the eigenvalues is \(1\). Rescale the corresponding eigenvector so that it sums to the total population \(E_0 + U_0\) to get the steady state \(\bar{X}\), and compare it to your \(X_T\) from Q1.1.
# modify hereQuestion 1.3
Show that the steady-state unemployment rate \(\bar{U} / (\bar{E} + \bar{U})\) is \(\alpha / (\alpha + \phi)\), and verify it numerically using your \(\bar{X}\) from Q1.2.
Answer:
(double click to edit your answer)
# modify here to verify numericallyQuestion 2
Consider five web pages \(A, B, C, D, E\) with the adjacency matrix below, in the convention of the lecture: row \(i\) is the page you are on and column \(j\) is \(1\) if it links to page \(j\), in the order \(A, B, C, D, E\).
M = np.array([[0, 1, 1, 1, 0],
[0, 0, 1, 0, 0],
[1, 0, 0, 0, 0],
[0, 0, 1, 0, 1],
[1, 0, 0, 0, 0]])Question 2.1
Count the incoming links of each page (i.e. the column sums of \(M\)). Then form the matrix \(S\) where each row of \(M\) is divided by the number of links on that page, and verify that every row of \(S\) sums to one.
# modify here
# S = ...Question 2.2
Start on page \(A\) for certain, i.e. \(r_0 = \begin{bmatrix} 1 & 0 & 0 & 0 & 0\end{bmatrix}^{\top}\), and iterate \(r_{t+1} = S^{\top} r_t\) for \(T = 50\) clicks. Report \(r_T\). Hint: simulate from Q1.1 works for any square matrix.
# modify hereQuestion 2.3
Find the eigenvector of \(S^{\top}\) with eigenvalue \(1\), rescale it to sum to one, and compare it to \(r_T\) from Q2.2. Which page has the highest rank? Is it the page with the most incoming links? Explain using the structure of the links.
# modify hereAnswer:
(double click to edit your answer)
Question 3
Generate a design matrix \(X \in {\mathbb{R}}^{200 \times 3}\) with independent standard normal entries and simulate \(y = X \beta + \epsilon\) for the \(\beta\) below with a small amount of noise.
N, K = 200, 3
X = np.random.randn(N, K)
beta = np.array([1.0, -2.0, 0.5])
y = X @ beta + 0.1 * np.random.randn(N)Question 3.1
Estimate \(\beta\) with scipy.linalg.lstsq and compare it to the true \(\beta\).
# modify hereQuestion 3.2
Estimate \(\beta\) again by forming the normal equations \((X^{\top} X) \beta = X^{\top} y\) and solving them with solve(..., assume_a="pos"). Compare to Q3.1.
# modify hereQuestion 3.3
Explain why \(X^{\top} X\) is positive definite for this \(X\), and what that implies about the solution to the least squares problem.
Answer:
(double click to edit your answer)
Question 4
Now append a fourth column to \(X\) which is a linear combination of the others.
X_col = np.column_stack([X, X[:, 0] + 2 * X[:, 1]])
print(X_col.shape)(200, 4)
Question 4.1
Calculate the rank of X_col and the condition number of \(X_{\text{col}}^{\top} X_{\text{col}}\). Compare to the same quantities for \(X\).
# modify hereQuestion 4.2
Try to solve the normal equations for X_col with solve, and then estimate with lstsq. Describe what happens in each case.
# modify hereQuestion 4.3
Explain why the individual parameters are not identified for X_col, and what the lstsq answer represents.
Answer:
(double click to edit your answer)