{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Expectations\n",
    "\n",
    "**Functions**\n",
    "\n",
    "`np.random.RandomState`, `RandomState.standard_normal`, `RandomState.standard_t`, `RandomState.chi2`,\n",
    "`np.exp`, `np.mean`, `np.std`, `scipy.integrate.quadrature`, `scipy.integrate.quad`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 11\n",
    "\n",
    "Compute $E\\left[X\\right]$, $E\\left[X^{2}\\right]$, $V\\left[X\\right]$ and the kurtosis of $X$ using Monte Carlo integration when $X$ is distributed:\n",
    "\n",
    "1. Standard Normal\n",
    "2. $N\\left(0.08,0.2^{2}\\right)$\n",
    "3. Students $t_{8}$\n",
    "4. $\\chi_{5}^{2}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:22.853548Z",
     "iopub.status.busy": "2023-09-28T12:33:22.853548Z",
     "iopub.status.idle": "2023-09-28T12:33:23.013031Z",
     "shell.execute_reply": "2023-09-28T12:33:23.013031Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "rs = np.random.RandomState(30092019)\n",
    "# More modern way, only for NumPy >= 1.17\n",
    "rs = np.random.default_rng(30092019)\n",
    "\n",
    "reps = 10000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Standard Normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.013031Z",
     "iopub.status.busy": "2023-09-28T12:33:23.013031Z",
     "iopub.status.idle": "2023-09-28T12:33:23.028054Z",
     "shell.execute_reply": "2023-09-28T12:33:23.028054Z"
    },
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[X] = 0.007148609933135336, E[X**2]=0.9905771409743099, V[X]=0.9905260383503338, K[X]=2.9809135609763335\n"
     ]
    }
   ],
   "source": [
    "x = rs.standard_normal(reps)\n",
    "mu = x.mean()\n",
    "mu2 = (x**2).mean()\n",
    "var = mu2 - mu**2\n",
    "kurt = np.mean(((x - mu) ** 4)) / var**2\n",
    "\n",
    "print(f\"E[X] = {mu}, E[X**2]={mu2}, V[X]={var}, K[X]={kurt}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$t_8$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.028054Z",
     "iopub.status.busy": "2023-09-28T12:33:23.028054Z",
     "iopub.status.idle": "2023-09-28T12:33:23.036516Z",
     "shell.execute_reply": "2023-09-28T12:33:23.036516Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[X] = 0.00942806090078673, E[X**2]=1.3236881368102587, V[X]=1.3235992484779098, K[X]=4.182677930856709\n"
     ]
    }
   ],
   "source": [
    "x = rs.standard_t(8, size=reps)\n",
    "mu = x.mean()\n",
    "mu2 = (x**2).mean()\n",
    "var = mu2 - mu**2\n",
    "kurt = np.mean(((x - mu) ** 4)) / var**2\n",
    "\n",
    "print(f\"E[X] = {mu}, E[X**2]={mu2}, V[X]={var}, K[X]={kurt}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\chi^2_5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.036516Z",
     "iopub.status.busy": "2023-09-28T12:33:23.036516Z",
     "iopub.status.idle": "2023-09-28T12:33:23.045098Z",
     "shell.execute_reply": "2023-09-28T12:33:23.045098Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[X] = 4.971433291129635, E[X**2]=34.524514021866466, V[X]=9.80936505371443, K[X]=4.943942341459135\n"
     ]
    }
   ],
   "source": [
    "x = rs.chisquare(5, size=reps)\n",
    "mu = x.mean()\n",
    "mu2 = (x**2).mean()\n",
    "var = mu2 - mu**2\n",
    "kurt = np.mean(((x - mu) ** 4)) / var**2\n",
    "\n",
    "print(f\"E[X] = {mu}, E[X**2]={mu2}, V[X]={var}, K[X]={kurt}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$N(8\\%, 20\\%^2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.045098Z",
     "iopub.status.busy": "2023-09-28T12:33:23.045098Z",
     "iopub.status.idle": "2023-09-28T12:33:23.053312Z",
     "shell.execute_reply": "2023-09-28T12:33:23.053312Z"
    },
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[X] = 0.08005584289746774, E[X**2]=0.04645194030328675, V[X]=0.04004300232126271, K[X]=3.039210738211553\n"
     ]
    }
   ],
   "source": [
    "x = rs.normal(0.08, 0.2, size=reps)\n",
    "mu = x.mean()\n",
    "mu2 = (x**2).mean()\n",
    "var = mu2 - mu**2\n",
    "kurt = np.mean(((x - mu) ** 4)) / var**2\n",
    "\n",
    "print(f\"E[X] = {mu}, E[X**2]={mu2}, V[X]={var}, K[X]={kurt}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Function are useful for avoiding many blocks of repetitive code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.053312Z",
     "iopub.status.busy": "2023-09-28T12:33:23.053312Z",
     "iopub.status.idle": "2023-09-28T12:33:23.120576Z",
     "shell.execute_reply": "2023-09-28T12:33:23.120576Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[X] = 4.997835444347266, E[X**2]=34.95855164293982, V[X]=9.980192514165985, K[X]=5.3632198733573935\n"
     ]
    }
   ],
   "source": [
    "def expectations(x):\n",
    "    mu = x.mean()\n",
    "    mu2 = (x**2).mean()\n",
    "    var = mu2 - mu**2\n",
    "    kurt = np.mean(((x - mu) ** 4)) / var**2\n",
    "\n",
    "    print(f\"E[X] = {mu}, E[X**2]={mu2}, V[X]={var}, K[X]={kurt}\")\n",
    "\n",
    "\n",
    "reps = 1000000\n",
    "expectations(rs.chisquare(5, reps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 12 \n",
    "\n",
    "1. Compute $E\\left[\\exp\\left(X\\right)\\right]$ when $X\\sim N\\left(0.08,0.2^{2}\\right)$.\n",
    "2. Compare this to the analytical result for a Log-Normal random variable.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.120576Z",
     "iopub.status.busy": "2023-09-28T12:33:23.120576Z",
     "iopub.status.idle": "2023-09-28T12:33:23.147484Z",
     "shell.execute_reply": "2023-09-28T12:33:23.147484Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[exp(X)]=1.1048650191534506\n"
     ]
    }
   ],
   "source": [
    "x = rs.normal(0.08, 0.2, size=reps)\n",
    "mu_exp = np.mean(np.exp(x))\n",
    "print(f\"E[exp(X)]={mu_exp}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.147484Z",
     "iopub.status.busy": "2023-09-28T12:33:23.147484Z",
     "iopub.status.idle": "2023-09-28T12:33:23.154214Z",
     "shell.execute_reply": "2023-09-28T12:33:23.154214Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Analytical = 1.1051709180756477\n"
     ]
    }
   ],
   "source": [
    "analytical = np.exp(0.08 + 0.2**2 / 2)\n",
    "print(f\"Analytical = {analytical}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 13\n",
    "\n",
    "Explore the role of uncertainty in Monte Carlo integration by increasing the number of simulations 300% relative to the base case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.154214Z",
     "iopub.status.busy": "2023-09-28T12:33:23.154214Z",
     "iopub.status.idle": "2023-09-28T12:33:23.210739Z",
     "shell.execute_reply": "2023-09-28T12:33:23.210739Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Std. Dev (base): 0.008612962366106433\n",
      "Std. Dev (4*base): 0.005737957392241197\n",
      "Ratio: 1.5010502479075196\n"
     ]
    }
   ],
   "source": [
    "base_reps = 10000\n",
    "\n",
    "x = rs.standard_normal((base_reps, 100))\n",
    "mus = x.mean(0)\n",
    "std = mus.std()\n",
    "\n",
    "x = rs.standard_normal((4 * base_reps, 100))\n",
    "mus_4 = x.mean(0)\n",
    "std_4 = mus_4.std()\n",
    "\n",
    "print(f\"Std. Dev (base): {std}\")\n",
    "print(f\"Std. Dev (4*base): {std_4}\")\n",
    "print(f\"Ratio: {std/std_4}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 14\n",
    "\n",
    "Compute the $N(8\\%, 20\\%^2)$  expectation in exercise 11 using quadrature.\n",
    "\n",
    "**Note**: This requires writing a function which will return $\\exp\\left(x\\right)\\times\\phi\\left(x\\right)$ where $\\phi\\left(x\\right)$ is the pdf evaluated at $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.210739Z",
     "iopub.status.busy": "2023-09-28T12:33:23.210739Z",
     "iopub.status.idle": "2023-09-28T12:33:23.977125Z",
     "shell.execute_reply": "2023-09-28T12:33:23.977125Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quadrature: 1.1051649239632353\n"
     ]
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "from scipy.integrate import quadrature, romberg\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    return np.exp(x) * stats.norm.pdf(x, 0.08, 0.2)\n",
    "\n",
    "\n",
    "res, err = quadrature(f, -5 * 0.2, 5 * 0.2)\n",
    "print(f\"Quadrature: {res}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-28T12:33:23.977125Z",
     "iopub.status.busy": "2023-09-28T12:33:23.977125Z",
     "iopub.status.idle": "2023-09-28T12:33:23.996222Z",
     "shell.execute_reply": "2023-09-28T12:33:23.996222Z"
    },
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quadrature (romberg): 1.1051649244195154\n"
     ]
    }
   ],
   "source": [
    "res = romberg(f, -5 * 0.2, 5 * 0.2)\n",
    "print(f\"Quadrature (romberg): {res}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 15 \n",
    "\n",
    "**Optional** (Much more challenging)\n",
    "\n",
    "Suppose log stock market returns are distributed according to a Students t with 8 degrees of\n",
    "freedom, mean 8% and volatility 20%. Utility maximizers hold a portfolio consisting of a\n",
    "risk-free asset paying 1% and the stock market. Assume that they are myopic and only care\n",
    "about next period wealth, so that \n",
    "\n",
    "$$U\\left(W_{t+1}\\right)=U\\left(\\exp\\left(r_{p}\\right)W_{t}\\right)$$\n",
    "\n",
    "and that $U\\left(W\\right)=\\frac{W^{1-\\gamma}}{1-\\gamma}$ is CRRA with risk aversion $\\gamma$.\n",
    "The portfolio return is $r_{p}=wr_{s}+\\left(1-w\\right)r_{f}$ where $s$ is for stock market\n",
    "and $f$ is for risk-free. A 4th order expansion of this utility around the expected wealth\n",
    "next period is\n",
    "\n",
    "$$E_{t}\\left[U\\left(W_{t+1}\\right)\\right]\\approx\\phi_{0}+\\phi_{1}\\mu_{1}^{\\prime}+\\phi_{2}\\mu_{2}^{\\prime}+\\phi_{3}\\mu_{3}^{\\prime}+\\phi_{4}\\mu_{4}^{\\prime}$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\phi_{j}=\\left(j!\\right)^{-1}U^{\\left(j\\right)}\\left(E_{t}\\left[W_{t+1}\\right]\\right),$$\n",
    "\n",
    "$$U^{(j)}=\\frac{\\partial^{j}U}{\\partial W^{j}},$$\n",
    "\n",
    "$$\\mu_{k}^{\\prime}=E_{t}\\left[\\left(r-\\mu\\right)_{p}^{k}\\right],$$\n",
    "\n",
    "and $\\mu=E_{t}\\left[r_{p}\\right]$. Use Monte Carlo integration to examine how the weight in\n",
    "the stock market varies as the risk aversion varies from 1.5 to 10. Note that when $\\gamma=1$, $U\\left(W\\right)=\\ln\\left(W\\right)$.\n",
    "Use $W_{t}=1$ without loss of generality since the portfolio problem is homogeneous of degree 0 in wealth."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}