import numpy as np
import time

N = 5_000_000
xs = [n * 1e-6 for n in range(N)]                 # a Python list of floats
x_n = np.arange(N, dtype=np.float64) * 1e-6       # the same numbers, as an array

t0 = time.perf_counter()
s_py = sum(x * x for x in xs)
t_py = time.perf_counter() - t0

t0 = time.perf_counter()
s_np = float(np.dot(x_n, x_n))                    # fused multiply-and-sum, via BLAS
t_np = time.perf_counter() - t0

print(f"Python loop: {t_py * 1000:8.1f} ms   sum = {s_py:.6f}")
print(f"NumPy (dot): {t_np * 1000:8.1f} ms   sum = {s_np:.6f}")
print(f"speedup:     {t_py / t_np:8.1f}x")
print(f"\nsame answer to {abs(s_py - s_np) / abs(s_py):.1e} relative   (tolerance 1e-12)")
assert abs(s_py - s_np) / abs(s_py) < 1e-12

Python loop:    471.0 ms   sum = 41666654.166667
NumPy (dot):     14.2 ms   sum = 41666654.166668
speedup:         33.1x

same answer to 1.8e-16 relative   (tolerance 1e-12)

# from a Python list, or a list of lists for higher dimension
x_i = np.array([1.0, 2.0, 3.0, 4.0])
M_i_j = np.array([[1, 2, 3],
                  [4, 5, 6]])
x_i, M_i_j

(array([1., 2., 3., 4.]),
 array([[1, 2, 3],
        [4, 5, 6]]))

# fixed-content factories
np.zeros(5), np.ones((2, 3)), np.full((2, 2), 3.14)

(array([0., 0., 0., 0., 0.]),
 array([[1., 1., 1.],
        [1., 1., 1.]]),
 array([[3.14, 3.14],
        [3.14, 3.14]]))

# regular sequences, and the identity
np.arange(0, 10, 2), np.linspace(0, 1, 5), np.eye(3)

(array([0, 2, 4, 6, 8]),
 array([0.  , 0.25, 0.5 , 0.75, 1.  ]),
 array([[1., 0., 0.],
        [0., 1., 0.],
        [0., 0., 1.]]))

# reshape: the same data, a new shape
np.arange(12).reshape(3, 4)

array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

I, J = 3, 4
M_i_j = np.arange(I * J).reshape(I, J)

print("shape:", M_i_j.shape)     # (I, J) -- and the name says so
print("ndim :", M_i_j.ndim)      # number of axes
print("size :", M_i_j.size)      # total number of elements
print("dtype:", M_i_j.dtype)     # element type

shape: (3, 4)
ndim : 2
size : 12
dtype: int32

n_i = np.array([1, 2, 3])              # int
x_i = np.array([0.5, 0.5, 0.5])        # float
print("promotion:", (n_i + x_i).dtype, "-- float wins")

print("explicit cast:", np.arange(6).reshape(2, 3).astype(np.float64).dtype)

# fixed-width integers wrap silently -- unlike Python's built-in int (fd01 section 5)
big = np.array([2 ** 62], dtype=np.int64)
print(f"\n2^62 * 4 as int64: {(big * 4)[0]}   <- silent overflow, no warning")
print(f"2^62 * 4 in pure Python: {2 ** 62 * 4}")

promotion: float64 -- float wins
explicit cast: float64

2^62 * 4 as int64: 0   <- silent overflow, no warning
2^62 * 4 in pure Python: 18446744073709551616

I, J = 4, 5
M_i_j = np.arange(I * J).reshape(I, J)
print(M_i_j)
print("\nM_i_j[1, 2] =", M_i_j[1, 2], "  M_i_j[-1, -1] =", M_i_j[-1, -1])

[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]
 [15 16 17 18 19]]

M_i_j[1, 2] = 7   M_i_j[-1, -1] = 19

print("M_i_j[:2, :3]:\n", M_i_j[:2, :3])
print("\nM_i_j[1:, ::2]  (rows from 1, every second column):\n", M_i_j[1:, ::2])

M_i_j[:2, :3]:
 [[0 1 2]
 [5 6 7]]

M_i_j[1:, ::2]  (rows from 1, every second column):
 [[ 5  7  9]
 [10 12 14]
 [15 17 19]]

M_i_j = np.arange(20).reshape(4, 5)
B_i_j = M_i_j[:2, :3]              # a VIEW of the top-left block
B_i_j[0, 0] = 999

print("mutating the slice changed the original:")
print(M_i_j[:2, :3])
print("share memory?", np.may_share_memory(M_i_j, B_i_j))

M_i_j = np.arange(20).reshape(4, 5)
B_i_j = M_i_j[:2, :3].copy()       # an explicit COPY
B_i_j[0, 0] = 999

print("\nwith .copy(), the original is untouched:", M_i_j[0, 0])
print("share memory?", np.may_share_memory(M_i_j, B_i_j))

mutating the slice changed the original:
[[999   1   2]
 [  5   6   7]]
share memory? True

with .copy(), the original is untouched: 0
share memory? False

x_i = np.array([10, 20, 30, 40, 50])
i_p = np.array([0, 2, 4, 4, 0])        # arbitrary positions, repeats allowed
x_p = x_i[i_p]
print("x_p =", x_p, "  shape", x_p.shape)

x_p = [10 30 50 50 10]   shape (5,)

x_i = np.arange(10)
mask_i = (x_i % 2 == 0) & (x_i > 3)
print("masked:", x_i[mask_i])

# a common pattern: filter, then assign in place
x_i = np.array([-2.0, 1.0, -3.0, 4.0, -1.0])
x_i[x_i < 0] = 0.0                     # the ReLU, in one line
print("after clipping negatives:", x_i)

# what happens if you forget the parentheses
try:
    _ = x_i % 2 == 0 and x_i > 3
except ValueError as err:
    print("\n'and' on arrays ->", type(err).__name__, ":", str(err)[:60], "...")

masked: [4 6 8]
after clipping negatives: [0. 1. 0. 4. 0.]

'and' on arrays -> ValueError : The truth value of an array with more than one element is am ...

x_i = np.linspace(0, 2 * np.pi, 5)
print("sin:", np.sin(x_i).round(6))

x_i = np.array([1.0, 2.0, 3.0])
y_i = np.array([10.0, 20.0, 30.0])
print("\nx+y:", x_i + y_i, "  x*y:", x_i * y_i, "  x**2:", x_i ** 2)

sin: [ 0.  1.  0. -1. -0.]

x+y: [11. 22. 33.]   x*y: [10. 40. 90.]   x**2: [1. 4. 9.]

x_j = np.array([1.0, 2.0, 3.0])        # (J,)  = (3,)
M_i_j = np.zeros((4, 3))               # (I, J) = (4, 3)
print("M_i_j + x_j broadcasts across rows:\n", M_i_j + x_j)

y_i = np.array([10, 20])               # (I,) = (2,)
x_j = np.array([1, 2, 3])              # (J,) = (3,)
print("\ny_i[:, None] + x_j[None, :] -> (I, J):\n", y_i[:, None] + x_j[None, :])

try:
    np.array([1, 2, 3]) + np.array([1, 2, 3, 4])
except ValueError as err:
    print("\nincompatible shapes ->", err)

M_i_j + x_j broadcasts across rows:
 [[1. 2. 3.]
 [1. 2. 3.]
 [1. 2. 3.]
 [1. 2. 3.]]

y_i[:, None] + x_j[None, :] -> (I, J):
 [[11 12 13]
 [21 22 23]]

incompatible shapes -> operands could not be broadcast together with shapes (3,) (4,)

