"""CSE 541, Autumn 2026 -- Homework 1 starter code (Problems 1 and 3). Fill in the functions marked TODO. Everything else is provided. Run `python hw1_starter.py` to produce the plots for 1.2 and 3.1-3.2 once the TODOs are done. Conventions ----------- R : np.ndarray of shape (T, d), gross returns r_t(i) > 0 an algorithm returns P : np.ndarray of shape (T, d), the portfolio p_t played on day t (row t must be computed from R[:t] only) comparator u : np.ndarray of shape (T,), index of the stock the comparator holds on day t """ import csv import numpy as np import matplotlib.pyplot as plt rng = np.random.default_rng(541) # ---------------------------------------------------------------- data ---- def load_stock_returns(path="stock_closes.csv"): """Daily closes -> gross returns r_t(i) = close_{t+1}(i) / close_t(i). Returns (tickers, R) with R of shape (T, d), entries > 0.""" with open(path) as f: rows = list(csv.reader(f)) tickers = rows[0][1:] closes = np.array([[float(x) for x in r[1:]] for r in rows[1:]]) return tickers, closes[1:] / closes[:-1] def synthetic_market(T=5000, d=8, S=3, drift=0.02, vol=0.02, seed=0): """Regime-switching market of gross returns. In each of S+1 regimes one stock is 'good': its daily log-return has mean +drift, all others have mean -drift/2; all have volatility `vol`. Returns (R, u) where R has shape (T, d) and u is the comparator holding the good stock in each regime (S switches).""" g = np.random.default_rng(seed) boundaries = np.linspace(0, T, S + 2).astype(int) logR = np.empty((T, d)) u = np.empty(T, dtype=int) good = -1 for k in range(S + 1): good = g.choice([i for i in range(d) if i != good]) lo, hi = boundaries[k], boundaries[k + 1] means = np.full(d, -drift / 2) means[good] = drift logR[lo:hi] = means + vol * g.normal(size=(hi - lo, d)) u[lo:hi] = good return np.exp(logR), u def piecewise_best_comparator(R, block): """Comparator holding the stock with the largest log-wealth within each block of `block` consecutive days.""" T = R.shape[0] u = np.empty(T, dtype=int) logR = np.log(R) for lo in range(0, T, block): hi = min(lo + block, T) u[lo:hi] = np.argmax(logR[lo:hi].sum(axis=0)) return u def num_switches(u): return int(np.sum(u[1:] != u[:-1])) # ---------------------------------------------------------- evaluation ---- def log_wealth(P, R): """Cumulative log-wealth sum_{s<=t} log of a portfolio sequence.""" return np.cumsum(np.log(np.sum(P * R, axis=1))) def comparator_log_wealth(R, u): return np.cumsum(np.log(R[np.arange(len(u)), u])) def dynamic_regret(P, R, u): """sum_t log r_t(u_t) - sum_t log (log-wealth of the switching comparator minus log-wealth of the algorithm).""" return comparator_log_wealth(R, u)[-1] - log_wealth(P, R)[-1] def fixed_share_bound(T, d, S, alpha=None): """The theorem's bound on the dynamic regret: part (i) for a constant alpha, part (ii) (alpha=None) for alpha_t = 1/t.""" if alpha is None: return (S + 1) * np.log(d) + (S + 1) * np.log(T) return (S + 1) * np.log(d) + S * np.log(1 / alpha) + (T - 1 - S) * np.log(1 / (1 - alpha)) # ---------------------------------------------------------- algorithms ---- def buy_and_hold(R): """Provided. Exponential weights with eta = 1 and losses -log r_t(i): w_{t+1}(i) = w_t(i) r_t(i), p_t = w_t / sum(w_t). Returns P of shape (T, d).""" T, d = R.shape P = np.empty((T, d)) w = np.ones(d) for t in range(T): P[t] = w / w.sum() w = w * R[t] w /= w.sum() # rescaling never changes p_t return P def fixed_share(R, alpha): """TODO (1.2): fixed share with log loss and eta = 1, in lecture's unnormalized notation: w_1 = 1; play p_t = w_t / sum(w_t); w~_{t+1} = w_t * r_t; w_{t+1} = (1-alpha) w~_{t+1} + (alpha/d) sum(w~_{t+1}). (Rescale w each day to avoid overflow; only p_t matters.) `alpha` is a float, or a function of the 1-indexed day, e.g. lambda t: 1/t, in which case alpha(t+1) is used when forming w_{t+1} from w_t (so the mixing applied after observing r_1 uses alpha(2), etc.).""" T, d = R.shape P = np.empty((T, d)) raise NotImplementedError def fund_of_funds(R, alpha): """TODO (1.1/1.2): keep one buy-and-hold fund B_s opened on every day s (portfolio p_t^{(s)} proportional to prod_{tau=s}^{t-1} r_tau, daily growth Z_t^{(s)} = ) and wealth shares nu_t over funds with nu_1(1) = 1 and nu_{t+1}(s) = (1-a) nu_t(s) Z_t^{(s)} / Z_t (s <= t), nu_{t+1}(t+1) = a, where Z_t = sum_s nu_t(s) Z_t^{(s)} and a = alpha(t+1) if alpha is callable else alpha. Play sum_s nu_t(s) p_t^{(s)}. O(d t) per day is fine.""" T, d = R.shape P = np.empty((T, d)) raise NotImplementedError # --------------------------------------------------------------- driver ---- def run_all(R, u, title, ax): T, d = R.shape S = num_switches(u) alpha_t = lambda t: 1.0 / t # parameter-free (theorem part ii) alpha_tuned = max(S, 1) / (T - 1) # needs S (theorem part i) algs = { "buy and hold": lambda: buy_and_hold(R), "fixed share, alpha_t = 1/t": lambda: fixed_share(R, alpha_t), "fund of funds, alpha_t = 1/t": lambda: fund_of_funds(R, alpha_t), f"fixed share, alpha = S/(T-1) = {alpha_tuned:.2g}": lambda: fixed_share(R, alpha_tuned), } results = {} stocks = np.cumsum(np.log(R), axis=0) # log-wealth of each single stock ax.plot(stocks, color="0.6", alpha=0.35, lw=0.8) ax.plot([], [], color="0.6", alpha=0.6, lw=0.8, label="individual stocks") ax.plot(comparator_log_wealth(R, u), "k--", label=f"comparator (S={S})") for name, fn in algs.items(): try: P = fn() except NotImplementedError: print(f" [{title}] {name}: not implemented yet") continue results[name] = P ax.plot(log_wealth(P, R), label=name) if "buy and hold" in name: print(f" [{title}] {name:38s} D-Reg_T = {dynamic_regret(P, R, u):8.3f}") else: bound = fixed_share_bound(T, d, S, None if "1/t" in name else alpha_tuned) print(f" [{title}] {name:38s} D-Reg_T = {dynamic_regret(P, R, u):8.3f} (bound {bound:.1f})") keys = [k for k in results if "1/t" in k] if len(keys) == 2: diff = np.max(np.abs(results[keys[0]] - results[keys[1]])) print(f" [{title}] max |p_t(fixed share) - p_t(fund of funds)| = {diff:.2e}") ax.set_title(title) ax.set_xlabel("day t") ax.set_ylabel("cumulative log-wealth") ax.legend(fontsize=8) # ===================================================================== # Problem 3: ERM vs online gradient descent # ===================================================================== def generate_stream(T=5000, d=10, mode="stationary", sigma=0.1, revolutions=3, seed=0): """Streaming linear regression. x_t uniform on the unit sphere (so E[x x^T] = I/d), y_t = + N(0, sigma^2), ||w*_t|| = 1. mode: 'stationary' (w*_t fixed), 'drift' (w*_t rotates in the (e1,e2) plane, `revolutions` full turns over T rounds), 'switch' (w* -> -w* at T/2). Returns X (T,d), Y (T,), Wstar (T,d).""" g = np.random.default_rng(seed) X = g.normal(size=(T, d)); X /= np.linalg.norm(X, axis=1, keepdims=True) w0 = g.normal(size=d); w0 /= np.linalg.norm(w0) Wstar = np.tile(w0, (T, 1)) if mode == "drift": theta = 2 * np.pi * revolutions * np.arange(T) / T Wstar = np.zeros((T, d)); Wstar[:, 0] = np.cos(theta); Wstar[:, 1] = np.sin(theta) elif mode == "switch": Wstar[T // 2:] *= -1 Y = np.sum(Wstar * X, axis=1) + sigma * g.normal(size=T) return X, Y, Wstar def excess_risk(W, wstar): """Population excess risk (1/2d)||w - w*||^2 for each row of W, for a fixed w*.""" d = W.shape[1] return np.sum((W - wstar) ** 2, axis=1) / (2 * d) def dynamic_regret_curve(W, X, Y, Wstar): """cumulative sum of l_t(w_t) - l_t(w*_t), with l_t(w) = (1/2)( - y_t)^2.""" alg = 0.5 * (np.sum(W * X, axis=1) - Y) ** 2 comp = 0.5 * (np.sum(Wstar * X, axis=1) - Y) ** 2 return np.cumsum(alg - comp) def erm_ftl(X, Y, rho=1e-3): """TODO (3.1): w_t = argmin_w sum_{s-y_s)^2 + (rho/2)||w||^2, re-solved every round (w_1 = 0). Return W of shape (T, d) with row t = w_t.""" raise NotImplementedError def ogd(X, Y, eta_fn, B=2.0): """TODO (3.1): projected OGD, w_{t+1} = Pi_{||w||<=B}(w_t - eta_fn(t) * g_t), g_t = ( - y_t) x_t, w_1 = 0, with eta_fn(t) the step on round t (1-indexed). Return (W, Wbar): the iterates and their running averages.""" raise NotImplementedError def windowed_erm(X, Y, W_len, rho=1e-3): """TODO (3.3): least squares on the most recent W_len rounds.""" raise NotImplementedError def problem3_figure(T=5000, d=10, sigma=0.1, seeds=range(10)): fig, axes = plt.subplots(1, 3, figsize=(15, 4.2)) methods = { "ERM": lambda X, Y: erm_ftl(X, Y), "OGD eta=d/t": lambda X, Y: ogd(X, Y, lambda t: d / t)[0], "OGD eta=d/t (avg)": lambda X, Y: ogd(X, Y, lambda t: d / t)[1], "OGD eta=0.05": lambda X, Y: ogd(X, Y, lambda t: 0.05)[0], "OGD eta=0.2": lambda X, Y: ogd(X, Y, lambda t: 0.2)[0], "OGD eta=0.5": lambda X, Y: ogd(X, Y, lambda t: 0.5)[0], } for ax, mode in zip(axes, ("stationary", "drift", "switch")): for name, fn in methods.items(): curves = [] for seed in seeds: X, Y, Wstar = generate_stream(T, d, mode, sigma, seed=seed) try: Wt = fn(X, Y) except NotImplementedError: Wt = None break curves.append(excess_risk(Wt, Wstar[0]) if mode == "stationary" else dynamic_regret_curve(Wt, X, Y, Wstar)) if not curves: print(f" [{mode}] {name}: not implemented yet"); continue ax.plot(np.mean(curves, axis=0), label=name) ax.set_title(mode); ax.set_xlabel("round t") if mode == "stationary": ax.set_xscale("log"); ax.set_yscale("log"); ax.set_ylabel("excess risk") else: ax.set_ylabel("cumulative dynamic regret") ax.legend(fontsize=7) fig.tight_layout(); fig.savefig("hw1_problem3.pdf"); print("wrote hw1_problem3.pdf") if __name__ == "__main__": # ---- Problem 1 ---- fig, axes = plt.subplots(2, 2, figsize=(11, 7.5)) for ax, S in zip(axes.flat[:3], (1, 3, 10)): R, u = synthetic_market(T=5000, d=8, S=S, seed=S) run_all(R, u, f"synthetic market, S = {S}", ax) tickers, R = load_stock_returns() u = piecewise_best_comparator(R, block=50) run_all(R, u, "stocks, best per 50 days", axes.flat[3]) fig.tight_layout() fig.savefig("hw1_problem1.pdf") print("wrote hw1_problem1.pdf") # ---- Problem 3 ---- problem3_figure()