black_litterman.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ae8d3af1",
   "metadata": {},
   "source": [
    "\n",
    "<a id='black-litterman'></a>\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "568e4ef6",
   "metadata": {},
   "source": [
    "# Two Modifications of Mean-Variance Portfolio Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7b12c76",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture describes extensions to the classical mean-variance portfolio theory summarized in our lecture [Elementary Asset Pricing Theory](https://python-advanced.quantecon.org/asset_pricing_lph.html).\n",
    "\n",
    "The classic theory described there assumes that a decision maker completely trusts the statistical model that he posits to govern the joint distribution of returns on a list of available assets.\n",
    "\n",
    "Both extensions described here put distrust of that statistical model into the mind of the decision maker.\n",
    "\n",
    "One is a model of Black and Litterman [[Black and Litterman, 1992](https://python-advanced.quantecon.org/zreferences.html#id68)] that imputes to the decision maker distrust of historically estimated mean returns but still complete trust of estimated covariances of returns.\n",
    "\n",
    "The second model also imputes to the decision maker doubts about his  statistical model, but now by saying that, because of that distrust,  the decision maker uses a version of  robust control theory described in this lecture [Robustness](https://python-advanced.quantecon.org/robustness.html).\n",
    "\n",
    "The famous **Black-Litterman** (1992) [[Black and Litterman, 1992](https://python-advanced.quantecon.org/zreferences.html#id68)] portfolio choice model was motivated by the finding that with high frequency or\n",
    "moderately high  frequency data, means are more difficult to estimate than\n",
    "variances.\n",
    "\n",
    "A model of **robust portfolio choice** that we’ll describe below also begins\n",
    "from the same starting point.\n",
    "\n",
    "To begin, we’ll take for granted that means are more difficult to\n",
    "estimate that covariances and will focus on how Black and Litterman, on\n",
    "the one hand, an robust control theorists, on the other, would recommend\n",
    "modifying the **mean-variance portfolio choice model** to take that into\n",
    "account.\n",
    "\n",
    "At the end of this lecture, we shall use some rates of convergence\n",
    "results and some simulations to verify how means are more difficult to\n",
    "estimate than variances.\n",
    "\n",
    "Among the ideas in play in this lecture will be\n",
    "\n",
    "- Mean-variance portfolio theory  \n",
    "- Bayesian approaches to estimating linear regressions  \n",
    "- A risk-sensitivity operator and its connection to robust control\n",
    "  theory  \n",
    "\n",
    "\n",
    "In summary, we’ll describe two ways to  modify the classic\n",
    "mean-variance portfolio choice model in ways designed to make its\n",
    "recommendations more plausible.\n",
    "\n",
    "Both of the adjustments that we describe are designed to confront a\n",
    "widely recognized embarrassment to mean-variance portfolio theory,\n",
    "namely, that it usually implies taking very extreme long-short portfolio\n",
    "positions.\n",
    "\n",
    "The two approaches build on a common and widespread hunch –\n",
    "that because it is much easier statistically to estimate covariances of\n",
    "excess returns than it is to estimate their means, it makes sense to\n",
    "adjust investors’ subjective beliefs about mean returns in order to render more plausible decisions.\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08052653",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stat\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact, FloatSlider"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "422026f9",
   "metadata": {},
   "source": [
    "## Mean-Variance Portfolio Choice\n",
    "\n",
    "A risk-free security earns one-period net return $ r_f $.\n",
    "\n",
    "An $ n \\times 1 $ vector of risky securities earns\n",
    "an $ n \\times 1 $ vector $ \\vec r - r_f {\\bf 1} $ of *excess\n",
    "returns*, where $ {\\bf 1} $ is an $ n \\times 1 $ vector of\n",
    "ones.\n",
    "\n",
    "The excess return vector is multivariate normal with mean $ \\mu $\n",
    "and covariance matrix $ \\Sigma $, which we express either as\n",
    "\n",
    "$$\n",
    "\\vec r - r_f {\\bf 1} \\sim {\\mathcal N}(\\mu, \\Sigma)\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "\\vec r - r_f {\\bf 1} = \\mu + C \\epsilon\n",
    "$$\n",
    "\n",
    "where $ \\epsilon \\sim {\\mathcal N}(0, I) $ is an $ n \\times 1 $\n",
    "random vector.\n",
    "\n",
    "Let $ w $ be an $ n \\times 1 $  vector of portfolio weights.\n",
    "\n",
    "A portfolio consisting $ w $ earns returns\n",
    "\n",
    "$$\n",
    "w' (\\vec r - r_f {\\bf 1}) \\sim {\\mathcal N}(w' \\mu, w' \\Sigma w )\n",
    "$$\n",
    "\n",
    "The **mean-variance portfolio choice problem** is to choose $ w $ to\n",
    "maximize\n",
    "\n",
    "\n",
    "<a id='equation-choice-problem'></a>\n",
    "$$\n",
    "U(\\mu,\\Sigma;w) = w'\\mu - \\frac{\\delta}{2} w' \\Sigma w \\tag{36.1}\n",
    "$$\n",
    "\n",
    "where $ \\delta > 0 $ is a risk-aversion parameter. The first-order\n",
    "condition for maximizing [(36.1)](#equation-choice-problem) with respect to the vector $ w $ is\n",
    "\n",
    "$$\n",
    "\\mu = \\delta \\Sigma w\n",
    "$$\n",
    "\n",
    "which implies the following design of a risky portfolio:\n",
    "\n",
    "\n",
    "<a id='equation-risky-portfolio'></a>\n",
    "$$\n",
    "w = (\\delta \\Sigma)^{-1} \\mu \\tag{36.2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76f084cf",
   "metadata": {},
   "source": [
    "## Estimating  Mean and Variance\n",
    "\n",
    "The key inputs into the portfolio choice model [(36.2)](#equation-risky-portfolio) are\n",
    "\n",
    "- estimates of the parameters $ \\mu, \\Sigma $ of the random excess\n",
    "  return vector$ (\\vec r - r_f {\\bf 1}) $  \n",
    "- the risk-aversion parameter $ \\delta $  \n",
    "\n",
    "\n",
    "A standard way of estimating $ \\mu $ is maximum-likelihood or least\n",
    "squares; that amounts to estimating $ \\mu $ by a sample mean of\n",
    "excess returns and estimating $ \\Sigma $ by a sample covariance\n",
    "matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b44155e9",
   "metadata": {},
   "source": [
    "## Black-Litterman Starting Point\n",
    "\n",
    "When estimates of $ \\mu $ and $ \\Sigma $ from historical\n",
    "sample means and covariances have been combined with **plausible** values\n",
    "of the risk-aversion parameter $ \\delta $ to compute an\n",
    "optimal portfolio from formula [(36.2)](#equation-risky-portfolio), a typical outcome has been\n",
    "$ w $’s with **extreme long and short positions**.\n",
    "\n",
    "A common reaction to these outcomes is that they are so implausible that a portfolio\n",
    "manager cannot recommend them to a customer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eca4baa",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(12)\n",
    "\n",
    "N = 10                                           # Number of assets\n",
    "T = 200                                          # Sample size\n",
    "\n",
    "# random market portfolio (sum is normalized to 1)\n",
    "w_m = np.random.rand(N)\n",
    "w_m = w_m / (w_m.sum())\n",
    "\n",
    "# True risk premia and variance of excess return (constructed\n",
    "# so that the Sharpe ratio is 1)\n",
    "μ = (np.random.randn(N) + 5)  /100      # Mean excess return (risk premium)\n",
    "S = np.random.randn(N, N)        # Random matrix for the covariance matrix\n",
    "V = S @ S.T           # Turn the random matrix into symmetric psd\n",
    "# Make sure that the Sharpe ratio is one\n",
    "Σ = V * (w_m @ μ)**2 / (w_m @ V @ w_m)\n",
    "\n",
    "# Risk aversion of market portfolio holder\n",
    "δ = 1 / np.sqrt(w_m @ Σ @ w_m)\n",
    "\n",
    "# Generate a sample of excess returns\n",
    "excess_return = stat.multivariate_normal(μ, Σ)\n",
    "sample = excess_return.rvs(T)\n",
    "\n",
    "# Estimate μ and Σ\n",
    "μ_est = sample.mean(0).reshape(N, 1)\n",
    "Σ_est = np.cov(sample.T)\n",
    "\n",
    "w = np.linalg.solve(δ * Σ_est, μ_est)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.set_title('Mean-variance portfolio weights recommendation and the market portfolio')\n",
    "ax.plot(np.arange(N)+1, w, 'o', c='k', label='$w$ (mean-variance)')\n",
    "ax.plot(np.arange(N)+1, w_m, 'o', c='r', label='$w_m$ (market portfolio)')\n",
    "ax.vlines(np.arange(N)+1, 0, w, lw=1)\n",
    "ax.vlines(np.arange(N)+1, 0, w_m, lw=1)\n",
    "ax.axhline(0, c='k')\n",
    "ax.axhline(-1, c='k', ls='--')\n",
    "ax.axhline(1, c='k', ls='--')\n",
    "ax.set_xlabel('Assets')\n",
    "ax.xaxis.set_ticks(np.arange(1, N+1, 1))\n",
    "plt.legend(numpoints=1, fontsize=11)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bb2d18a",
   "metadata": {},
   "source": [
    "Black and Litterman’s responded to this situation in the following way:\n",
    "\n",
    "- They continue to accept [(36.2)](#equation-risky-portfolio) as a good model for choosing an optimal\n",
    "  portfolio $ w $.  \n",
    "- They want to continue to allow the customer to express his or her\n",
    "  risk tolerance by setting $ \\delta $.  \n",
    "- Leaving $ \\Sigma $ at its maximum-likelihood value, they push\n",
    "  $ \\mu $ away from its maximum-likelihood value in a way designed to make\n",
    "  portfolio choices that are more plausible in terms of conforming to\n",
    "  what most people actually do.  \n",
    "\n",
    "\n",
    "In particular, given $ \\Sigma $ and a plausible value\n",
    "of $ \\delta $, Black and Litterman reverse engineered a vector\n",
    "$ \\mu_{BL} $ of mean excess returns that makes the $ w $\n",
    "implied by formula [(36.2)](#equation-risky-portfolio) equal the **actual** market portfolio\n",
    "$ w_m $, so that\n",
    "\n",
    "$$\n",
    "w_m = (\\delta \\Sigma)^{-1} \\mu_{BL}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37b6a407",
   "metadata": {},
   "source": [
    "## Details\n",
    "\n",
    "Let’s define\n",
    "\n",
    "$$\n",
    "w_m' \\mu \\equiv ( r_m - r_f)\n",
    "$$\n",
    "\n",
    "as the (scalar) excess return on the market portfolio $ w_m $.\n",
    "\n",
    "Define\n",
    "\n",
    "$$\n",
    "\\sigma^2 = w_m' \\Sigma w_m\n",
    "$$\n",
    "\n",
    "as the variance of the excess return on the market portfolio\n",
    "$ w_m $.\n",
    "\n",
    "Define\n",
    "\n",
    "$$\n",
    "{\\bf SR}_m = \\frac{ r_m - r_f}{\\sigma}\n",
    "$$\n",
    "\n",
    "as the **Sharpe-ratio** on the market portfolio $ w_m $.\n",
    "\n",
    "Let $ \\delta_m $ be the value of the risk aversion parameter that\n",
    "induces an investor to hold the market portfolio in light of the optimal\n",
    "portfolio choice rule [(36.2)](#equation-risky-portfolio).\n",
    "\n",
    "Evidently, portfolio rule [(36.2)](#equation-risky-portfolio) then implies that\n",
    "$ r_m - r_f = \\delta_m \\sigma^2 $ or\n",
    "\n",
    "$$\n",
    "\\delta_m = \\frac{r_m - r_f}{\\sigma^2}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "\\delta_m = \\frac{{\\bf SR}_m}{\\sigma}\n",
    "$$\n",
    "\n",
    "Following the Black-Litterman philosophy, our first step will be to back\n",
    "a value of $ \\delta_m $ from\n",
    "\n",
    "- an estimate of the Sharpe-ratio, and  \n",
    "- our maximum likelihood estimate of $ \\sigma $ drawn from our\n",
    "  estimates or $ w_m $ and $ \\Sigma $  \n",
    "\n",
    "\n",
    "The second key Black-Litterman step is then to use this value of\n",
    "$ \\delta $ together with the maximum likelihood estimate of\n",
    "$ \\Sigma $ to deduce a $ \\mu_{\\bf BL} $ that verifies\n",
    "portfolio rule [(36.2)](#equation-risky-portfolio) at the market portfolio $ w = w_m $\n",
    "\n",
    "$$\n",
    "\\mu_m = \\delta_m \\Sigma w_m\n",
    "$$\n",
    "\n",
    "The starting point of the Black-Litterman portfolio choice model is thus\n",
    "a pair $ (\\delta_m, \\mu_m) $ that tells the customer to hold the\n",
    "market portfolio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4abd55a0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Observed mean excess market return\n",
    "r_m = w_m @ μ_est\n",
    "\n",
    "# Estimated variance of the market portfolio\n",
    "σ_m = w_m @ Σ_est @ w_m\n",
    "\n",
    "# Sharpe-ratio\n",
    "sr_m = r_m / np.sqrt(σ_m)\n",
    "\n",
    "# Risk aversion of market portfolio holder\n",
    "d_m = r_m / σ_m\n",
    "\n",
    "# Derive \"view\" which would induce the market portfolio\n",
    "μ_m = (d_m * Σ_est @ w_m).reshape(N, 1)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.set_title(r'Difference between $\\hat{\\mu}$ (estimate) and $\\mu_{BL}$ (market implied)')\n",
    "ax.plot(np.arange(N)+1, μ_est, 'o', c='k', label='$\\hat{\\mu}$')\n",
    "ax.plot(np.arange(N)+1, μ_m, 'o', c='r', label='$\\mu_{BL}$')\n",
    "ax.vlines(np.arange(N) + 1, μ_m, μ_est, lw=1)\n",
    "ax.axhline(0, c='k', ls='--')\n",
    "ax.set_xlabel('Assets')\n",
    "ax.xaxis.set_ticks(np.arange(1, N+1, 1))\n",
    "plt.legend(numpoints=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acc49438",
   "metadata": {},
   "source": [
    "## Adding Views\n",
    "\n",
    "Black and Litterman start with a baseline customer who asserts that he\n",
    "or she shares the **market’s views**, which means that he or she\n",
    "believes that excess returns are governed by\n",
    "\n",
    "\n",
    "<a id='equation-excess-returns'></a>\n",
    "$$\n",
    "\\vec r - r_f {\\bf 1} \\sim {\\mathcal N}( \\mu_{BL}, \\Sigma) \\tag{36.3}\n",
    "$$\n",
    "\n",
    "Black and Litterman would advise that customer to hold the market\n",
    "portfolio of risky securities.\n",
    "\n",
    "Black and Litterman then imagine a consumer who would like to express a\n",
    "view that differs from the market’s.\n",
    "\n",
    "The consumer wants appropriately to\n",
    "mix his view with the market’s before using [(36.2)](#equation-risky-portfolio) to choose a portfolio.\n",
    "\n",
    "Suppose that the customer’s view is expressed by a hunch that rather\n",
    "than [(36.3)](#equation-excess-returns), excess returns are governed by\n",
    "\n",
    "$$\n",
    "\\vec r - r_f {\\bf 1} \\sim {\\mathcal N}( \\hat \\mu, \\tau \\Sigma)\n",
    "$$\n",
    "\n",
    "where $ \\tau > 0 $ is a scalar parameter that determines how the\n",
    "decision maker wants to mix his view $ \\hat \\mu $ with the market’s\n",
    "view $ \\mu_{\\bf BL} $.\n",
    "\n",
    "Black and Litterman would then use a formula like the following one to\n",
    "mix the views $ \\hat \\mu $ and $ \\mu_{\\bf BL} $\n",
    "\n",
    "\n",
    "<a id='equation-mix-views'></a>\n",
    "$$\n",
    "\\tilde \\mu = (\\Sigma^{-1} + (\\tau \\Sigma)^{-1})^{-1} (\\Sigma^{-1} \\mu_{BL}  + (\\tau \\Sigma)^{-1} \\hat \\mu) \\tag{36.4}\n",
    "$$\n",
    "\n",
    "Black and Litterman would then advise the customer to hold the portfolio\n",
    "associated with these views implied by rule [(36.2)](#equation-risky-portfolio):\n",
    "\n",
    "$$\n",
    "\\tilde w = (\\delta \\Sigma)^{-1} \\tilde \\mu\n",
    "$$\n",
    "\n",
    "This portfolio $ \\tilde w $ will deviate from the\n",
    "portfolio $ w_{BL} $ in amounts that depend on the mixing parameter\n",
    "$ \\tau $.\n",
    "\n",
    "If $ \\hat \\mu $ is the maximum likelihood estimator\n",
    "and $ \\tau $ is chosen heavily to weight this view, then the\n",
    "customer’s portfolio will involve big short-long positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93694ab6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def black_litterman(λ, μ1, μ2, Σ1, Σ2):\n",
    "    \"\"\"\n",
    "    This function calculates the Black-Litterman mixture\n",
    "    mean excess return and covariance matrix\n",
    "    \"\"\"\n",
    "    Σ1_inv = np.linalg.inv(Σ1)\n",
    "    Σ2_inv = np.linalg.inv(Σ2)\n",
    "\n",
    "    μ_tilde = np.linalg.solve(Σ1_inv + λ * Σ2_inv,\n",
    "                              Σ1_inv @ μ1 + λ * Σ2_inv @ μ2)\n",
    "    return μ_tilde\n",
    "\n",
    "τ = 1\n",
    "μ_tilde = black_litterman(1, μ_m, μ_est, Σ_est, τ * Σ_est)\n",
    "\n",
    "# The Black-Litterman recommendation for the portfolio weights\n",
    "w_tilde = np.linalg.solve(δ * Σ_est, μ_tilde)\n",
    "\n",
    "τ_slider = FloatSlider(min=0.05, max=10, step=0.5, value=τ)\n",
    "\n",
    "@interact(τ=τ_slider)\n",
    "def BL_plot(τ):\n",
    "    μ_tilde = black_litterman(1, μ_m, μ_est, Σ_est, τ * Σ_est)\n",
    "    w_tilde = np.linalg.solve(δ * Σ_est, μ_tilde)\n",
    "\n",
    "    fig, ax = plt.subplots(1, 2, figsize=(16, 6))\n",
    "    ax[0].plot(np.arange(N)+1, μ_est, 'o', c='k',\n",
    "               label=r'$\\hat{\\mu}$ (subj view)')\n",
    "    ax[0].plot(np.arange(N)+1, μ_m, 'o', c='r',\n",
    "               label=r'$\\mu_{BL}$ (market)')\n",
    "    ax[0].plot(np.arange(N)+1, μ_tilde, 'o', c='y',\n",
    "               label=r'$\\tilde{\\mu}$ (mixture)')\n",
    "    ax[0].vlines(np.arange(N)+1, μ_m, μ_est, lw=1)\n",
    "    ax[0].axhline(0, c='k', ls='--')\n",
    "    ax[0].set(xlim=(0, N+1), xlabel='Assets',\n",
    "              title=r'Relationship between $\\hat{\\mu}$, $\\mu_{BL}$, and  $ \\tilde{\\mu}$')\n",
    "    ax[0].xaxis.set_ticks(np.arange(1, N+1, 1))\n",
    "    ax[0].legend(numpoints=1)\n",
    "\n",
    "    ax[1].set_title('Black-Litterman portfolio weight recommendation')\n",
    "    ax[1].plot(np.arange(N)+1, w, 'o', c='k', label=r'$w$ (mean-variance)')\n",
    "    ax[1].plot(np.arange(N)+1, w_m, 'o', c='r', label=r'$w_{m}$ (market, BL)')\n",
    "    ax[1].plot(np.arange(N)+1, w_tilde, 'o', c='y',\n",
    "               label=r'$\\tilde{w}$ (mixture)')\n",
    "    ax[1].vlines(np.arange(N)+1, 0, w, lw=1)\n",
    "    ax[1].vlines(np.arange(N)+1, 0, w_m, lw=1)\n",
    "    ax[1].axhline(0, c='k')\n",
    "    ax[1].axhline(-1, c='k', ls='--')\n",
    "    ax[1].axhline(1, c='k', ls='--')\n",
    "    ax[1].set(xlim=(0, N+1), xlabel='Assets',\n",
    "              title='Black-Litterman portfolio weight recommendation')\n",
    "    ax[1].xaxis.set_ticks(np.arange(1, N+1, 1))\n",
    "    ax[1].legend(numpoints=1)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74b4fb99",
   "metadata": {},
   "source": [
    "## Bayesian Interpretation\n",
    "\n",
    "Consider the following Bayesian interpretation of the Black-Litterman\n",
    "recommendation.\n",
    "\n",
    "The prior belief over the mean excess returns is consistent with the\n",
    "market portfolio and is given by\n",
    "\n",
    "$$\n",
    "\\mu \\sim \\mathcal{N}(\\mu_{BL}, \\Sigma)\n",
    "$$\n",
    "\n",
    "Given a particular realization of the mean excess returns\n",
    "$ \\mu $ one observes the average excess returns $ \\hat \\mu $\n",
    "on the market according to the distribution\n",
    "\n",
    "$$\n",
    "\\hat \\mu \\mid \\mu, \\Sigma \\sim \\mathcal{N}(\\mu, \\tau\\Sigma)\n",
    "$$\n",
    "\n",
    "where $ \\tau $ is typically small capturing the idea that the\n",
    "variation in the mean is smaller than the variation of the individual\n",
    "random variable.\n",
    "\n",
    "Given the realized excess returns one should then update the prior over\n",
    "the mean excess returns according to Bayes rule.\n",
    "\n",
    "The corresponding\n",
    "posterior over mean excess returns is normally distributed with mean\n",
    "\n",
    "$$\n",
    "(\\Sigma^{-1} + (\\tau \\Sigma)^{-1})^{-1} (\\Sigma^{-1}\\mu_{BL}   + (\\tau \\Sigma)^{-1} \\hat \\mu)\n",
    "$$\n",
    "\n",
    "The covariance matrix is\n",
    "\n",
    "$$\n",
    "(\\Sigma^{-1} + (\\tau \\Sigma)^{-1})^{-1}\n",
    "$$\n",
    "\n",
    "Hence, the Black-Litterman recommendation is consistent with the Bayes\n",
    "update of the prior over the mean excess returns in light of the\n",
    "realized average excess returns on the market."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f800a0d",
   "metadata": {},
   "source": [
    "## Curve Decolletage\n",
    "\n",
    "Consider two independent “competing” views on the excess market returns\n",
    "\n",
    "$$\n",
    "\\vec r_e  \\sim {\\mathcal N}( \\mu_{BL}, \\Sigma)\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\vec r_e \\sim {\\mathcal N}( \\hat{\\mu}, \\tau\\Sigma)\n",
    "$$\n",
    "\n",
    "A special feature of the multivariate normal random variable\n",
    "$ Z $ is that its density function depends only on the (Euclidiean)\n",
    "length of its realization $ z $.\n",
    "\n",
    "Formally, let the\n",
    "$ k $-dimensional random vector be\n",
    "\n",
    "$$\n",
    "Z\\sim \\mathcal{N}(\\mu, \\Sigma)\n",
    "$$\n",
    "\n",
    "then\n",
    "\n",
    "$$\n",
    "\\bar{Z} \\equiv \\Sigma(Z-\\mu)\\sim \\mathcal{N}(\\mathbf{0}, I)\n",
    "$$\n",
    "\n",
    "and so the points where the density takes the same value can be\n",
    "described by the ellipse\n",
    "\n",
    "\n",
    "<a id='equation-ellipse'></a>\n",
    "$$\n",
    "\\bar z \\cdot \\bar z =  (z - \\mu)'\\Sigma^{-1}(z - \\mu) = \\bar d \\tag{36.5}\n",
    "$$\n",
    "\n",
    "where $ \\bar d\\in\\mathbb{R}_+ $ denotes the (transformation) of a\n",
    "particular density value.\n",
    "\n",
    "The curves defined by equation [(36.5)](#equation-ellipse) can be\n",
    "labeled as iso-likelihood ellipses\n",
    "\n",
    "> **Remark:** More generally there is a class of density functions\n",
    "that possesses this feature, i.e.\n",
    "\n",
    "$$\n",
    "\\exists g: \\mathbb{R}_+ \\mapsto \\mathbb{R}_+ \\ \\ \\text{ and } \\ \\ c \\geq 0,\n",
    "  \\ \\ \\text{s.t.  the density } \\ \\ f \\ \\ \\text{of} \\ \\ Z  \\ \\\n",
    "  \\text{ has the form } \\quad f(z) = c g(z\\cdot z)\n",
    "$$\n",
    "\n",
    "This property is called **spherical symmetry** (see p 81. in Leamer\n",
    "(1978) [[Leamer, 1978](https://python-advanced.quantecon.org/zreferences.html#id67)]).\n",
    "\n",
    "\n",
    "In our specific example, we can use the pair\n",
    "$ (\\bar d_1, \\bar d_2) $ as being two “likelihood” values for which\n",
    "the corresponding iso-likelihood ellipses in the excess return space are\n",
    "given by\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "(\\vec r_e - \\mu_{BL})'\\Sigma^{-1}(\\vec r_e - \\mu_{BL}) &= \\bar d_1 \\\\\n",
    "(\\vec r_e - \\hat \\mu)'\\left(\\tau \\Sigma\\right)^{-1}(\\vec r_e - \\hat \\mu) &= \\bar d_2\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Notice that for particular $ \\bar d_1 $ and $ \\bar d_2 $ values\n",
    "the two ellipses have a tangency point.\n",
    "\n",
    "These tangency points, indexed\n",
    "by the pairs $ (\\bar d_1, \\bar d_2) $, characterize points\n",
    "$ \\vec r_e $ from which there exists no deviation where one can\n",
    "increase the likelihood of one view without decreasing the likelihood of\n",
    "the other view.\n",
    "\n",
    "The pairs $ (\\bar d_1, \\bar d_2) $ for which there\n",
    "is such a point outlines a curve in the excess return space. This curve\n",
    "is reminiscent of the Pareto curve in an Edgeworth-box setting.\n",
    "\n",
    "Dickey (1975) [[Dickey, 1975](https://python-advanced.quantecon.org/zreferences.html#id60)] calls it a *curve decolletage*.\n",
    "\n",
    "Leamer (1978) [[Leamer, 1978](https://python-advanced.quantecon.org/zreferences.html#id67)] calls it an *information contract curve* and\n",
    "describes it by the following program: maximize the likelihood of one\n",
    "view, say the Black-Litterman recommendation while keeping the\n",
    "likelihood of the other view at least at a prespecified constant\n",
    "$ \\bar d_2 $\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    " \\bar d_1(\\bar d_2) &\\equiv \\max_{\\vec r_e} \\ \\ (\\vec r_e - \\mu_{BL})'\\Sigma^{-1}(\\vec r_e - \\mu_{BL}) \\\\\n",
    "\\text{subject to }  \\quad  &(\\vec r_e - \\hat\\mu)'(\\tau\\Sigma)^{-1}(\\vec r_e - \\hat \\mu) \\geq \\bar d_2\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Denoting the multiplier on the constraint by $ \\lambda $, the\n",
    "first-order condition is\n",
    "\n",
    "$$\n",
    "2(\\vec r_e - \\mu_{BL} )'\\Sigma^{-1} + \\lambda 2(\\vec r_e - \\hat\\mu)'(\\tau\\Sigma)^{-1} = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "which defines the *information contract curve* between\n",
    "$ \\mu_{BL} $ and $ \\hat \\mu $\n",
    "\n",
    "\n",
    "<a id='equation-info-curve'></a>\n",
    "$$\n",
    "\\vec r_e = (\\Sigma^{-1} + \\lambda (\\tau \\Sigma)^{-1})^{-1} (\\Sigma^{-1} \\mu_{BL} + \\lambda (\\tau \\Sigma)^{-1}\\hat \\mu ) \\tag{36.6}\n",
    "$$\n",
    "\n",
    "Note that if $ \\lambda = 1 $, [(36.6)](#equation-info-curve) is equivalent with [(36.4)](#equation-mix-views) and it\n",
    "identifies one point on the information contract curve.\n",
    "\n",
    "Furthermore, because $ \\lambda $ is a function of the minimum likelihood\n",
    "$ \\bar d_2 $ on the RHS of the constraint, by varying\n",
    "$ \\bar d_2 $ (or $ \\lambda $ ), we can trace out the whole curve\n",
    "as the figure below illustrates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "358e9567",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(1987102)\n",
    "\n",
    "N = 2                                           # Number of assets\n",
    "T = 200                                         # Sample size\n",
    "τ = 0.8\n",
    "\n",
    "# Random market portfolio (sum is normalized to 1)\n",
    "w_m = np.random.rand(N)\n",
    "w_m = w_m / (w_m.sum())\n",
    "\n",
    "μ = (np.random.randn(N) + 5) / 100\n",
    "S = np.random.randn(N, N)\n",
    "V = S @ S.T\n",
    "Σ = V * (w_m @ μ)**2 / (w_m @ V @ w_m)\n",
    "\n",
    "excess_return = stat.multivariate_normal(μ, Σ)\n",
    "sample = excess_return.rvs(T)\n",
    "\n",
    "μ_est = sample.mean(0).reshape(N, 1)\n",
    "Σ_est = np.cov(sample.T)\n",
    "\n",
    "σ_m = w_m @ Σ_est @ w_m\n",
    "d_m = (w_m @ μ_est) / σ_m\n",
    "μ_m = (d_m * Σ_est @ w_m).reshape(N, 1)\n",
    "\n",
    "N_r1, N_r2 = 100, 100\n",
    "r1 = np.linspace(-0.04, .1, N_r1)\n",
    "r2 = np.linspace(-0.02, .15, N_r2)\n",
    "\n",
    "λ_grid = np.linspace(.001, 20, 100)\n",
    "curve = np.asarray([black_litterman(λ, μ_m, μ_est, Σ_est,\n",
    "                                    τ * Σ_est).flatten() for λ in λ_grid])\n",
    "\n",
    "λ_slider = FloatSlider(min=.1, max=7, step=.5, value=1)\n",
    "\n",
    "@interact(λ=λ_slider)\n",
    "def decolletage(λ):\n",
    "    dist_r_BL = stat.multivariate_normal(μ_m.squeeze(), Σ_est)\n",
    "    dist_r_hat = stat.multivariate_normal(μ_est.squeeze(), τ * Σ_est)\n",
    "\n",
    "    X, Y = np.meshgrid(r1, r2)\n",
    "    Z_BL = np.zeros((N_r1, N_r2))\n",
    "    Z_hat = np.zeros((N_r1, N_r2))\n",
    "\n",
    "    for i in range(N_r1):\n",
    "        for j in range(N_r2):\n",
    "            Z_BL[i, j] = dist_r_BL.pdf(np.hstack([X[i, j], Y[i, j]]))\n",
    "            Z_hat[i, j] = dist_r_hat.pdf(np.hstack([X[i, j], Y[i, j]]))\n",
    "\n",
    "    μ_tilde = black_litterman(λ, μ_m, μ_est, Σ_est, τ * Σ_est).flatten()\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(10, 6))\n",
    "    ax.contourf(X, Y, Z_hat, cmap='viridis', alpha =.4)\n",
    "    ax.contourf(X, Y, Z_BL, cmap='viridis', alpha =.4)\n",
    "    ax.contour(X, Y, Z_BL, [dist_r_BL.pdf(μ_tilde)], cmap='viridis', alpha=.9)\n",
    "    ax.contour(X, Y, Z_hat, [dist_r_hat.pdf(μ_tilde)], cmap='viridis', alpha=.9)\n",
    "    ax.scatter(μ_est[0], μ_est[1])\n",
    "    ax.scatter(μ_m[0], μ_m[1])\n",
    "    ax.scatter(μ_tilde[0], μ_tilde[1], c='k', s=20*3)\n",
    "\n",
    "    ax.plot(curve[:, 0], curve[:, 1], c='k')\n",
    "    ax.axhline(0, c='k', alpha=.8)\n",
    "    ax.axvline(0, c='k', alpha=.8)\n",
    "    ax.set_xlabel(r'Excess return on the first asset, $r_{e, 1}$')\n",
    "    ax.set_ylabel(r'Excess return on the second asset, $r_{e, 2}$')\n",
    "    ax.text(μ_est[0] + 0.003, μ_est[1], r'$\\hat{\\mu}$')\n",
    "    ax.text(μ_m[0] + 0.003, μ_m[1] + 0.005, r'$\\mu_{BL}$')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b612c5ba",
   "metadata": {},
   "source": [
    "Note that the line that connects the two points\n",
    "$ \\hat \\mu $ and $ \\mu_{BL} $ is linear, which comes from the\n",
    "fact that the covariance matrices of the two competing distributions\n",
    "(views) are proportional to each other.\n",
    "\n",
    "To illustrate the fact that this is not necessarily the case, consider\n",
    "another example using the same parameter values, except that the “second\n",
    "view” constituting the constraint has covariance matrix\n",
    "$ \\tau I $ instead of $ \\tau \\Sigma $.\n",
    "\n",
    "This leads to the\n",
    "following figure, on which the curve connecting $ \\hat \\mu $\n",
    "and $ \\mu_{BL} $ are bending"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "474ef5a6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "λ_grid = np.linspace(.001, 20000, 1000)\n",
    "curve = np.asarray([black_litterman(λ, μ_m, μ_est, Σ_est,\n",
    "                                    τ * np.eye(N)).flatten() for λ in λ_grid])\n",
    "\n",
    "λ_slider = FloatSlider(min=5, max=1500, step=100, value=200)\n",
    "\n",
    "@interact(λ=λ_slider)\n",
    "def decolletage(λ):\n",
    "    dist_r_BL = stat.multivariate_normal(μ_m.squeeze(), Σ_est)\n",
    "    dist_r_hat = stat.multivariate_normal(μ_est.squeeze(), τ * np.eye(N))\n",
    "\n",
    "    X, Y = np.meshgrid(r1, r2)\n",
    "    Z_BL = np.zeros((N_r1, N_r2))\n",
    "    Z_hat = np.zeros((N_r1, N_r2))\n",
    "\n",
    "    for i in range(N_r1):\n",
    "        for j in range(N_r2):\n",
    "            Z_BL[i, j] = dist_r_BL.pdf(np.hstack([X[i, j], Y[i, j]]))\n",
    "            Z_hat[i, j] = dist_r_hat.pdf(np.hstack([X[i, j], Y[i, j]]))\n",
    "\n",
    "    μ_tilde = black_litterman(λ, μ_m, μ_est, Σ_est, τ * np.eye(N)).flatten()\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(10, 6))\n",
    "    ax.contourf(X, Y, Z_hat, cmap='viridis', alpha=.4)\n",
    "    ax.contourf(X, Y, Z_BL, cmap='viridis', alpha=.4)\n",
    "    ax.contour(X, Y, Z_BL, [dist_r_BL.pdf(μ_tilde)], cmap='viridis', alpha=.9)\n",
    "    ax.contour(X, Y, Z_hat, [dist_r_hat.pdf(μ_tilde)], cmap='viridis', alpha=.9)\n",
    "    ax.scatter(μ_est[0], μ_est[1])\n",
    "    ax.scatter(μ_m[0], μ_m[1])\n",
    "\n",
    "    ax.scatter(μ_tilde[0], μ_tilde[1], c='k', s=20*3)\n",
    "\n",
    "    ax.plot(curve[:, 0], curve[:, 1], c='k')\n",
    "    ax.axhline(0, c='k', alpha=.8)\n",
    "    ax.axvline(0, c='k', alpha=.8)\n",
    "    ax.set_xlabel(r'Excess return on the first asset, $r_{e, 1}$')\n",
    "    ax.set_ylabel(r'Excess return on the second asset, $r_{e, 2}$')\n",
    "    ax.text(μ_est[0] + 0.003, μ_est[1], r'$\\hat{\\mu}$')\n",
    "    ax.text(μ_m[0] + 0.003, μ_m[1] + 0.005, r'$\\mu_{BL}$')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf558668",
   "metadata": {},
   "source": [
    "## Black-Litterman Recommendation as Regularization\n",
    "\n",
    "First, consider the OLS regression\n",
    "\n",
    "$$\n",
    "\\min_{\\beta} \\Vert X\\beta - y \\Vert^2\n",
    "$$\n",
    "\n",
    "which yields the solution\n",
    "\n",
    "$$\n",
    "\\hat{\\beta}_{OLS} = (X'X)^{-1}X'y\n",
    "$$\n",
    "\n",
    "A common performance measure of estimators is the *mean squared error\n",
    "(MSE)*.\n",
    "\n",
    "An estimator is “good” if its MSE is relatively small. Suppose\n",
    "that $ \\beta_0 $ is the “true” value of the coefficient, then the MSE\n",
    "of the OLS estimator is\n",
    "\n",
    "$$\n",
    "\\text{mse}(\\hat \\beta_{OLS}, \\beta_0) := \\mathbb E \\Vert \\hat \\beta_{OLS} - \\beta_0\\Vert^2 =\n",
    "\\underbrace{\\mathbb E \\Vert \\hat \\beta_{OLS} - \\mathbb E\n",
    "\\beta_{OLS}\\Vert^2}_{\\text{variance}} +\n",
    "\\underbrace{\\Vert \\mathbb E \\hat\\beta_{OLS} - \\beta_0\\Vert^2}_{\\text{bias}}\n",
    "$$\n",
    "\n",
    "From this decomposition, one can see that in order for the MSE to be\n",
    "small, both the bias and the variance terms must be small.\n",
    "\n",
    "For example,\n",
    "consider the case when $ X $ is a $ T $-vector of ones (where\n",
    "$ T $ is the sample size), so $ \\hat\\beta_{OLS} $ is simply the\n",
    "sample average, while $ \\beta_0\\in \\mathbb{R} $ is defined by the\n",
    "true mean of $ y $.\n",
    "\n",
    "In this example the MSE is\n",
    "\n",
    "$$\n",
    "\\text{mse}(\\hat \\beta_{OLS}, \\beta_0) = \\underbrace{\\frac{1}{T^2}\n",
    "\\mathbb E \\left(\\sum_{t=1}^{T} (y_{t}- \\beta_0)\\right)^2 }_{\\text{variance}} +\n",
    "\\underbrace{0}_{\\text{bias}}\n",
    "$$\n",
    "\n",
    "However, because there is a trade-off between the estimator’s bias and\n",
    "variance, there are cases when by permitting a small bias we can\n",
    "substantially reduce the variance so overall the MSE gets smaller.\n",
    "\n",
    "A typical scenario when this proves to be useful is when the number of\n",
    "coefficients to be estimated is large relative to the sample size.\n",
    "\n",
    "In these cases, one approach to handle the bias-variance trade-off is the\n",
    "so called *Tikhonov regularization*.\n",
    "\n",
    "A general form with regularization matrix $ \\Gamma $ can be written as\n",
    "\n",
    "$$\n",
    "\\min_{\\beta} \\Big\\{ \\Vert X\\beta - y \\Vert^2 + \\Vert \\Gamma (\\beta - \\tilde \\beta) \\Vert^2 \\Big\\}\n",
    "$$\n",
    "\n",
    "which yields the solution\n",
    "\n",
    "$$\n",
    "\\hat{\\beta}_{Reg} = (X'X + \\Gamma'\\Gamma)^{-1}(X'y + \\Gamma'\\Gamma\\tilde \\beta)\n",
    "$$\n",
    "\n",
    "Substituting the value of $ \\hat{\\beta}_{OLS} $ yields\n",
    "\n",
    "$$\n",
    "\\hat{\\beta}_{Reg} = (X'X + \\Gamma'\\Gamma)^{-1}(X'X\\hat{\\beta}_{OLS} + \\Gamma'\\Gamma\\tilde \\beta)\n",
    "$$\n",
    "\n",
    "Often, the regularization matrix takes the form\n",
    "$ \\Gamma = \\lambda I $ with $ \\lambda>0 $\n",
    "and $ \\tilde \\beta = \\mathbf{0} $.\n",
    "\n",
    "Then the Tikhonov regularization is equivalent to what is called *ridge regression* in statistics.\n",
    "\n",
    "To illustrate how this estimator addresses the bias-variance trade-off,\n",
    "we compute the MSE of the ridge estimator\n",
    "\n",
    "$$\n",
    "\\text{mse}(\\hat \\beta_{\\text{ridge}}, \\beta_0) = \\underbrace{\\frac{1}{(T+\\lambda)^2}\n",
    "\\mathbb E \\left(\\sum_{t=1}^{T} (y_{t}- \\beta_0)\\right)^2 }_{\\text{variance}} +\n",
    "\\underbrace{\\left(\\frac{\\lambda}{T+\\lambda}\\right)^2 \\beta_0^2}_{\\text{bias}}\n",
    "$$\n",
    "\n",
    "The ridge regression shrinks the coefficients of the estimated vector\n",
    "towards zero relative to the OLS estimates thus reducing the variance\n",
    "term at the cost of introducing a “small” bias.\n",
    "\n",
    "However, there is nothing special about the zero vector.\n",
    "\n",
    "When $ \\tilde \\beta \\neq \\mathbf{0} $ shrinkage occurs in the direction\n",
    "of $ \\tilde \\beta $.\n",
    "\n",
    "Now, we can give a regularization interpretation of the Black-Litterman\n",
    "portfolio recommendation.\n",
    "\n",
    "To this end, first  simplify the equation [(36.4)](#equation-mix-views) that characterizes the Black-Litterman recommendation\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\tilde \\mu &= (\\Sigma^{-1} + (\\tau \\Sigma)^{-1})^{-1} (\\Sigma^{-1}\\mu_{BL}  + (\\tau \\Sigma)^{-1}\\hat \\mu) \\\\\n",
    "&= (1 + \\tau^{-1})^{-1}\\Sigma \\Sigma^{-1} (\\mu_{BL}  + \\tau ^{-1}\\hat \\mu) \\\\\n",
    "&= (1 + \\tau^{-1})^{-1} ( \\mu_{BL}  + \\tau ^{-1}\\hat \\mu)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "In our case, $ \\hat \\mu $ is the estimated mean excess returns of\n",
    "securities. This could be written as a vector autoregression where\n",
    "\n",
    "- $ y $ is the stacked vector of observed excess returns of size\n",
    "  $ (N T\\times 1) $ – $ N $ securities and $ T $\n",
    "  observations.  \n",
    "- $ X = \\sqrt{T^{-1}}(I_{N} \\otimes \\iota_T) $ where $ I_N $\n",
    "  is the identity matrix and $ \\iota_T $ is a column vector of\n",
    "  ones.  \n",
    "\n",
    "\n",
    "Correspondingly, the OLS regression of $ y $ on $ X $ would\n",
    "yield the mean excess returns as coefficients.\n",
    "\n",
    "With $ \\Gamma = \\sqrt{\\tau T^{-1}}(I_{N} \\otimes \\iota_T) $ we can\n",
    "write the regularized version of the mean excess return estimation\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\hat{\\beta}_{Reg} &= (X'X + \\Gamma'\\Gamma)^{-1}(X'X\\hat{\\beta}_{OLS} + \\Gamma'\\Gamma\\tilde \\beta) \\\\\n",
    "&= (1 + \\tau)^{-1}X'X (X'X)^{-1} (\\hat \\beta_{OLS}  + \\tau \\tilde \\beta) \\\\\n",
    "&= (1 + \\tau)^{-1} (\\hat \\beta_{OLS}  + \\tau \\tilde \\beta) \\\\\n",
    "&= (1 + \\tau^{-1})^{-1} ( \\tau^{-1}\\hat \\beta_{OLS}  +  \\tilde \\beta)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Given that\n",
    "$ \\hat \\beta_{OLS} = \\hat \\mu $ and $ \\tilde \\beta = \\mu_{BL} $\n",
    "in the Black-Litterman model, we have the following interpretation of the\n",
    "model’s recommendation.\n",
    "\n",
    "The estimated (personal) view of the mean excess returns,\n",
    "$ \\hat{\\mu} $ that would lead to extreme short-long positions are\n",
    "“shrunk” towards the conservative market view, $ \\mu_{BL} $, that\n",
    "leads to the more conservative market portfolio.\n",
    "\n",
    "So the Black-Litterman procedure results in a recommendation that is a\n",
    "compromise between the conservative market portfolio and the more\n",
    "extreme portfolio that is implied by estimated “personal” views."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95bebd5d",
   "metadata": {},
   "source": [
    "## A Robust Control Operator\n",
    "\n",
    "The Black-Litterman approach is partly inspired by the econometric\n",
    "insight that it is easier to estimate covariances of excess returns than\n",
    "the means.\n",
    "\n",
    "That is what gave Black and Litterman license to adjust\n",
    "investors’ perception of mean excess returns while not tampering with\n",
    "the covariance matrix of excess returns.\n",
    "\n",
    "The robust control theory is another approach that also hinges on\n",
    "adjusting mean excess returns but not covariances.\n",
    "\n",
    "Associated with a robust control problem is what Hansen and Sargent [[Hansen and Sargent, 2001](https://python-advanced.quantecon.org/zreferences.html#id59)], [[Hansen and Sargent, 2008](https://python-advanced.quantecon.org/zreferences.html#id157)] call\n",
    "a $ {\\sf T} $ operator.\n",
    "\n",
    "Let’s define the $ {\\sf T} $ operator as it applies to the problem\n",
    "at hand.\n",
    "\n",
    "Let $ x $ be an $ n \\times 1 $ Gaussian random vector with mean\n",
    "vector $ \\mu $ and covariance matrix $ \\Sigma = C C' $. This\n",
    "means that $ x $ can be represented as\n",
    "\n",
    "$$\n",
    "x = \\mu + C \\epsilon\n",
    "$$\n",
    "\n",
    "where $ \\epsilon \\sim {\\mathcal N}(0,I) $.\n",
    "\n",
    "Let $ \\phi(\\epsilon) $ denote the associated standardized Gaussian\n",
    "density.\n",
    "\n",
    "Let $ m(\\epsilon,\\mu) $ be a **likelihood ratio**, meaning that it\n",
    "satisfies\n",
    "\n",
    "- $ m(\\epsilon, \\mu) > 0 $  \n",
    "- $ \\int m(\\epsilon,\\mu) \\phi(\\epsilon) d \\epsilon =1 $  \n",
    "\n",
    "\n",
    "That is, $ m(\\epsilon, \\mu) $ is a non-negative random variable with\n",
    "mean 1.\n",
    "\n",
    "Multiplying $ \\phi(\\epsilon) $\n",
    "by the likelihood ratio $ m(\\epsilon, \\mu) $ produces a distorted distribution for\n",
    "$ \\epsilon $, namely\n",
    "\n",
    "$$\n",
    "\\tilde \\phi(\\epsilon) = m(\\epsilon,\\mu) \\phi(\\epsilon)\n",
    "$$\n",
    "\n",
    "The next concept that we need is the **entropy** of the distorted\n",
    "distribution $ \\tilde \\phi $ with respect to $ \\phi $.\n",
    "\n",
    "**Entropy** is defined as\n",
    "\n",
    "$$\n",
    "{\\rm ent} = \\int \\log m(\\epsilon,\\mu) m(\\epsilon,\\mu) \\phi(\\epsilon) d \\epsilon\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "{\\rm ent} = \\int \\log m(\\epsilon,\\mu) \\tilde \\phi(\\epsilon) d \\epsilon\n",
    "$$\n",
    "\n",
    "That is, relative entropy is the expected value of the likelihood ratio\n",
    "$ m $ where the expectation is taken with respect to the twisted\n",
    "density $ \\tilde \\phi $.\n",
    "\n",
    "Relative entropy is non-negative. It is a measure of the discrepancy\n",
    "between two probability distributions.\n",
    "\n",
    "As such, it plays an important\n",
    "role in governing the behavior of statistical tests designed to\n",
    "discriminate one probability distribution from another.\n",
    "\n",
    "We are ready to define the $ {\\sf T} $ operator.\n",
    "\n",
    "Let $ V(x) $ be a value function.\n",
    "\n",
    "Define\n",
    "\n",
    "$$\n",
    "\\begin{aligned} {\\sf T}\\left(V(x)\\right) & = \\min_{m(\\epsilon,\\mu)} \\int m(\\epsilon,\\mu)[V(\\mu + C \\epsilon) + \\theta \\log m(\\epsilon,\\mu) ] \\phi(\\epsilon) d \\epsilon \\cr\n",
    "                        & = - \\log \\theta \\int \\exp \\left( \\frac{- V(\\mu + C \\epsilon)}{\\theta} \\right) \\phi(\\epsilon) d \\epsilon \\end{aligned}\n",
    "$$\n",
    "\n",
    "This asserts that $ {\\sf T} $ is an indirect utility function for a\n",
    "minimization problem in which an **adversary** chooses a distorted\n",
    "probability distribution $ \\tilde \\phi $ to lower expected utility,\n",
    "subject to a penalty term that gets bigger the larger is relative\n",
    "entropy.\n",
    "\n",
    "Here the penalty parameter\n",
    "\n",
    "$$\n",
    "\\theta \\in [\\underline \\theta, +\\infty]\n",
    "$$\n",
    "\n",
    "is a robustness parameter when it is $ +\\infty $, there is no scope for the minimizing agent to distort the distribution,\n",
    "so no robustness to alternative distributions is acquired.\n",
    "\n",
    "As $ \\theta $ is lowered, more robustness is achieved.\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">The $ {\\sf T} $ operator is sometimes called a\n",
    "*risk-sensitivity* operator.\n",
    "\n",
    "We shall apply $ {\\sf T} $ to the special case of a linear value\n",
    "function $ w'(\\vec r - r_f 1) $\n",
    "where $ \\vec r - r_f 1 \\sim {\\mathcal N}(\\mu,\\Sigma) $ or\n",
    "$ \\vec r - r_f {\\bf 1} = \\mu + C \\epsilon $ and\n",
    "$ \\epsilon \\sim {\\mathcal N}(0,I) $.\n",
    "\n",
    "The associated worst-case distribution of $ \\epsilon $ is Gaussian\n",
    "with mean $ v =-\\theta^{-1} C' w $ and covariance matrix $ I $\n",
    "\n",
    "(When the value function is affine, the worst-case distribution distorts\n",
    "the mean vector of $ \\epsilon $ but not the covariance matrix\n",
    "of $ \\epsilon $).\n",
    "\n",
    "For utility function argument $ w'(\\vec r - r_f 1) $\n",
    "\n",
    "$$\n",
    "{\\sf T} ( \\vec r - r_f {\\bf 1}) = w' \\mu + \\zeta - \\frac{1}{2 \\theta} w' \\Sigma w\n",
    "$$\n",
    "\n",
    "and entropy is\n",
    "\n",
    "$$\n",
    "\\frac{v'v}{2} = \\frac{1}{2\\theta^2}  w' C C' w\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "000da59f",
   "metadata": {},
   "source": [
    "## A Robust Mean-Variance Portfolio Model\n",
    "\n",
    "According to criterion [(36.1)](#equation-choice-problem), the mean-variance portfolio choice problem\n",
    "chooses $ w $ to maximize\n",
    "\n",
    "$$\n",
    "E [w ( \\vec r - r_f {\\bf 1})]] - {\\rm var} [ w ( \\vec r - r_f {\\bf 1}) ]\n",
    "$$\n",
    "\n",
    "which equals\n",
    "\n",
    "$$\n",
    "w'\\mu - \\frac{\\delta}{2} w' \\Sigma w\n",
    "$$\n",
    "\n",
    "A robust decision maker can be modeled as replacing the mean return\n",
    "$ E [w ( \\vec r - r_f {\\bf 1})] $ with the risk-sensitive criterion\n",
    "\n",
    "$$\n",
    "{\\sf T} [w ( \\vec r - r_f {\\bf 1})] = w' \\mu - \\frac{1}{2 \\theta} w' \\Sigma w\n",
    "$$\n",
    "\n",
    "that comes from replacing the mean $ \\mu $ of $ \\vec r - r\\_f {\\bf 1} $ with the worst-case mean\n",
    "\n",
    "$$\n",
    "\\mu - \\theta^{-1} \\Sigma w\n",
    "$$\n",
    "\n",
    "Notice how the worst-case mean vector depends on the portfolio\n",
    "$ w $.\n",
    "\n",
    "The operator $ {\\sf T} $ is the indirect utility function that\n",
    "emerges from solving a problem in which an agent who chooses\n",
    "probabilities does so in order to minimize the expected utility of a\n",
    "maximizing agent (in our case, the maximizing agent chooses portfolio\n",
    "weights $ w $).\n",
    "\n",
    "The robust version of the mean-variance portfolio choice problem is then\n",
    "to choose a portfolio $ w $ that maximizes\n",
    "\n",
    "$$\n",
    "{\\sf T} [w ( \\vec r - r_f {\\bf 1})] - \\frac{\\delta}{2} w' \\Sigma w\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-robust-mean-variance'></a>\n",
    "$$\n",
    "w' (\\mu - \\theta^{-1} \\Sigma w ) - \\frac{\\delta}{2} w' \\Sigma w \\tag{36.7}\n",
    "$$\n",
    "\n",
    "The minimizer of [(36.7)](#equation-robust-mean-variance) is\n",
    "\n",
    "$$\n",
    "w_{\\rm rob} = \\frac{1}{\\delta + \\gamma } \\Sigma^{-1} \\mu\n",
    "$$\n",
    "\n",
    "where $ \\gamma \\equiv \\theta^{-1} $ is sometimes called the\n",
    "risk-sensitivity parameter.\n",
    "\n",
    "An increase in the risk-sensitivity parameter $ \\gamma $ shrinks the\n",
    "portfolio weights toward zero in the same way that an increase in risk\n",
    "aversion does."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb1044d9",
   "metadata": {},
   "source": [
    "## Appendix\n",
    "\n",
    "We want to illustrate the “folk theorem” that with high or moderate\n",
    "frequency data, it is more difficult to estimate means than variances.\n",
    "\n",
    "In order to operationalize this statement, we take two analog\n",
    "estimators:\n",
    "\n",
    "- sample average: $ \\bar X_N = \\frac{1}{N}\\sum_{i=1}^{N} X_i $  \n",
    "- sample variance:\n",
    "  $ S_N = \\frac{1}{N-1}\\sum_{t=1}^{N} (X_i - \\bar X_N)^2 $  \n",
    "\n",
    "\n",
    "to estimate the unconditional mean and unconditional variance of the\n",
    "random variable $ X $, respectively.\n",
    "\n",
    "To measure the “difficulty of estimation”, we use *mean squared error*\n",
    "(MSE), that is the average squared difference between the estimator and\n",
    "the true value.\n",
    "\n",
    "Assuming that the process $ \\{X_i\\} $is ergodic,\n",
    "both analog estimators are known to converge to their true values as the\n",
    "sample size $ N $ goes to infinity.\n",
    "\n",
    "More precisely\n",
    "for all $ \\varepsilon > 0 $\n",
    "\n",
    "$$\n",
    "\\lim_{N\\to \\infty} \\ \\ P\\left\\{ \\left |\\bar X_N - \\mathbb E X \\right| > \\varepsilon \\right\\} = 0 \\quad \\quad\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\lim_{N\\to \\infty} \\ \\ P \\left\\{ \\left| S_N - \\mathbb V X \\right| > \\varepsilon \\right\\} = 0\n",
    "$$\n",
    "\n",
    "A necessary condition for these convergence results is that the\n",
    "associated MSEs vanish as $ N $ goes to infinity, or in other words,\n",
    "\n",
    "$$\n",
    "\\text{MSE}(\\bar X_N, \\mathbb E X) = o(1) \\quad \\quad  \\text{and} \\quad \\quad \\text{MSE}(S_N, \\mathbb V X) = o(1)\n",
    "$$\n",
    "\n",
    "Even if the MSEs converge to zero, the associated rates might be\n",
    "different. Looking at the limit of the *relative MSE* (as the sample\n",
    "size grows to infinity)\n",
    "\n",
    "$$\n",
    "\\frac{\\text{MSE}(S_N, \\mathbb V X)}{\\text{MSE}(\\bar X_N, \\mathbb E X)} = \\frac{o(1)}{o(1)} \\underset{N \\to \\infty}{\\to} B\n",
    "$$\n",
    "\n",
    "can inform us about the relative (asymptotic) rates.\n",
    "\n",
    "We will show that in general, with dependent data, the limit\n",
    "$ B $ depends on the sampling frequency.\n",
    "\n",
    "In particular, we find\n",
    "that the rate of convergence of the variance estimator is less sensitive\n",
    "to increased sampling frequency than the rate of convergence of the mean\n",
    "estimator.\n",
    "\n",
    "Hence, we can expect the relative asymptotic\n",
    "rate, $ B $, to get smaller with higher frequency data,\n",
    "illustrating that “it is more difficult to estimate means than\n",
    "variances”.\n",
    "\n",
    "That is, we need significantly more data to obtain a given\n",
    "precision of the mean estimate than for our variance estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ffc6947",
   "metadata": {},
   "source": [
    "## Special Case – IID Sample\n",
    "\n",
    "We start our analysis with the benchmark case of IID data.\n",
    "\n",
    "Consider a\n",
    "sample of size $ N $ generated by the following IID process,\n",
    "\n",
    "$$\n",
    "X_i \\sim \\mathcal{N}(\\mu, \\sigma^2)\n",
    "$$\n",
    "\n",
    "Taking $ \\bar X_N $ to estimate the mean, the MSE is\n",
    "\n",
    "$$\n",
    "\\text{MSE}(\\bar X_N, \\mu) = \\frac{\\sigma^2}{N}\n",
    "$$\n",
    "\n",
    "Taking $ S_N $ to estimate the variance, the MSE is\n",
    "\n",
    "$$\n",
    "\\text{MSE}(S_N, \\sigma^2) = \\frac{2\\sigma^4}{N-1}\n",
    "$$\n",
    "\n",
    "Both estimators are unbiased and hence the MSEs reflect the\n",
    "corresponding variances of the estimators.\n",
    "\n",
    "Furthermore, both MSEs are\n",
    "$ o(1) $ with a (multiplicative) factor of difference in their rates\n",
    "of convergence:\n",
    "\n",
    "$$\n",
    "\\frac{\\text{MSE}(S_N, \\sigma^2)}{\\text{MSE}(\\bar X_N, \\mu)} = \\frac{N2\\sigma^2}{N-1} \\quad \\underset{N \\to \\infty}{\\to} \\quad 2\\sigma^2\n",
    "$$\n",
    "\n",
    "We are interested in how this (asymptotic) relative rate of convergence\n",
    "changes as increasing sampling frequency puts dependence into the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f55f8b7f",
   "metadata": {},
   "source": [
    "## Dependence and Sampling Frequency\n",
    "\n",
    "To investigate how sampling frequency affects relative rates of\n",
    "convergence, we assume that the data are generated by a mean-reverting\n",
    "continuous time process of the form\n",
    "\n",
    "$$\n",
    "dX_t = -\\kappa (X_t -\\mu)dt + \\sigma dW_t\\quad\\quad\n",
    "$$\n",
    "\n",
    "where $ \\mu $is the unconditional mean, $ \\kappa > 0 $ is a\n",
    "persistence parameter, and $ \\{W_t\\} $ is a standardized Brownian\n",
    "motion.\n",
    "\n",
    "Observations arising from this system in particular discrete periods\n",
    "$ \\mathcal T(h) \\equiv \\{nh : n \\in \\mathbb Z \\} $with$ h>0 $\n",
    "can be described by the following process\n",
    "\n",
    "$$\n",
    "X_{t+1} = (1 - \\exp(-\\kappa h))\\mu + \\exp(-\\kappa h)X_t + \\epsilon_{t, h}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\epsilon_{t, h} \\sim \\mathcal{N}(0, \\Sigma_h) \\quad \\text{with}\\quad \\Sigma_h = \\frac{\\sigma^2(1-\\exp(-2\\kappa h))}{2\\kappa}\n",
    "$$\n",
    "\n",
    "We call $ h $ the *frequency* parameter, whereas $ n $\n",
    "represents the number of *lags* between observations.\n",
    "\n",
    "Hence, the effective distance between two observations $ X_t $ and\n",
    "$ X_{t+n} $ in the discrete time notation is equal\n",
    "to $ h\\cdot n $ in terms of the underlying continuous time process.\n",
    "\n",
    "Straightforward calculations show that the autocorrelation function for\n",
    "the stochastic process $ \\{X_{t}\\}_{t\\in \\mathcal T(h)} $ is\n",
    "\n",
    "$$\n",
    "\\Gamma_h(n) \\equiv \\text{corr}(X_{t + h n}, X_t) = \\exp(-\\kappa h n)\n",
    "$$\n",
    "\n",
    "and the auto-covariance function is\n",
    "\n",
    "$$\n",
    "\\gamma_h(n) \\equiv \\text{cov}(X_{t + h n}, X_t) = \\frac{\\exp(-\\kappa h n)\\sigma^2}{2\\kappa} .\n",
    "$$\n",
    "\n",
    "It follows that if $ n=0 $, the unconditional variance is given\n",
    "by $ \\gamma_h(0) = \\frac{\\sigma^2}{2\\kappa} $ irrespective of the\n",
    "sampling frequency.\n",
    "\n",
    "The following figure illustrates how the dependence between the\n",
    "observations is related to the sampling frequency\n",
    "\n",
    "- For any given $ h $, the autocorrelation converges to zero as we increase the distance – $ n $– between the observations. This represents the “weak dependence” of the $ X $ process.  \n",
    "- 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.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37151cf3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "μ = .0\n",
    "κ = .1\n",
    "σ = .5\n",
    "var_uncond = σ**2 / (2 * κ)\n",
    "\n",
    "n_grid = np.linspace(0, 40, 100)\n",
    "autocorr_h1 = np.exp(-κ * n_grid * 1)\n",
    "autocorr_h2 = np.exp(-κ * n_grid * 2)\n",
    "autocorr_h5 = np.exp(-κ * n_grid * 5)\n",
    "autocorr_h1000 = np.exp(-κ * n_grid * 1e8)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.plot(n_grid, autocorr_h1, label=r'$h=1$', c='darkblue', lw=2)\n",
    "ax.plot(n_grid, autocorr_h2, label=r'$h=2$', c='darkred', lw=2)\n",
    "ax.plot(n_grid, autocorr_h5, label=r'$h=5$', c='orange', lw=2)\n",
    "ax.plot(n_grid, autocorr_h1000, label=r'\"$h=\\infty$\"', c='darkgreen', lw=2)\n",
    "ax.legend()\n",
    "ax.grid()\n",
    "ax.set(title=r'Autocorrelation functions, $\\Gamma_h(n)$',\n",
    "       xlabel=r'Lags between observations, $n$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17ebaf44",
   "metadata": {},
   "source": [
    "## Frequency and the Mean Estimator\n",
    "\n",
    "Consider again the AR(1) process generated by discrete sampling with\n",
    "frequency $ h $. Assume that we have a sample of size $ N $ and\n",
    "we would like to estimate the unconditional mean – in our case the true\n",
    "mean is $ \\mu $.\n",
    "\n",
    "Again, the sample average is an unbiased estimator of the unconditional\n",
    "mean\n",
    "\n",
    "$$\n",
    "\\mathbb{E}[\\bar X_N] = \\frac{1}{N}\\sum_{i = 1}^N \\mathbb{E}[X_i] = \\mathbb{E}[X_0] = \\mu\n",
    "$$\n",
    "\n",
    "The variance of the sample mean is given by\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbb{V}\\left(\\bar X_N\\right) &= \\mathbb{V}\\left(\\frac{1}{N}\\sum_{i = 1}^N X_i\\right) \\\\\n",
    "&= \\frac{1}{N^2} \\left(\\sum_{i = 1}^N \\mathbb{V}(X_i) + 2 \\sum_{i = 1}^{N-1} \\sum_{s = i+1}^N \\text{cov}(X_i, X_s) \\right) \\\\\n",
    "&= \\frac{1}{N^2} \\left( N \\gamma(0) + 2 \\sum_{i=1}^{N-1} i \\cdot \\gamma\\left(h\\cdot (N - i)\\right) \\right) \\\\\n",
    "&= \\frac{1}{N^2} \\left( N \\frac{\\sigma^2}{2\\kappa} + 2 \\sum_{i=1}^{N-1} i \\cdot \\exp(-\\kappa h (N - i)) \\frac{\\sigma^2}{2\\kappa} \\right)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "It is explicit in the above equation that time dependence in the data\n",
    "inflates the variance of the mean estimator through the covariance\n",
    "terms.\n",
    "\n",
    "Moreover, as we can see, a higher sampling frequency—smaller\n",
    "$ h $—makes all the covariance terms larger, everything else being\n",
    "fixed.\n",
    "\n",
    "This implies a relatively slower rate of convergence of the\n",
    "sample average for high-frequency data.\n",
    "\n",
    "Intuitively, stronger dependence across observations for high-frequency data reduces the\n",
    "“information content” of each observation relative to the IID case.\n",
    "\n",
    "We can upper bound the variance term in the following way\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbb{V}(\\bar X_N) &= \\frac{1}{N^2} \\left( N \\sigma^2 + 2 \\sum_{i=1}^{N-1} i \\cdot \\exp(-\\kappa h (N - i)) \\sigma^2 \\right) \\\\\n",
    "&\\leq \\frac{\\sigma^2}{2\\kappa N} \\left(1 + 2 \\sum_{i=1}^{N-1} \\cdot \\exp(-\\kappa h (i)) \\right) \\\\\n",
    "&= \\underbrace{\\frac{\\sigma^2}{2\\kappa N}}_{\\text{IID  case}} \\left(1 + 2 \\frac{1 - \\exp(-\\kappa h)^{N-1}}{1 - \\exp(-\\kappa h)} \\right)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Asymptotically, the term $ \\exp(-\\kappa h)^{N-1} $ vanishes and the\n",
    "dependence in the data inflates the benchmark IID variance by a factor\n",
    "of\n",
    "\n",
    "$$\n",
    "\\left(1 + 2 \\frac{1}{1 - \\exp(-\\kappa h)} \\right)\n",
    "$$\n",
    "\n",
    "This long run factor is larger the higher is the frequency (the smaller\n",
    "is $ h $).\n",
    "\n",
    "Therefore, we expect the asymptotic relative MSEs, $ B $, to change\n",
    "with time-dependent data. We just saw that the mean estimator’s rate is\n",
    "roughly changing by a factor of\n",
    "\n",
    "$$\n",
    "\\left(1 + 2 \\frac{1}{1 - \\exp(-\\kappa h)} \\right)\n",
    "$$\n",
    "\n",
    "Unfortunately, the variance estimator’s MSE is harder to derive.\n",
    "\n",
    "Nonetheless, we can approximate it by using (large sample) simulations,\n",
    "thus getting an idea about how the asymptotic relative MSEs changes in\n",
    "the sampling frequency $ h $ relative to the IID case that we\n",
    "compute in closed form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c4270cb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def sample_generator(h, N, M):\n",
    "    ϕ = (1 - np.exp(-κ * h)) * μ\n",
    "    ρ = np.exp(-κ * h)\n",
    "    s = σ**2 * (1 - np.exp(-2 * κ * h)) / (2 * κ)\n",
    "\n",
    "    mean_uncond = μ\n",
    "    std_uncond = np.sqrt(σ**2 / (2 * κ))\n",
    "\n",
    "    ε_path = stat.norm(0, np.sqrt(s)).rvs((M, N))\n",
    "\n",
    "    y_path = np.zeros((M, N + 1))\n",
    "    y_path[:, 0] = stat.norm(mean_uncond, std_uncond).rvs(M)\n",
    "\n",
    "    for i in range(N):\n",
    "        y_path[:, i + 1] = ϕ + ρ * y_path[:, i] + ε_path[:, i]\n",
    "\n",
    "    return y_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8652f3d2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Generate large sample for different frequencies\n",
    "N_app, M_app = 1000, 30000        # Sample size, number of simulations\n",
    "h_grid = np.linspace(.1, 80, 30)\n",
    "\n",
    "var_est_store = []\n",
    "mean_est_store = []\n",
    "labels = []\n",
    "\n",
    "for h in h_grid:\n",
    "    labels.append(h)\n",
    "    sample = sample_generator(h, N_app, M_app)\n",
    "    mean_est_store.append(np.mean(sample, 1))\n",
    "    var_est_store.append(np.var(sample, 1))\n",
    "\n",
    "var_est_store = np.array(var_est_store)\n",
    "mean_est_store = np.array(mean_est_store)\n",
    "\n",
    "# Save mse of estimators\n",
    "mse_mean = np.var(mean_est_store, 1) + (np.mean(mean_est_store, 1) - μ)**2\n",
    "mse_var = np.var(var_est_store, 1) \\\n",
    "          + (np.mean(var_est_store, 1) - var_uncond)**2\n",
    "\n",
    "benchmark_rate = 2 * var_uncond       # IID case\n",
    "\n",
    "# Relative MSE for large samples\n",
    "rate_h = mse_var / mse_mean\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.plot(h_grid, rate_h, c='darkblue', lw=2,\n",
    "        label=r'large sample relative MSE, $B(h)$')\n",
    "ax.axhline(benchmark_rate, c='k', ls='--', label=r'IID benchmark')\n",
    "ax.set_title('Relative MSE for large samples as a function of sampling \\\n",
    "              frequency \\n MSE($S_N$) relative to MSE($\\\\bar X_N$)')\n",
    "ax.set_xlabel('Sampling frequency, $h$')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2356b960",
   "metadata": {},
   "source": [
    "The above figure illustrates the relationship between the asymptotic\n",
    "relative MSEs and the sampling frequency\n",
    "\n",
    "- We can see that with low-frequency data – large values of $ h $\n",
    "  – the ratio of asymptotic rates approaches the IID case.  \n",
    "- As $ h $ gets smaller – the higher the frequency – the relative\n",
    "  performance of the variance estimator is better in the sense that the\n",
    "  ratio of asymptotic rates gets smaller. That is, as the time\n",
    "  dependence gets more pronounced, the rate of convergence of the mean\n",
    "  estimator’s MSE deteriorates more than that of the variance\n",
    "  estimator.  "
   ]
  }
 ],
 "metadata": {
  "date": 1724218393.5508492,
  "filename": "black_litterman.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Two Modifications of Mean-Variance Portfolio Theory"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}