rng = np.random.default_rng(5)         # seeded once for the lecture; the API is section 14

I, J, K = 5, 7, 2
P_i_k = rng.normal(size=(I, K))
Q_j_k = rng.normal(size=(J, K))

diff_i_j_k = P_i_k[:, None, :] - Q_j_k[None, :, :]     # (I, 1, K) - (1, J, K) -> (I, J, K)
D_i_j = np.sqrt((diff_i_j_k ** 2).sum(axis=-1))        # collapse k -> (I, J)

print("shapes:", P_i_k.shape, Q_j_k.shape, "->", diff_i_j_k.shape, "->", D_i_j.shape)

# verification: an explicit double loop must give the same matrix
D_check_i_j = np.empty((I, J))
for i in range(I):
    for j in range(J):
        D_check_i_j[i, j] = np.sqrt(((P_i_k[i] - Q_j_k[j]) ** 2).sum())

gap = np.abs(D_i_j - D_check_i_j).max()
print(f"max |broadcast - loop| = {gap:.2e}   (tolerance 1e-12)")
assert gap < 1e-12
print("check passed.")

shapes: (5, 2) (7, 2) -> (5, 7, 2) -> (5, 7)
max |broadcast - loop| = 0.00e+00   (tolerance 1e-12)
check passed.

I, J = 3, 4
M_i_j = np.arange(I * J).reshape(I, J)
print(M_i_j)
print("\nsum()       =", M_i_j.sum(), "        (everything)")
print("sum(axis=0) =", M_i_j.sum(axis=0), "  shape", M_i_j.sum(axis=0).shape, "-> an s_j")
print("sum(axis=1) =", M_i_j.sum(axis=1), "      shape", M_i_j.sum(axis=1).shape, "-> an s_i")

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

sum()       = 66         (everything)
sum(axis=0) = [12 15 18 21]   shape (4,) -> an s_j
sum(axis=1) = [ 6 22 38]       shape (3,) -> an s_i

mean_i_1 = M_i_j.mean(axis=1, keepdims=True)      # (I, 1)
demeaned_i_j = M_i_j - mean_i_1                   # broadcasts back to (I, J)

print("mean_i_1 shape:", mean_i_1.shape)
print("row means after demeaning:", demeaned_i_j.mean(axis=1).round(12))
assert np.allclose(demeaned_i_j.mean(axis=1), 0.0, atol=1e-12)
print("check passed: every row now has mean zero (tolerance 1e-12).")

mean_i_1 shape: (3, 1)
row means after demeaning: [0. 0. 0.]
check passed: every row now has mean zero (tolerance 1e-12).

A_i_j = np.array([[3.0, 2.0], [1.0, 2.0]])
b_i = np.array([7.0, 5.0])

x_j = np.linalg.solve(A_i_j, b_i)
residual = np.abs(A_i_j @ x_j - b_i).max()

print("x =", x_j)
print(f"max |A x - b| = {residual:.2e}   (tolerance 1e-12)")
assert residual < 1e-12
print("check passed.")

x = [1. 2.]
max |A x - b| = 0.00e+00   (tolerance 1e-12)
check passed.

N, T, K = 4, 5, 3
face_n = np.array([1000.0, 1000.0, 5000.0, 1000.0])
coupon_n = np.array([50.0, 80.0, 250.0, 100.0])
r_k = np.array([0.02, 0.05, 0.08])

# cash flows, shape (N, T): a coupon each year, plus the face value at maturity
C_n_t = np.tile(coupon_n[:, None], (1, T))
C_n_t[:, -1] += face_n

# the price table, shape (K, T): D_k_t[k, t] is the price at rate r_k of a dollar at t+1
period_t = np.arange(1, T + 1)
D_k_t = (1 + r_k[:, None]) ** (-period_t[None, :])

PV_k_n = D_k_t @ C_n_t.T                                # (K, T) @ (T, N) -> (K, N)

print("price table D_k_t (rows = rates, cols = dates):\n", D_k_t.round(4))
print("\ncash flows C_n_t (rows = bonds, cols = dates):\n", C_n_t)
print("\nPV_k_n (rows = rates, cols = bonds):\n", PV_k_n.round(2))

price table D_k_t (rows = rates, cols = dates):
 [[0.9804 0.9612 0.9423 0.9238 0.9057]
 [0.9524 0.907  0.8638 0.8227 0.7835]
 [0.9259 0.8573 0.7938 0.735  0.6806]]

cash flows C_n_t (rows = bonds, cols = dates):
 [[  50.   50.   50.   50. 1050.]
 [  80.   80.   80.   80. 1080.]
 [ 250.  250.  250.  250. 5250.]
 [ 100.  100.  100.  100. 1100.]]

PV_k_n (rows = rates, cols = bonds):
 [[1141.4  1282.81 5707.02 1377.08]
 [1000.   1129.88 5000.   1216.47]
 [ 880.22 1000.   4401.09 1079.85]]

# verification against the fd01 loop, bond by bond and rate by rate
PV_check_k_n = np.empty((K, N))
for k in range(K):
    for n in range(N):
        PV_check_k_n[k, n] = sum(C_n_t[n, t] / (1 + r_k[k]) ** (t + 1) for t in range(T))

gap = np.abs(PV_k_n - PV_check_k_n).max()
print(f"max |vectorized - loop| = {gap:.2e}   (tolerance 1e-9)")
assert gap < 1e-9

# and the economics: value falls as the discount rate rises, for every bond
assert np.all(np.diff(PV_k_n, axis=0) < 0)
print("check passed: 12 present values, one matrix product, and PV is decreasing in r.")

max |vectorized - loop| = 9.09e-13   (tolerance 1e-9)
check passed: 12 present values, one matrix product, and PV is decreasing in r.

rng = np.random.default_rng(5)

N, K = 200, 3
X_n_k = rng.normal(size=(N, K))
beta_k = np.array([1.0, -0.5, 2.0])
y_n = X_n_k @ beta_k + rng.normal(scale=0.5, size=N)

betahat_k = np.linalg.solve(X_n_k.T @ X_n_k, X_n_k.T @ y_n)      # never inv()

# independent route: the least-squares driver, a different algorithm (QR, not normal equations)
betahat_lstsq_k = np.linalg.lstsq(X_n_k, y_n, rcond=None)[0]

print("true beta     :", beta_k)
print("solve()       :", betahat_k.round(6))
print("lstsq()       :", betahat_lstsq_k.round(6))

gap = np.abs(betahat_k - betahat_lstsq_k).max()
print(f"\nmax |solve - lstsq| = {gap:.2e}   (tolerance 1e-9)")
assert gap < 1e-9

# the residuals must be orthogonal to the regressors -- the defining property of OLS
resid_n = y_n - X_n_k @ betahat_k
orth = np.abs(X_n_k.T @ resid_n).max()
print(f"max |X' e|          = {orth:.2e}   (tolerance 1e-9)")
assert orth < 1e-9
print("check passed: two algorithms agree, and the normal equations hold.")

true beta     : [ 1.  -0.5  2. ]
solve()       : [ 0.946551 -0.5541    2.004688]
lstsq()       : [ 0.946551 -0.5541    2.004688]

max |solve - lstsq| = 1.78e-15   (tolerance 1e-9)
max |X' e|          = 3.11e-13   (tolerance 1e-9)
check passed: two algorithms agree, and the normal equations hold.

A_i_j = np.array([[1, 2],
                  [3, 4]])
Id_k_l = np.eye(2, dtype=int)

