Use vectorization and numba jitting for optimization

Author: kp992

lectures/black_litterman.md

@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -79,15 +81,14 @@ adjust investors' subjective beliefs about mean returns in order to render more
 
 Let's start with some imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 import numpy as np
 import scipy.stats as stat
 import matplotlib.pyplot as plt
+from numba import jit
 from ipywidgets import interact, FloatSlider
 ```
 
-
-
 ## Mean-Variance Portfolio Choice
 
 A risk-free security earns one-period net return $r_f$.
@@ -169,7 +170,7 @@ $w$'s with **extreme long and short positions**.
 A common reaction to these outcomes is that they are so implausible that a portfolio
 manager cannot recommend them to a customer.
 
-```{code-cell} python3
+```{code-cell} ipython3
 np.random.seed(12)
 
 N = 10                                           # Number of assets
@@ -300,7 +301,7 @@ The starting point of the Black-Litterman portfolio choice model is thus
 a pair $(\delta_m, \mu_m)$ that tells the customer to hold the
 market portfolio.
 
-```{code-cell} python3
+```{code-cell} ipython3
 # Observed mean excess market return
 r_m = w_m @ μ_est
 
@@ -316,11 +317,12 @@ d_m = r_m / σ_m
 # Derive "view" which would induce the market portfolio
 μ_m = (d_m * Σ_est @ w_m).reshape(N, 1)
 
+x = np.arange(N) + 1
 fig, ax = plt.subplots(figsize=(8, 5))
 ax.set_title(r'Difference between $\hat{\mu}$ (estimate) and $\mu_{BL}$ (market implied)')
-ax.plot(np.arange(N)+1, μ_est, 'o', c='k', label='$\hat{\mu}$')
-ax.plot(np.arange(N)+1, μ_m, 'o', c='r', label='$\mu_{BL}$')
-ax.vlines(np.arange(N) + 1, μ_m, μ_est, lw=1)
+ax.plot(x, μ_est, 'o', c='k', label='$\hat{\mu}$')
+ax.plot(x, μ_m, 'o', c='r', label='$\mu_{BL}$')
+ax.vlines(x, μ_m, μ_est, lw=1)
 ax.axhline(0, c='k', ls='--')
 ax.set_xlabel('Assets')
 ax.xaxis.set_ticks(np.arange(1, N+1, 1))
@@ -384,7 +386,7 @@ If $\hat \mu$ is the maximum likelihood estimator
 and $\tau$ is chosen heavily to weight this view, then the
 customer's portfolio will involve big short-long positions.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def black_litterman(λ, μ1, μ2, Σ1, Σ2):
     """
     This function calculates the Black-Litterman mixture
@@ -607,7 +609,7 @@ $\bar d_2$ on the RHS of the constraint, by varying
 $\bar d_2$ (or $\lambda$ ), we can trace out the whole curve
 as the figure below illustrates.
 
-```{code-cell} python3
+```{code-cell} ipython3
 np.random.seed(1987102)
 
 N = 2                                           # Number of assets
@@ -649,13 +651,9 @@ def decolletage(λ):
     dist_r_hat = stat.multivariate_normal(μ_est.squeeze(), τ * Σ_est)
 
     X, Y = np.meshgrid(r1, r2)
-    Z_BL = np.zeros((N_r1, N_r2))
-    Z_hat = np.zeros((N_r1, N_r2))
-
-    for i in range(N_r1):
-        for j in range(N_r2):
-            Z_BL[i, j] = dist_r_BL.pdf(np.hstack([X[i, j], Y[i, j]]))
-            Z_hat[i, j] = dist_r_hat.pdf(np.hstack([X[i, j], Y[i, j]]))
+    XY = np.stack((X, Y), axis=-1)
+    Z_BL = dist_r_BL.pdf(XY) 
+    Z_hat = dist_r_hat.pdf(XY)
 
     μ_tilde = black_litterman(λ, μ_m, μ_est, Σ_est, τ * Σ_est).flatten()
 
@@ -692,7 +690,7 @@ This leads to the
 following figure, on which the curve connecting $\hat \mu$
 and $\mu_{BL}$ are bending
 
-```{code-cell} python3
+```{code-cell} ipython3
 λ_grid = np.linspace(.001, 20000, 1000)
 curve = np.asarray([black_litterman(λ, μ_m, μ_est, Σ_est,
                                     τ * np.eye(N)).flatten() for λ in λ_grid])
@@ -705,13 +703,9 @@ def decolletage(λ):
     dist_r_hat = stat.multivariate_normal(μ_est.squeeze(), τ * np.eye(N))
 
     X, Y = np.meshgrid(r1, r2)
-    Z_BL = np.zeros((N_r1, N_r2))
-    Z_hat = np.zeros((N_r1, N_r2))
-
-    for i in range(N_r1):
-        for j in range(N_r2):
-            Z_BL[i, j] = dist_r_BL.pdf(np.hstack([X[i, j], Y[i, j]]))
-            Z_hat[i, j] = dist_r_hat.pdf(np.hstack([X[i, j], Y[i, j]]))
+    XY = np.stack((X, Y), axis=-1)
+    Z_BL = dist_r_BL.pdf(XY) 
+    Z_hat = dist_r_hat.pdf(XY)
 
     μ_tilde = black_litterman(λ, μ_m, μ_est, Σ_est, τ * np.eye(N)).flatten()
 
@@ -1247,7 +1241,7 @@ observations is related to the sampling frequency
 
 - Moreover, for a fixed lag length, $n$, the dependence vanishes as the sampling frequency goes to infinity. In fact, letting $h$ go to $\infty$ gives back the case of IID data.
 
-```{code-cell} python3
+```{code-cell} ipython3
 μ = .0
 κ = .1
 σ = .5
@@ -1346,7 +1340,8 @@ thus getting an idea about how the asymptotic relative MSEs changes in
 the sampling frequency $h$ relative to the IID case that we
 compute in closed form.
 
-```{code-cell} python3
+```{code-cell} ipython3
+@jit
 def sample_generator(h, N, M):
     ϕ = (1 - np.exp(-κ * h)) * μ
     ρ = np.exp(-κ * h)
@@ -1355,18 +1350,18 @@ def sample_generator(h, N, M):
     mean_uncond = μ
     std_uncond = np.sqrt(σ**2 / (2 * κ))
 
-    ε_path = stat.norm(0, np.sqrt(s)).rvs((M, N))
+    ε_path = np.random.normal(0, np.sqrt(s), (M, N))
 
     y_path = np.zeros((M, N + 1))
-    y_path[:, 0] = stat.norm(mean_uncond, std_uncond).rvs(M)
+    y_path[:, 0] = np.random.normal(mean_uncond, std_uncond, M)
 
     for i in range(N):
         y_path[:, i + 1] = ϕ + ρ * y_path[:, i] + ε_path[:, i]
 
     return y_path
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 # Generate large sample for different frequencies
 N_app, M_app = 1000, 30000        # Sample size, number of simulations
 h_grid = np.linspace(.1, 80, 30)