print("A (x) I  -- copies of I scaled by the entries of A:\n", np.kron(A_i_j, Id_k_l))
print("\nI (x) A  -- copies of A down the diagonal:\n", np.kron(Id_k_l, A_i_j))

A (x) I  -- copies of I scaled by the entries of A:
 [[1 0 2 0]
 [0 1 0 2]
 [3 0 4 0]
 [0 3 0 4]]

I (x) A  -- copies of A down the diagonal:
 [[1 2 0 0]
 [3 4 0 0]
 [0 0 1 2]
 [0 0 3 4]]

rng = np.random.default_rng(11)
I, J, K, L, M, Nn = 3, 4, 5, 2, 2, 3
A_i_j = rng.normal(size=(I, J))
B_k_l = rng.normal(size=(K, L))
C_j_m = rng.normal(size=(J, M))
D_l_n = rng.normal(size=(L, Nn))

# mixed product: (A (x) B)(C (x) D) = (AC) (x) (BD)
lhs = np.kron(A_i_j, B_k_l) @ np.kron(C_j_m, D_l_n)
rhs = np.kron(A_i_j @ C_j_m, B_k_l @ D_l_n)
gap_mixed = np.abs(lhs - rhs).max()
print(f"mixed product      : max gap {gap_mixed:.2e}")
assert gap_mixed < 1e-10

# transpose: (A (x) B)' = A' (x) B'
gap_tr = np.abs(np.kron(A_i_j, B_k_l).T - np.kron(A_i_j.T, B_k_l.T)).max()
print(f"transpose identity : max gap {gap_tr:.2e}")
assert gap_tr < 1e-10
print("both hold (tolerance 1e-10).")

mixed product      : max gap 1.33e-15
transpose identity : max gap 0.00e+00
both hold (tolerance 1e-10).

X_j_k = rng.normal(size=(J, K))

vec_col = lambda M_a_b: M_a_b.T.reshape(-1)      # column stacking: vec
vec_C = lambda M_a_b: M_a_b.reshape(-1)          # row-major: vec_C, the house convention

# (12.1) column-stacking:  vec(A X B) = (B' (x) A) vec(X)
gap_col = np.abs(vec_col(A_i_j @ X_j_k @ B_k_l)
                 - np.kron(B_k_l.T, A_i_j) @ vec_col(X_j_k)).max()

# (12.2) row-major:        vec_C(A X B) = (A (x) B') vec_C(X)
gap_row = np.abs(vec_C(A_i_j @ X_j_k @ B_k_l)
                 - np.kron(A_i_j, B_k_l.T) @ vec_C(X_j_k)).max()

print(f"(12.1) column-stacking vec  : max gap {gap_col:.2e}")
print(f"(12.2) row-major   vec_C    : max gap {gap_row:.2e}")
assert gap_col < 1e-10 and gap_row < 1e-10

# Now the dangerous part: pair the column-stacking Kronecker matrix with a
# row-major vector. Does it complain?
wrong = np.kron(B_k_l.T, A_i_j) @ vec_C(X_j_k)
right = np.kron(A_i_j, B_k_l.T) @ vec_C(X_j_k)

print(f"\nwrong pairing: shape {wrong.shape}  -- no exception raised")
print(f"right pairing: shape {right.shape}")
print(f"same shape: {wrong.shape == right.shape}    same values: {np.allclose(wrong, right)}")
print(f"max |wrong - right| = {np.abs(wrong - right).max():.3f}")

(12.1) column-stacking vec  : max gap 8.88e-16
(12.2) row-major   vec_C    : max gap 8.88e-16

wrong pairing: shape (6,)  -- no exception raised
right pairing: shape (6,)
same shape: True    same values: False
max |wrong - right| = 9.327

n = 60
A_i_j = rng.normal(size=(n, n))
B_k_l = rng.normal(size=(n, n))
X_j_k = rng.normal(size=(n, n))

t0 = time.perf_counter()
big = np.kron(A_i_j, B_k_l.T)               # (n^2, n^2) = (3600, 3600)
y_kron = big @ vec_C(X_j_k)
t_kron = time.perf_counter() - t0

t0 = time.perf_counter()
y_direct = vec_C(A_i_j @ X_j_k @ B_k_l)
t_direct = time.perf_counter() - t0

gap = np.abs(y_kron - y_direct).max()
print(f"same answer: max gap {gap:.2e}   (tolerance 1e-8)")
assert gap < 1e-8
print(f"materialized Kronecker: {t_kron * 1000:8.2f} ms   ({big.nbytes / 1e6:.0f} MB)")
print(f"vec trick (A @ X @ B) : {t_direct * 1000:8.2f} ms")
print(f"speedup: {t_kron / t_direct:.0f}x, at {big.nbytes / 1e6:.0f} MB less memory")

same answer: max gap 4.26e-13   (tolerance 1e-8)
materialized Kronecker:    77.24 ms   (104 MB)
vec trick (A @ X @ B) :     1.22 ms
speedup: 63x, at 104 MB less memory

X = ["analyst", "engineer", "manager"]
Y = ["bank", "startup", "university", "government"]
nbx, nby = len(X), len(Y)

Phi_x_y = np.array([[10.0, 9.0, 5.0, 4.0],       # the surplus table of fd02 section 12
                    [ 9.0, 6.0, 4.0, 3.0],
                    [ 3.0, 4.0, 8.0, 7.0]])

Phi_xy = Phi_x_y.reshape(-1)                     # vec_C

# the fd02 index map, rebuilt by hand, must agree with reshape(-1)
idx_x_y = {(i, j): i * nby + j for i in range(nbx) for j in range(nby)}
gap = max(abs(Phi_xy[idx_x_y[(i, j)]] - Phi_x_y[i, j])
          for i in range(nbx) for j in range(nby))

print(f"Phi_x_y.shape = {Phi_x_y.shape},  Phi_xy.shape = {Phi_xy.shape}")
print(f"Phi_xy = {Phi_xy}")
print(f"\nmax |reshape(-1) - hand-built index map| = {gap:.1e}   (exact; tolerance 0)")
assert gap == 0.0
print("check passed: reshape(-1) IS the fd02 index map.")

Phi_x_y.shape = (3, 4),  Phi_xy.shape = (12,)
Phi_xy = [10.  9.  5.  4.  9.  6.  4.  3.  3.  4.  8.  7.]

max |reshape(-1) - hand-built index map| = 0.0e+00   (exact; tolerance 0)
check passed: reshape(-1) IS the fd02 index map.

ones_x = np.ones((1, nbx))
ones_y = np.ones((1, nby))

A_x_xy = np.kron(np.eye(nbx), ones_y)            # (nbx, nbx*nby): sums over y
A_y_xy = np.kron(ones_x, np.eye(nby))            # (nby, nbx*nby): sums over x
A_z_xy = np.vstack([A_x_xy, A_y_xy])             # (nbx+nby, nbx*nby)

print(f"A_x_xy {A_x_xy.shape}:\n{A_x_xy.astype(int)}")
print(f"\nA_y_xy {A_y_xy.shape}:\n{A_y_xy.astype(int)}")
print(f"\nstacked A_z_xy: {A_z_xy.shape}, {int(A_z_xy.sum())} nonzeros out of {A_z_xy.size}")

A_x_xy (3, 12):
[[1 1 1 1 0 0 0 0 0 0 0 0]
 [0 0 0 0 1 1 1 1 0 0 0 0]
 [0 0 0 0 0 0 0 0 1 1 1 1]]

A_y_xy (4, 12):
[[1 0 0 0 1 0 0 0 1 0 0 0]
 [0 1 0 0 0 1 0 0 0 1 0 0]
 [0 0 1 0 0 0 1 0 0 0 1 0]
 [0 0 0 1 0 0 0 1 0 0 0 1]]

stacked A_z_xy: (7, 12), 24 nonzeros out of 84

# verification 1: the operator must reproduce the margins of an arbitrary mu
mu_x_y = rng.uniform(size=(nbx, nby))
margins = A_z_xy @ mu_x_y.reshape(-1)

gap_rows = np.abs(margins[:nbx] - mu_x_y.sum(axis=1)).max()
gap_cols = np.abs(margins[nbx:] - mu_x_y.sum(axis=0)).max()
print(f"row margins : max gap {gap_rows:.2e}")
print(f"col margins : max gap {gap_cols:.2e}   (tolerance 1e-12)")
assert max(gap_rows, gap_cols) < 1e-12
print("check passed: A_z_xy reads off both families of margins.")

row margins : max gap 0.00e+00
col margins : max gap 0.00e+00   (tolerance 1e-12)
check passed: A_z_xy reads off both families of margins.

u_x = rng.normal(size=nbx)
v_y = rng.normal(size=nby)

lhs_xy = A_z_xy.T @ np.concatenate([u_x, v_y])       # (nbx*nby,)
rhs_x_y = u_x[:, None] + v_y[None, :]                # (nbx, nby) by broadcasting

gap = np.abs(lhs_xy - rhs_x_y.reshape(-1)).max()
print(f"max |A' [u; v]  -  (u_x + v_y)| = {gap:.2e}   (tolerance 1e-12)")
assert gap < 1e-12
print("check passed: the transpose of the constraint matrix generates the dual constraint.")

max |A' [u; v]  -  (u_x + v_y)| = 0.00e+00   (tolerance 1e-12)
check passed: the transpose of the constraint matrix generates the dual constraint.

# the optimal matching found by enumeration in fd02: analyst->startup,
# engineer->bank, manager->university. Its objective must be the 26.0 found there.
mu_star_x_y = np.zeros((nbx, nby))
mu_star_x_y[0, 1] = mu_star_x_y[1, 0] = mu_star_x_y[2, 2] = 1.0

objective = Phi_xy @ mu_star_x_y.reshape(-1)
margins = A_z_xy @ mu_star_x_y.reshape(-1)

print(f"objective <Phi, mu*> = {objective}   (fd02 found 26.0)")
assert objective == 26.0
print(f"worker margins n = {margins[:nbx].astype(int)}   (one job each)")
print(f"firm margins   m = {margins[nbx:].astype(int)}   (government unfilled)")
print("\ncheck passed: the matrix formulation reproduces the fd02 result exactly.")

objective <Phi, mu*> = 26.0   (fd02 found 26.0)
worker margins n = [1 1 1]   (one job each)
firm margins   m = [1 1 1 0]   (government unfilled)

check passed: the matrix formulation reproduces the fd02 result exactly.

import math

rng = np.random.default_rng(seed=42)

print("normal(0,1), 5 draws :", rng.normal(loc=0.0, scale=1.0, size=5).round(4))
print("uniform on [0,1)     :", rng.random(size=5).round(4))
print("integers in {0,...,9}:", rng.integers(low=0, high=10, size=10))

normal(0,1), 5 draws : [ 0.3047 -1.04    0.7505  0.9406 -1.951 ]
uniform on [0,1)     : [0.9756 0.7611 0.7861 0.1281 0.4504]
integers in {0,...,9}: [5 3 1 9 7 6 4 8 5 4]

# the same seed gives the same stream -- this is the whole point
a_i = np.random.default_rng(seed=2025).normal(size=4)
b_i = np.random.default_rng(seed=2025).normal(size=4)
c_i = np.random.default_rng(seed=2026).normal(size=4)

print("seed 2025:", a_i.round(6))
print("seed 2025:", b_i.round(6), "  identical:", np.array_equal(a_i, b_i))
print("seed 2026:", c_i.round(6), "  identical:", np.array_equal(a_i, c_i))
assert np.array_equal(a_i, b_i) and not np.array_equal(a_i, c_i)
print("\ncheck passed: the stream is a deterministic function of the seed.")

seed 2025: [-2.221254  0.026    -0.538969 -1.129193]
seed 2025: [-2.221254  0.026    -0.538969 -1.129193]   identical: True
seed 2026: [-0.793122  0.240571 -1.896326  1.395772]   identical: False

check passed: the stream is a deterministic function of the seed.

rng = np.random.default_rng(0)

# a correlated bivariate normal: K = 2 coordinates, N draws
Sigma_k_k = np.array([[1.0, 0.8],
                      [0.8, 2.0]])
N = 100_000
X_n_k = rng.multivariate_normal(mean=np.zeros(2), cov=Sigma_k_k, size=N)

print("shape:", X_n_k.shape, " -- N draws of a K-vector")
print("\nsample mean       :", X_n_k.mean(axis=0).round(4), " (true: [0, 0])")
print("sample covariance :\n", np.cov(X_n_k.T).round(4))
print("\ntrue covariance   :\n", Sigma_k_k)

gap = np.abs(np.cov(X_n_k.T) - Sigma_k_k).max()
tol = 4 / math.sqrt(N) * Sigma_k_k.max()          # a few Monte Carlo standard errors
print(f"\nmax |sample - true| = {gap:.4f}   (sampling tolerance {tol:.4f})")
assert gap < tol
print("check passed, at a tolerance set by the sampling error -- not by machine precision.")

shape: (100000, 2)  -- N draws of a K-vector

sample mean       : [-0.0004 -0.003 ]  (true: [0, 0])
sample covariance :
 [[0.9969 0.8008]
 [0.8008 2.014 ]]

true covariance   :
 [[1.  0.8]
 [0.8 2. ]]

max |sample - true| = 0.0140   (sampling tolerance 0.0253)
check passed, at a tolerance set by the sampling error -- not by machine precision.

def simulate_returns(rng, T, mu=0.0005, sigma=0.01):
    """Daily log-returns under a Gaussian model. The generator is an argument."""
    return rng.normal(loc=mu, scale=sigma, size=T)

rng = np.random.default_rng(42)
r_t = simulate_returns(rng, T=252)                # one year of trading days

print(f"mean daily return {r_t.mean():+.5f}, daily volatility {r_t.std():.5f}")
print(f"annualized: {252 * r_t.mean():+.3%} drift, {math.sqrt(252) * r_t.std():.1%} volatility")

# calling it twice with the SAME generator gives different draws -- the state advances
rng = np.random.default_rng(42)
first, second = simulate_returns(rng, 3), simulate_returns(rng, 3)
print(f"\nfirst call : {first.round(5)}")
print(f"second call: {second.round(5)}   (different: the generator advanced)")

mean daily return +0.00002, daily volatility 0.00936
annualized: +0.411% drift, 14.9% volatility

first call : [ 0.00355 -0.0099   0.008  ]
second call: [ 0.00991 -0.01901 -0.01252]   (different: the generator advanced)

rng = np.random.default_rng(2025)
N = 2_000_000

X_n_k = rng.random(size=(N, 2))                   # N draws of a K = 2 vector
inside_n = (X_n_k ** 2).sum(axis=1) < 1.0         # a boolean array: the indicator f(X)

pi_hat = 4 * inside_n.mean()
se = 4 * inside_n.std(ddof=1) / math.sqrt(N)      # standard error of the ESTIMATE

print(f"pi_hat  = {pi_hat:.6f} +/- {se:.6f}   (true {math.pi:.6f})")
print(f"error   = {abs(pi_hat - math.pi):.6f}")
print(f"error / se = {abs(pi_hat - math.pi) / se:.2f}   (should be O(1); |z| > 3 would be alarming)")
assert abs(pi_hat - math.pi) < 4 * se
print("\ncheck passed: the error is within four standard errors, as it should be.")

pi_hat  = 3.140212 +/- 0.001162   (true 3.141593)
error   = 0.001381
error / se = 1.19   (should be O(1); |z| > 3 would be alarming)

check passed: the error is within four standard errors, as it should be.

import matplotlib.pyplot as plt          # plotting proper is fd06

rng = np.random.default_rng(0)
N = 5_000
n_n = np.arange(1, N + 1)

xnorm_n = rng.normal(loc=0.0, scale=1.0, size=N)
xexp_n = rng.exponential(scale=1.0, size=N)       # mean 1, sd 1

runnorm_n = np.cumsum(xnorm_n) / n_n
runexp_n = np.cumsum(xexp_n) / n_n

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(runnorm_n, lw=1, label=r"$N(0,1)$, mean 0")
ax.plot(runexp_n, lw=1, label=r"$\mathrm{Exp}(1)$, mean 1")
ax.plot(n_n, 2 / np.sqrt(n_n), "k:", lw=1, label=r"$\pm 2/\sqrt{n}$ band")
ax.plot(n_n, -2 / np.sqrt(n_n), "k:", lw=1)
ax.plot(n_n, 1 + 2 / np.sqrt(n_n), "k:", lw=1)
ax.plot(n_n, 1 - 2 / np.sqrt(n_n), "k:", lw=1)
ax.set_xlabel("n"); ax.set_ylabel("running mean")
ax.set_title("The law of large numbers, with the theoretical band")
ax.legend(); plt.tight_layout(); plt.show()

# A single path wanders in and out of the band, so counting how long it stays inside
# tests nothing. The rate is a statement about the DISTRIBUTION of the sample mean,
# so verify it across many independent samples: sd(mean of n draws) should be sigma/sqrt(n).
rng = np.random.default_rng(3)
R_paths = 2_000

print(f"{'n':>7}{'sd of sample mean':>20}{'sigma / sqrt(n)':>18}{'ratio':>9}")
print("-" * 54)
worst = 0.0
for n_draws in (100, 1_000, 10_000):
    means_r = rng.exponential(scale=1.0, size=(R_paths, n_draws)).mean(axis=1)
    observed = means_r.std(ddof=1)
    predicted = 1.0 / math.sqrt(n_draws)          # Exp(1) has sigma = 1
    worst = max(worst, abs(observed / predicted - 1))
    print(f"{n_draws:>7}{observed:>20.6f}{predicted:>18.6f}{observed / predicted:>9.4f}")

print("-" * 54)
print(f"worst relative departure = {worst:.1%}   (tolerance 5%, itself a sampling limit)")
assert worst < 0.05
print("check passed: the sample mean's dispersion falls exactly like 1/sqrt(n).")

      n   sd of sample mean   sigma / sqrt(n)    ratio
------------------------------------------------------
    100            0.100978          0.100000   1.0098
   1000            0.032282          0.031623   1.0209

  10000            0.009867          0.010000   0.9867
------------------------------------------------------
worst relative departure = 2.1%   (tolerance 5%, itself a sampling limit)
check passed: the sample mean's dispersion falls exactly like 1/sqrt(n).

rng = np.random.default_rng(7)

ns = [1, 5, 30, 100]
R = 20_000                                        # replications per panel
mu_exp, sigma_exp = 1.0, 1.0

x_p = np.linspace(-4, 4, 200)
phi_p = np.exp(-x_p ** 2 / 2) / math.sqrt(2 * math.pi)

fig, axes = plt.subplots(1, 4, figsize=(14, 3.2), sharey=True)
skews = []
for ax, n in zip(axes, ns):
    X_r_n = rng.exponential(scale=1.0, size=(R, n))
    Z_r = (X_r_n.mean(axis=1) - mu_exp) / (sigma_exp / math.sqrt(n))
    skews.append(((Z_r - Z_r.mean()) ** 3).mean() / Z_r.std() ** 3)
    ax.hist(Z_r, bins=60, density=True, alpha=0.7)
    ax.plot(x_p, phi_p, "k--", lw=1.2, label=r"$N(0,1)$")
    ax.set_xlim(-4, 4); ax.set_title(f"n = {n}"); ax.set_xlabel(r"$Z_n$")
axes[0].set_ylabel("density"); axes[-1].legend()
plt.tight_layout(); plt.show()

# the skewness of Z_n must decay like 2/sqrt(n) for Exp(1), whose skewness is 2
print(f"{'n':>5}{'observed skew':>16}{'predicted 2/sqrt(n)':>22}{'ratio':>9}")
print("-" * 52)
for n, sk in zip(ns, skews):
    predicted = 2 / math.sqrt(n)
    print(f"{n:>5}{sk:>16.4f}{predicted:>22.4f}{sk / predicted:>9.3f}")

ratios = [sk / (2 / math.sqrt(n)) for n, sk in zip(ns, skews)]
print("-" * 52)
print(f"ratios all within 10% of 1: {all(abs(r - 1) < 0.10 for r in ratios)}")
assert all(abs(r - 1) < 0.10 for r in ratios)
print("check passed: the departure from normality vanishes at the predicted rate.")

    n   observed skew   predicted 2/sqrt(n)    ratio
----------------------------------------------------
    1          1.9675                2.0000    0.984
    5          0.8722                0.8944    0.975
   30          0.3814                0.3651    1.044
  100          0.2037                0.2000    1.019
----------------------------------------------------
ratios all within 10% of 1: True
check passed: the departure from normality vanishes at the predicted rate.

rng = np.random.default_rng(123)

# estimate E[X^2] for X ~ N(0,1); the true value is 1
N = 100_000
fX_n = rng.normal(size=N) ** 2

mu_hat = fX_n.mean()
se = fX_n.std(ddof=1) / math.sqrt(N)
lo, hi = mu_hat - 1.96 * se, mu_hat + 1.96 * se

print(f"estimate   = {mu_hat:.5f}")
print(f"std error  = {se:.5f}")
print(f"95% CI     = ({lo:.5f}, {hi:.5f})")
print(f"true value = 1.00000   -- covered: {lo < 1.0 < hi}")

estimate   = 0.99962
std error  = 0.00448
95% CI     = (0.99084, 1.00840)
true value = 1.00000   -- covered: True

# does a 95% interval actually cover 95% of the time? build 2000 of them and count.
rng = np.random.default_rng(9)
R, n_small = 2_000, 500

X_r_n = rng.normal(size=(R, n_small)) ** 2
mu_r = X_r_n.mean(axis=1)
se_r = X_r_n.std(axis=1, ddof=1) / math.sqrt(n_small)
covered_r = (mu_r - 1.96 * se_r < 1.0) & (1.0 < mu_r + 1.96 * se_r)

rate = covered_r.mean()
se_rate = math.sqrt(rate * (1 - rate) / R)
print(f"coverage of the nominal 95% interval: {rate:.1%} +/- {1.96 * se_rate:.1%} over {R} replications")
assert abs(rate - 0.95) < 4 * se_rate + 0.01
print("check passed: the interval covers at approximately its nominal rate.")

coverage of the nominal 95% interval: 94.2% +/- 1.0% over 2000 replications
check passed: the interval covers at approximately its nominal rate.

# estimate E[max(X, 0)] for X ~ N(0,1); the true value is 1/sqrt(2 pi)
true_value = 1 / math.sqrt(2 * math.pi)
rng = np.random.default_rng(0)
N = 200_000

X_n = rng.normal(size=N)                                  # standard Monte Carlo
f_n = np.maximum(X_n, 0.0)
mu_std, se_std = f_n.mean(), f_n.std(ddof=1) / math.sqrt(N)

M = N // 2                                                # antithetic: M pairs, same budget
X_m = rng.normal(size=M)
f_m = 0.5 * (np.maximum(X_m, 0.0) + np.maximum(-X_m, 0.0))
mu_anti, se_anti = f_m.mean(), f_m.std(ddof=1) / math.sqrt(M)

print(f"true value                    : {true_value:.6f}")
print(f"standard   (N = {N:,})   : {mu_std:.6f} +/- {se_std:.6f}")
print(f"antithetic (M = {M:,} pairs): {mu_anti:.6f} +/- {se_anti:.6f}")
print(f"\nvariance ratio at equal budget: {(se_std / se_anti) ** 2:.2f}x")
assert abs(mu_anti - true_value) < 4 * se_anti
print("check passed: the antithetic estimator is unbiased and strictly more precise.")

true value                    : 0.398942
standard   (N = 200,000)   : 0.399525 +/- 0.001306
antithetic (M = 100,000 pairs): 0.399208 +/- 0.000954

variance ratio at equal budget: 1.87x
check passed: the antithetic estimator is unbiased and strictly more precise.

mu_c, sigma_c, a = 1.0, 0.30, 2.0        # mean consumption, its sd, risk aversion

Eu_exact = -math.exp(-a * mu_c + 0.5 * a ** 2 * sigma_c ** 2)
CE_exact = mu_c - 0.5 * a * sigma_c ** 2
premium_exact = 0.5 * a * sigma_c ** 2

rng = np.random.default_rng(42)
N = 400_000
C_n = mu_c + sigma_c * rng.standard_normal(N)
u_n = -np.exp(-a * C_n)

Eu_mc = u_n.mean()
se_mc = u_n.std(ddof=1) / math.sqrt(N)
CE_mc = -math.log(-Eu_mc) / a
premium_mc = mu_c - CE_mc

print(f"{'':<22}{'closed form':>14}{'Monte Carlo':>16}{'MC std error':>15}")
print("-" * 68)
print(f"{'E[u(C)]':<22}{Eu_exact:>14.6f}{Eu_mc:>16.6f}{se_mc:>15.6f}")
print(f"{'certainty equivalent':<22}{CE_exact:>14.6f}{CE_mc:>16.6f}")
print(f"{'risk premium':<22}{premium_exact:>14.6f}{premium_mc:>16.6f}")
print("-" * 68)
z = abs(Eu_mc - Eu_exact) / se_mc
print(f"|MC - exact| / se = {z:.2f} standard errors")
assert z < 4
print("check passed: the simulation reproduces the closed form to within sampling error.")

                         closed form     Monte Carlo   MC std error
--------------------------------------------------------------------
E[u(C)]                    -0.162026       -0.162033       0.000169
certainty equivalent        0.910000        0.909976
risk premium                0.090000        0.090024
--------------------------------------------------------------------
|MC - exact| / se = 0.05 standard errors
check passed: the simulation reproduces the closed form to within sampling error.

# antithetic variates: u is monotone in C, hence in the underlying normal draw
rng = np.random.default_rng(42)
M = N // 2
Z_m = rng.standard_normal(M)
upair_m = 0.5 * (-np.exp(-a * (mu_c + sigma_c * Z_m)) - np.exp(-a * (mu_c - sigma_c * Z_m)))

Eu_anti = upair_m.mean()
se_anti = upair_m.std(ddof=1) / math.sqrt(M)

print(f"standard   : E[u] = {Eu_mc:.6f} +/- {se_mc:.6f}")
print(f"antithetic : E[u] = {Eu_anti:.6f} +/- {se_anti:.6f}")
print(f"variance ratio at equal budget: {(se_mc / se_anti) ** 2:.2f}x")
assert abs(Eu_anti - Eu_exact) < 4 * se_anti

standard   : E[u] = -0.162033 +/- 0.000169
antithetic : E[u] = -0.162205 +/- 0.000094
variance ratio at equal budget: 3.19x

# your answer here

# your answer here

# your answer here

# your answer here

# your answer here

rng = np.random.default_rng(5)

X_i_j = rng.normal(size=(5, 5))
M_i_j = X_i_j @ X_i_j.T                       # symmetric positive semi-definite

eigvals, eigvecs = np.linalg.eigh(M_i_j)      # ascending
lambda_1, lambda_2 = eigvals[-1], eigvals[-2]
w_1 = eigvecs[:, -1]
predicted_rate = abs(lambda_2 / lambda_1)

v_i = rng.normal(size=5)
v_i /= np.linalg.norm(v_i)

errors = []
for _ in range(80):
    errors.append(np.linalg.norm(v_i - np.dot(v_i, w_1) * w_1))
    Mv_i = M_i_j @ v_i
    v_i = Mv_i / np.linalg.norm(Mv_i)

lambda_hat = v_i @ M_i_j @ v_i
gap = abs(lambda_hat - lambda_1)

print(f"power iteration lambda = {lambda_hat:.12f}")
print(f"np.linalg.eigh  lambda = {lambda_1:.12f}")
print(f"gap = {gap:.2e}   (tolerance 1e-10)")
assert gap < 1e-10
print(f"eigenvector gap = {np.linalg.norm(np.abs(v_i) - np.abs(w_1)):.2e}")

power iteration lambda = 10.271536301373
np.linalg.eigh  lambda = 10.271536301373
gap = 3.55e-15   (tolerance 1e-10)
eigenvector gap = 7.54e-16

# Verify the predicted contraction rate. The rate is asymptotic, so it is measured
# only in the clean band: past the transient (error < 1e-2) and above the rounding
# floor (error > 1e-10). Rows outside the band are printed but not asserted on.
print(f"predicted |lambda_2 / lambda_1| = {predicted_rate:.9f}\n")
print(f"{'k':>4}{'error':>12}{'ratio':>14}{'gap vs theory':>16}   in band")
print("-" * 58)
worst = 0.0
for k in range(len(errors) - 1):
    if errors[k + 1] < 1e-10:
        break
    ratio = errors[k + 1] / errors[k]
    gap = abs(ratio - predicted_rate)
    in_band = errors[k] < 1e-2
    if in_band:
        worst = max(worst, gap)
    if k < 4 or k % 4 == 0:
        print(f"{k:>4}{errors[k]:>12.2e}{ratio:>14.9f}{gap:>16.1e}   {in_band}")

print("-" * 58)
print(f"worst gap inside the band = {worst:.1e}   (tolerance 1e-4)")
assert worst < 1e-4
print("check passed: the error contracts at |lambda_2 / lambda_1|.")

predicted |lambda_2 / lambda_1| = 0.574154668

   k       error         ratio   gap vs theory   in band
----------------------------------------------------------
   0    9.57e-01   0.584325926         1.0e-02   False
   1    5.59e-01   0.606706966         3.3e-02   False
   2    3.39e-01   0.596604927         2.2e-02   False
   3    2.03e-01   0.582155206         8.0e-03   False
   4    1.18e-01   0.576840561         2.7e-03   False
   8    1.29e-02   0.574186686         3.2e-05   False
  12    1.40e-03   0.574155046         3.8e-07   True
  16    1.52e-04   0.574154672         4.5e-09   True
  20    1.66e-05   0.574154668         5.4e-11   True
  24    1.80e-06   0.574154668         2.1e-11   True
  28    1.96e-07   0.574154668         4.7e-11   True
  32    2.12e-08   0.574154670         2.0e-09   True
  36    2.31e-09   0.574154665         2.3e-09   True
  40    2.51e-10   0.574154720         5.2e-08   True
----------------------------------------------------------
worst gap inside the band = 1.1e-05   (tolerance 1e-4)
check passed: the error contracts at |lambda_2 / lambda_1|.

rng = np.random.default_rng(5)

I = J = 4
K = L = 3
A_i_j = rng.normal(size=(I, J))
B_k_l = rng.normal(size=(K, L))
C_i_l = rng.normal(size=(I, L))

# route 1: vectorize with the HOUSE identity (12.2), vec_C(A X B) = (A (x) B') vec_C(X)
Kron_il_jk = np.kron(A_i_j, B_k_l.T)                  # (I*L, J*K) = (12, 12)
x_jk = np.linalg.solve(Kron_il_jk, C_i_l.reshape(-1))
X_kron_j_k = x_jk.reshape(J, K)

# route 2: no Kronecker product at all -- two ordinary solves
#   A X B = C  =>  X B = A^{-1} C  =>  B' X' = (A^{-1} C)'
X_direct_j_k = np.linalg.solve(B_k_l.T, np.linalg.solve(A_i_j, C_i_l).T).T

gap = np.abs(X_kron_j_k - X_direct_j_k).max()
resid = np.abs(A_i_j @ X_direct_j_k @ B_k_l - C_i_l).max()

print(f"Kronecker matrix shape: {Kron_il_jk.shape}")
print(f"max |route 1 - route 2| = {gap:.2e}   (tolerance 1e-9)")
print(f"max |A X B - C|         = {resid:.2e}   (tolerance 1e-10)")
assert gap < 1e-9 and resid < 1e-10
print("check passed: both routes solve the matrix equation.")

Kronecker matrix shape: (12, 12)
max |route 1 - route 2| = 3.55e-14   (tolerance 1e-9)
max |A X B - C|         = 2.66e-15   (tolerance 1e-10)
check passed: both routes solve the matrix equation.

# the cost, at the n = 60 of section 12
n = 60
kron_side = n * n                                  # the Kronecker system is n^2 x n^2
flops_kron = (2 / 3) * kron_side ** 3              # LU on an (n^2 x n^2) matrix
flops_direct = 2 * (2 / 3) * n ** 3                # two LU factorizations of n x n

print(f"n = {n}:")
print(f"  Kronecker matrix : {kron_side} x {kron_side} = {kron_side ** 2 / 1e6:.1f}M entries, "
      f"{kron_side ** 2 * 8 / 1e9:.2f} GB in float64")
print(f"  solve on it      : ~{flops_kron:.3e} flops")
print(f"  two n x n solves : ~{flops_direct:.3e} flops")
print(f"  ratio            : {flops_kron / flops_direct:.3e}x")

n = 60:
  Kronecker matrix : 3600 x 3600 = 13.0M entries, 0.10 GB in float64
  solve on it      : ~3.110e+10 flops
  two n x n solves : ~2.880e+05 flops
  ratio            : 1.080e+05x

rng = np.random.default_rng(5)

# the exact standard deviation of the Monte Carlo integrand 4 * 1{inside}
p_in = math.pi / 4
sigma_f = 4 * math.sqrt(p_in * (1 - p_in))

print(f"{'N':>10}{'Leibniz err':>14}{'N x err':>10}"
      f"{'MC err':>12}{'MC std err':>13}{'|err| / se':>12}")
print("-" * 71)
for N_ex in (10 ** 3, 10 ** 5, 10 ** 7):
    k_n = np.arange(N_ex)
    err_leib = abs(4 * ((-1.0) ** k_n / (2 * k_n + 1)).sum() - math.pi)

    U_n_k = rng.random(size=(N_ex, 2))
    err_mc = abs(4 * ((U_n_k ** 2).sum(axis=1) < 1.0).mean() - math.pi)
    se_mc = sigma_f / math.sqrt(N_ex)

    print(f"{N_ex:>10}{err_leib:>14.2e}{N_ex * err_leib:>10.4f}"
          f"{err_mc:>12.2e}{se_mc:>13.2e}{err_mc / se_mc:>12.2f}")

print("-" * 71)
print("Leibniz: N x error is exactly ~1, so error = 1/N.")
print("Monte Carlo: each error is O(1) standard errors, and se = sigma_f / sqrt(N).")

         N   Leibniz err   N x err      MC err   MC std err  |err| / se
-----------------------------------------------------------------------
      1000      1.00e-03    1.0000    8.24e-02     5.19e-02        1.59
    100000      1.00e-05    1.0000    1.87e-03     5.19e-03        0.36

  10000000      1.00e-07    1.0000    7.83e-05     5.19e-04        0.15
-----------------------------------------------------------------------
Leibniz: N x error is exactly ~1, so error = 1/N.
Monte Carlo: each error is O(1) standard errors, and se = sigma_f / sqrt(N).

mu_c, sigma_c = 1.0, 0.30
a_grid = np.array([0.5, 1.0, 2.0, 4.0])
N_ex = 100_000

def premium(a, Z_n):
    """Risk premium implied by draws Z_n of a standard normal."""
    u_n = -np.exp(-a * (mu_c + sigma_c * Z_n))
    return mu_c + math.log(-u_n.mean()) / a

rng = np.random.default_rng(5)
fresh = [premium(a, rng.standard_normal(N_ex)) for a in a_grid]      # new draws each time

rng = np.random.default_rng(5)
Z_common_n = rng.standard_normal(N_ex)                               # ONE set of draws
common = [premium(a, Z_common_n) for a in a_grid]

exact = 0.5 * a_grid * sigma_c ** 2

print(f"{'a':>6}{'exact':>12}{'fresh draws':>15}{'common draws':>15}")
print("-" * 48)
for a, ex, fr, co in zip(a_grid, exact, fresh, common):
    print(f"{a:>6.1f}{ex:>12.6f}{fr:>15.6f}{co:>15.6f}")
print("-" * 48)
print(f"max |fresh  - exact| = {np.abs(np.array(fresh) - exact).max():.2e}")
print(f"max |common - exact| = {np.abs(np.array(common) - exact).max():.2e}")

     a       exact    fresh draws   common draws
------------------------------------------------
   0.5    0.022500       0.022167       0.022167
   1.0    0.045000       0.044489       0.044592
   2.0    0.090000       0.090534       0.089377
   4.0    0.180000       0.178464       0.178481
------------------------------------------------
max |fresh  - exact| = 1.54e-03
max |common - exact| = 1.52e-03

# the comparative static: pi(4) - pi(0.5), estimated 200 times under each scheme
R_rep = 200
rng = np.random.default_rng(11)

diff_fresh = np.empty(R_rep)
diff_common = np.empty(R_rep)
for r in range(R_rep):
    diff_fresh[r] = premium(4.0, rng.standard_normal(N_ex)) - premium(0.5, rng.standard_normal(N_ex))
    Z_n = rng.standard_normal(N_ex)
    diff_common[r] = premium(4.0, Z_n) - premium(0.5, Z_n)

exact_diff = 0.5 * (4.0 - 0.5) * sigma_c ** 2

print(f"exact difference pi(4) - pi(0.5) = {exact_diff:.6f}\n")
print(f"{'scheme':<20}{'mean':>12}{'std dev':>12}{'bias':>12}")
print("-" * 56)
print(f"{'fresh draws':<20}{diff_fresh.mean():>12.6f}{diff_fresh.std(ddof=1):>12.6f}"
      f"{diff_fresh.mean() - exact_diff:>12.2e}")
print(f"{'common draws':<20}{diff_common.mean():>12.6f}{diff_common.std(ddof=1):>12.6f}"
      f"{diff_common.mean() - exact_diff:>12.2e}")
print("-" * 56)
ratio = (diff_fresh.std(ddof=1) / diff_common.std(ddof=1)) ** 2
print(f"variance ratio (fresh / common) = {ratio:.1f}x")
assert ratio > 1.5
print("check passed: common random numbers sharpen the comparative static.")

exact difference pi(4) - pi(0.5) = 0.157500

scheme                      mean     std dev        bias
--------------------------------------------------------
fresh draws             0.157456    0.001730   -4.38e-05
common draws            0.157465    0.001012   -3.46e-05
--------------------------------------------------------
variance ratio (fresh / common) = 2.9x
check passed: common random numbers sharpen the comparative static.

rng = np.random.default_rng(5)
N_ex = 100_000
n_n = np.arange(1, N_ex + 1)

x_cauchy_n = rng.standard_cauchy(size=N_ex)
run_mean_n = np.cumsum(x_cauchy_n) / n_n

fig, axes = plt.subplots(1, 2, figsize=(12, 3.6))
axes[0].plot(run_mean_n, lw=0.8)
axes[0].axhline(0.0, color="k", ls="--", lw=0.8)
axes[0].set_title("Cauchy: the running MEAN never settles")
axes[0].set_xlabel("n"); axes[0].set_ylabel("running mean")

run_median_n = np.array([np.median(x_cauchy_n[:m]) for m in
                         np.linspace(10, N_ex, 300).astype(int)])
axes[1].plot(np.linspace(10, N_ex, 300), run_median_n, lw=1, color="C1")
axes[1].axhline(0.0, color="k", ls="--", lw=0.8)
axes[1].set_title("Cauchy: the running MEDIAN converges to 0")
axes[1].set_xlabel("n"); axes[1].set_ylabel("running median")
plt.tight_layout(); plt.show()

print(f"largest |x| drawn: {np.abs(x_cauchy_n).max():,.0f}")
print(f"running mean at n = 1000, 10000, 100000: "
      f"{run_mean_n[999]:.3f}, {run_mean_n[9999]:.3f}, {run_mean_n[-1]:.3f}")

largest |x| drawn: 38,457
running mean at n = 1000, 10000, 100000: -38.417, -3.997, -0.489

# the "standard error", computed as if the CLT applied, across three seeds
print(f"{'seed':>6}{'sample mean':>16}{'sample sd':>16}{'std error (fake)':>18}")
print("-" * 54)
for seed in (1, 2, 3):
    x_n = np.random.default_rng(seed).standard_cauchy(size=N_ex)
    se_fake = x_n.std(ddof=1) / math.sqrt(N_ex)
    print(f"{seed:>6}{x_n.mean():>16.4f}{x_n.std(ddof=1):>16.4f}{se_fake:>16.4f}")
print("-" * 54)
sds = [float(np.random.default_rng(sd).standard_cauchy(size=N_ex).std(ddof=1)) for sd in (1, 2, 3)]
print(f"the sample SD itself varies by a factor of {max(sds) / min(sds):.1f} across seeds --")
print("it is not converging to anything, because the variance is infinite.")

# the median, by contrast, is stable across seeds
medians = [float(np.median(np.random.default_rng(seed).standard_cauchy(size=N_ex)))
           for seed in (1, 2, 3)]
print(f"\nsample medians across seeds: {[round(m, 4) for m in medians]}")
assert max(abs(m) for m in medians) < 0.05
print("check passed: the median is stable at 0; the mean is not stable at anything.")

  seed     sample mean       sample sd  std error (fake)
------------------------------------------------------
     1         -0.0155        135.5096          0.4285
     2          0.6581        222.7595          0.7044
     3         -1.8533        620.4898          1.9622
------------------------------------------------------
the sample SD itself varies by a factor of 4.6 across seeds --
it is not converging to anything, because the variance is infinite.

sample medians across seeds: [-0.0076, -0.005, -0.005]
check passed: the median is stable at 0; the mean is not stable at anything.

name	shape	meaning
`x_i`	$(I,)$	a vector indexed by $i$
`M_i_j`	$(I, J)$	a matrix with axes $(i, j)$
`P_x_y_xp`	$(n_x, n_y, n_x)$	transition array with axes $(x, y, x')$; `xp` denotes $x'$
`mu_xy`	$(n_x n_y,)$	the row-major flattening of `mu_x_y`
`mean_i_1`	$(I, 1)$	a reduction with `keepdims=True`; the `1` marks the kept axis

dtype	bytes	range / precision
`int32`, `int64`	4, 8	$\pm 2^{31}$, $\pm 2^{63}$, exact
`float32`, `float64`	4, 8	IEEE 754 single, double
`bool`	1	`True` / `False`
`complex64`, `complex128`	8, 16	two floats

Operation	NumPy
Matrix product $AB$	`A @ B`
Solve $Ax = b$	`np.linalg.solve(A, b)`
Eigendecomposition	`np.linalg.eig` (general), `eigh` (symmetric)
SVD $A = U\Sigma V^\top$	`np.linalg.svd`
Least squares	`np.linalg.lstsq`
Determinant, trace, norm	`np.linalg.det`, `np.trace`, `np.linalg.norm`
Matrix inverse	don't: see below

draw	call
uniform on $[0,1)$	`rng.random(size)`
normal $\mathcal{N}(\mu,\sigma^2)$	`rng.normal(loc, scale, size)`
multivariate normal	`rng.multivariate_normal(mean, cov, size)`
exponential, mean $\theta$	`rng.exponential(scale, size)`
discrete uniform	`rng.integers(low, high, size)`
a categorical draw	`rng.choice(a, size, p=...)`
a random permutation	`rng.permutation(n)`

NumPy: arrays, vectorization, and simulation¶

Alfred Galichon (NYU)¶

'math+econ+code' masterclass series: fundamentals of research in python¶

With python code examples¶

Learning objectives¶

References¶

1. Motivation¶

Part I: arrays and vectorization¶

2. Creating arrays¶

3. Attributes, dtypes, and the naming convention¶

4. Indexing and slicing¶

5. Fancy indexing¶

6. Universal functions¶

7. Broadcasting¶

8. Reductions and the axis keyword¶

9. Linear algebra¶

10. Worked example: a bond portfolio, vectorized¶

11. Worked example: OLS by hand¶

12. The Kronecker product, and two vectorization conventions¶

Two conventions, two identities¶

13. Worked example: assembling an optimal-transport constraint matrix¶

Part II: randomness and simulation¶

14. Pseudorandom numbers, and the Generator API¶

15. Sampling from distributions¶

16. Reproducibility, in practice¶

17. Monte Carlo estimation, and the $N^{-1/2}$ rate¶

18. Worked example: $\pi$ by Monte Carlo¶

19. The law of large numbers, empirically¶

20. The central limit theorem, empirically¶

21. Reporting Monte Carlo: standard errors and confidence intervals¶

22. Variance reduction: antithetic variates¶

23. Worked example: the risk premium of an uncertain income¶

24. Summary¶

25. Exercises¶

26. Further directions¶

27. Solutions to the exercises¶

Solution to Exercise 1: Power iteration, and its rate¶

Solution to Exercise 2: The vec trick at work¶

Solution to Exercise 3: Leibniz against Monte Carlo¶

Solution to Exercise 4: Risk aversion, and common random numbers¶

Solution to Exercise 5: When the law of large numbers fails¶

8. Reductions and the `axis` keyword¶

14. Pseudorandom numbers, and the `Generator` API¶