{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5a7b6130-0706-4016-b0e6-a9bcd6fc2ee1",
   "metadata": {},
   "source": [
    "# <center> Python 4 Finance - Phd Course - Roma 3 University<center>\n",
    "## <center> Financial Time Series <center>\n",
    "### <center> Simulation & risk: MC, bootstrap, (m)VaR, backtesting; ShR testing  <center>\n",
    "\n",
    "This version: April 17, 2026"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d335e60b-4944-441c-8bbe-d631cef2571b",
   "metadata": {},
   "source": [
    "## 6.3 Risk Measures\n",
    "\n",
    "Assets are evaluated according to their profitability (i.e., returns) and their risk. But how do we measure the risk? Here is a list of popular measures you can use. \n",
    "\n",
    "For a given period $M$ of size m, considering an asset $i$ we have: \n",
    "\n",
    "1. **Volatility**: \n",
    "    $$\\sigma_i = \\sqrt{\\frac{1}{m}\\sum_{i=1}^{m}(r_i - \\bar{r})^2},$$\n",
    "    where: \n",
    "        . $r_i$ is the return vector of returns of the asset $i$\n",
    "        . $\\bar{r}$ is the average return computed over the period $M$.\n",
    "    A higher volatility is interpreted as a higher risk of loss.\n",
    "     \n",
    "2. **Beta**: \n",
    "    $$\\beta_i = \\frac{Cov(r_i, r_{mkt})}{Var(r_{mkt})},$$\n",
    "    with $r_{mkt}$ the excess of returns of the market.\n",
    "    \n",
    "    As we saw, beta is a measure of an asset's sensitivity to market movements. See the beta interpretation above.\n",
    "    \n",
    "3. **Value at Risk** (**VaR**)\n",
    "\n",
    "    $$VaR_\\alpha = F^{-1}_{r_i}(\\alpha),$$\n",
    "    where $F^{-1}_{r_i}(\\alpha)$ represents the inverse cumulative distribution function (CDF) of the return distribution, evaluated at the confidence level $\\alpha \\in (0,1)$.\n",
    "\n",
    "    The VaR is defined as a measure of the maximum potential loss an asset could experience within a specified time frame, at a certain level of confidence $1-\\alpha$ (e.g., 99.5\\%). \n",
    "\n",
    "4. **Expected Shortfall (ES)** or **Conditional Value at Risk** (**CVaR**):\n",
    "\n",
    "$$\\text{ES}_\\alpha = \\frac{1}{1-\\alpha}\\int_{-\\infty}^{\\text{VaR}_\\alpha}r_if(r_i)dr_i $$\n",
    "It measures the average loss beyond the VaR. It provides a more comprehensive risk assessment than VaR.\n",
    "\n",
    "\n",
    "Following, you can find a chunck of code with the definition of all these measures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5ace2f91-3c77-4cf4-a951-b6fb0100a8aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. volatility\n",
    "def volatility(returns):\n",
    "    return np.sqrt(np.mean((returns - np.mean(returns)) ** 2))\n",
    "\n",
    "# 2. CAPM beta\n",
    "def beta(asset_returns, market_returns):\n",
    "    covariance = np.cov(asset_returns, market_returns)[0, 1] # take the element in pos 0,1\n",
    "    market_variance = np.var(market_returns)\n",
    "    return covariance / market_variance\n",
    "\n",
    "# 3. VaR\n",
    "def VaR(returns, alpha):\n",
    "    return np.percentile(returns, alpha)\n",
    "\n",
    "# 4. ES\n",
    "def ES(returns, alpha):\n",
    "    var_alpha = VaR(returns, alpha)\n",
    "    return np.mean(returns[returns < var_alpha])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a1f1803-7ac7-4596-9a47-db8fb66489a9",
   "metadata": {},
   "source": [
    "Let me download a financial time series here and download the Fama-French factors (what follows in the next two blocks is just a copy-paste from the AP class). \n",
    "\n",
    "Let's download the time series of an uranium-based company: [Cameco Corp.](https://www.cameco.com/) (CCJ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "adae5897-7457-4d3b-86e8-7e5135b41f14",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "YF.download() has changed argument auto_adjust default to True\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[*********************100%***********************]  1 of 1 completed\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Ticker</th>\n",
       "      <th>CCJ.d</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1996-03-14</th>\n",
       "      <td>5.983490</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-03-15</th>\n",
       "      <td>5.886761</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-03-18</th>\n",
       "      <td>5.914394</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-03-19</th>\n",
       "      <td>5.886761</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-03-20</th>\n",
       "      <td>5.914394</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Ticker         CCJ.d\n",
       "Date                \n",
       "1996-03-14  5.983490\n",
       "1996-03-15  5.886761\n",
       "1996-03-18  5.914394\n",
       "1996-03-19  5.886761\n",
       "1996-03-20  5.914394"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import yfinance as yf\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# download CCJ\n",
    "CCJ = yf.download(\"CCJ\", start=\"1990-12-31\", end=\"2026-04-17\")\n",
    "CCJ_d = CCJ[\"Close\"]\n",
    "CCJ_d = CCJ_d.rename(columns={CCJ_d.columns[0]: \"CCJ.d\"})\n",
    "CCJ_d.index = pd.to_datetime(CCJ_d.index)\n",
    "CCJ_d.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9e7fc37e-e27d-4680-9344-7c059827c22c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 103712 characters]",
      "text/plain": [
       "<Figure size 1920x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "\n",
    "# Create a wider figure\n",
    "fig, ax = plt.subplots(figsize=(16, 5), dpi=120, constrained_layout=True)\n",
    "\n",
    "# Plot closing prices\n",
    "ax.plot(CCJ_d.index, CCJ_d, linewidth=2)\n",
    "\n",
    "# Set one tick per year (January)\n",
    "ax.xaxis.set_major_locator(mdates.MonthLocator(bymonth=1))\n",
    "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n",
    "ax.set_title(\"CCJ Closing Prices\")\n",
    "ax.grid(True, alpha=0.3)\n",
    "plt.setp(ax.get_xticklabels(), rotation=45, ha=\"center\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "afdadfe2-4503-4624-9805-95acc80757e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monthly log returns\n",
    "monthly_prices = CCJ_d.resample(\"ME\").last()\n",
    "lret_m = np.log(monthly_prices / monthly_prices.shift(1)).dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1b81aa27-4742-4d69-bc2a-7cc3e5659cba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Mkt-RF</th>\n",
       "      <th>SMB</th>\n",
       "      <th>HML</th>\n",
       "      <th>RMW</th>\n",
       "      <th>CMA</th>\n",
       "      <th>RF</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1996-01-01</th>\n",
       "      <td>0.0226</td>\n",
       "      <td>-0.0264</td>\n",
       "      <td>0.0046</td>\n",
       "      <td>0.0026</td>\n",
       "      <td>0.0230</td>\n",
       "      <td>0.0043</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-02-01</th>\n",
       "      <td>0.0133</td>\n",
       "      <td>0.0177</td>\n",
       "      <td>-0.0105</td>\n",
       "      <td>0.0016</td>\n",
       "      <td>-0.0189</td>\n",
       "      <td>0.0039</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-03-01</th>\n",
       "      <td>0.0074</td>\n",
       "      <td>0.0156</td>\n",
       "      <td>0.0042</td>\n",
       "      <td>0.0144</td>\n",
       "      <td>-0.0095</td>\n",
       "      <td>0.0039</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-04-01</th>\n",
       "      <td>0.0205</td>\n",
       "      <td>0.0470</td>\n",
       "      <td>-0.0397</td>\n",
       "      <td>-0.0020</td>\n",
       "      <td>-0.0203</td>\n",
       "      <td>0.0046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-05-01</th>\n",
       "      <td>0.0236</td>\n",
       "      <td>0.0323</td>\n",
       "      <td>-0.0090</td>\n",
       "      <td>0.0003</td>\n",
       "      <td>-0.0033</td>\n",
       "      <td>0.0042</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Mkt-RF     SMB     HML     RMW     CMA      RF\n",
       "Date                                                      \n",
       "1996-01-01  0.0226 -0.0264  0.0046  0.0026  0.0230  0.0043\n",
       "1996-02-01  0.0133  0.0177 -0.0105  0.0016 -0.0189  0.0039\n",
       "1996-03-01  0.0074  0.0156  0.0042  0.0144 -0.0095  0.0039\n",
       "1996-04-01  0.0205  0.0470 -0.0397 -0.0020 -0.0203  0.0046\n",
       "1996-05-01  0.0236  0.0323 -0.0090  0.0003 -0.0033  0.0042"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "url = \"https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_5_Factors_2x3_CSV.zip\"\n",
    "\n",
    "# Read raw CSV, skipping metadata\n",
    "ff5 = pd.read_csv(url, skiprows=3)\n",
    "\n",
    "# Keep only rows where the first column looks like YYYYMM (6 digits)\n",
    "ff5 = ff5[ff5.iloc[:,0].astype(str).str.strip().str.len() == 6]\n",
    "\n",
    "# Rename first column\n",
    "ff5 = ff5.rename(columns={ff5.columns[0]: \"Date\"})\n",
    "\n",
    "# Convert YYYYMM → datetime\n",
    "ff5[\"Date\"] = pd.to_datetime(ff5[\"Date\"].astype(str), format=\"%Y%m\")\n",
    "\n",
    "# Set index\n",
    "ff5 = ff5.set_index(\"Date\")\n",
    "\n",
    "# Restrict to 1990 onwards\n",
    "ff5 = ff5.loc[\"1996-01\":]\n",
    "\n",
    "# we have to convert all the column in numeric values\n",
    "ff5 = ff5.apply(lambda s: pd.to_numeric(s.astype(str).str.strip(), errors=\"coerce\"))\n",
    "\n",
    "# now divide factor by 100, they are expressed in % terms\n",
    "ff5 = ff5 /100\n",
    "\n",
    "ff5.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2a082b85-ece2-4bc3-9ebb-42a34df9b972",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CCJ.d</th>\n",
       "      <th>Mkt-RF</th>\n",
       "      <th>SMB</th>\n",
       "      <th>HML</th>\n",
       "      <th>RMW</th>\n",
       "      <th>CMA</th>\n",
       "      <th>RF</th>\n",
       "      <th>e_i</th>\n",
       "      <th>e_mkt</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1996-04-30</th>\n",
       "      <td>0.018824</td>\n",
       "      <td>0.0205</td>\n",
       "      <td>0.0470</td>\n",
       "      <td>-0.0397</td>\n",
       "      <td>-0.0020</td>\n",
       "      <td>-0.0203</td>\n",
       "      <td>0.0046</td>\n",
       "      <td>0.014224</td>\n",
       "      <td>0.0205</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-05-31</th>\n",
       "      <td>-0.057569</td>\n",
       "      <td>0.0236</td>\n",
       "      <td>0.0323</td>\n",
       "      <td>-0.0090</td>\n",
       "      <td>0.0003</td>\n",
       "      <td>-0.0033</td>\n",
       "      <td>0.0042</td>\n",
       "      <td>-0.061769</td>\n",
       "      <td>0.0236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-06-30</th>\n",
       "      <td>-0.071577</td>\n",
       "      <td>-0.0113</td>\n",
       "      <td>-0.0367</td>\n",
       "      <td>0.0233</td>\n",
       "      <td>0.0350</td>\n",
       "      <td>0.0120</td>\n",
       "      <td>0.0040</td>\n",
       "      <td>-0.075577</td>\n",
       "      <td>-0.0113</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-07-31</th>\n",
       "      <td>-0.005333</td>\n",
       "      <td>-0.0596</td>\n",
       "      <td>-0.0378</td>\n",
       "      <td>0.0519</td>\n",
       "      <td>0.0290</td>\n",
       "      <td>0.0268</td>\n",
       "      <td>0.0045</td>\n",
       "      <td>-0.009833</td>\n",
       "      <td>-0.0596</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996-08-31</th>\n",
       "      <td>0.059680</td>\n",
       "      <td>0.0277</td>\n",
       "      <td>0.0258</td>\n",
       "      <td>-0.0082</td>\n",
       "      <td>-0.0038</td>\n",
       "      <td>-0.0247</td>\n",
       "      <td>0.0041</td>\n",
       "      <td>0.055580</td>\n",
       "      <td>0.0277</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               CCJ.d  Mkt-RF     SMB     HML     RMW     CMA      RF  \\\n",
       "Date                                                                   \n",
       "1996-04-30  0.018824  0.0205  0.0470 -0.0397 -0.0020 -0.0203  0.0046   \n",
       "1996-05-31 -0.057569  0.0236  0.0323 -0.0090  0.0003 -0.0033  0.0042   \n",
       "1996-06-30 -0.071577 -0.0113 -0.0367  0.0233  0.0350  0.0120  0.0040   \n",
       "1996-07-31 -0.005333 -0.0596 -0.0378  0.0519  0.0290  0.0268  0.0045   \n",
       "1996-08-31  0.059680  0.0277  0.0258 -0.0082 -0.0038 -0.0247  0.0041   \n",
       "\n",
       "                 e_i   e_mkt  \n",
       "Date                          \n",
       "1996-04-30  0.014224  0.0205  \n",
       "1996-05-31 -0.061769  0.0236  \n",
       "1996-06-30 -0.075577 -0.0113  \n",
       "1996-07-31 -0.009833 -0.0596  \n",
       "1996-08-31  0.055580  0.0277  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# align dates: Ken French monthly dates are typically at the beginning\n",
    "# of the month (YYYY-MM-01). We shift them to end-of-month to match r_i.\n",
    "ff5.index = ff5.index.to_period(\"M\").to_timestamp(\"M\")\n",
    "\n",
    "# join CCJ monthly returns (r_i) with factor data\n",
    "# r_i is already end-of-month\n",
    "XX = lret_m.join(ff5, how=\"inner\").dropna()\n",
    "\n",
    "# compute excess returns\n",
    "# Ford excess return = Ford return - risk-free rate\n",
    "XX[\"e_i\"] = XX[\"CCJ.d\"] - XX[\"RF\"]\n",
    "\n",
    "# market excess return is already (Mkt - RF)\n",
    "XX[\"e_mkt\"] = XX[\"Mkt-RF\"]\n",
    "\n",
    "XX.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "2706d629-a8b1-402e-a205-9ad59a07bb7d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>RM</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Volatility</th>\n",
       "      <td>0.112451</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Beta</th>\n",
       "      <td>1.030300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>VaR (0.5%)</th>\n",
       "      <td>-0.319363</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ES (0.5%)</th>\n",
       "      <td>-0.319590</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  RM\n",
       "Volatility  0.112451\n",
       "Beta        1.030300\n",
       "VaR (0.5%) -0.319363\n",
       "ES (0.5%)  -0.319590"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute these risk measures for our returns\n",
    "alpha = 0.005\n",
    "e = XX['e_i']\n",
    "mkt = XX['Mkt-RF']\n",
    "\n",
    "vol = volatility(e)\n",
    "Beta = beta(e,mkt)\n",
    "VaR_995 = VaR(e,alpha)\n",
    "ES_995 = ES(e,alpha)\n",
    "\n",
    "risk_table = pd.DataFrame({\n",
    "    \"RM\": [vol, Beta, VaR_995, ES_995]\n",
    "}, index =  [\"Volatility\", \"Beta\", \"VaR (0.5%)\", \"ES (0.5%)\"],\n",
    ")\n",
    "risk_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f156af58-469b-472f-800c-56367fff6645",
   "metadata": {},
   "source": [
    "## 6.4 Risk-Adjusted Performance Measures\n",
    "\n",
    "Assets are evaluated using a combination of returns and profitability. This is often referred as risk-weighted performance. \n",
    "\n",
    "In this section you can find a list of the most popular performance metrics. \n",
    "\n",
    "For a given asset $i$ observed over a period $M$ of size $m$, we define: \n",
    "\n",
    "1. **Sharpe Ratio (SR)**:\n",
    "\n",
    "    $$SR = \\frac{r_i - r_f}{\\sigma_i},$$\n",
    "    where r_f is the risk-free rate.\n",
    "    \n",
    "    It compares an asset's average return to its risk, measured by volatility. A higher value indicates a better risk-return tradeoff.\n",
    "2. **Sortino Ratio (SoR)**:\n",
    "    $$\\frac{r_i - r_f}{\\sigma^{d}_i},$$\n",
    "    where $\\sigma^{d}_i$ is the downside risk defined as: \n",
    "    \n",
    "    $$\\sigma^{d}_i = \\sqrt{\\frac{1}{m}\\sum_{i=1}^{m}(\\min(0, r_i - \\bar{r}))^2},$$\n",
    "    namely the standard deviation of the negative returns of the asset $i$, representing the risk associated with potential losses.\n",
    "    \n",
    "    Similar to the SR, but it uses only the std deviation of negative returns to measure risk. This makes it more suitable for assessing downside risk.\n",
    "    \n",
    "3. **Max Drawdown** (**MDD**):\n",
    "\n",
    "    $$MDD_M = \\frac{P_M-L_M}{P_M},$$\n",
    "    where: \n",
    "    - $P_M$ is the peak value before the largest drop occurred during M and\n",
    "    - $L_M$ is the lowest value before new high established.\n",
    "\n",
    "    It represents the maximum loss experienced by an asset relative to its historical peak, measured as the percentage decline from the peak to the lowest point. In other words, it is the maximum percentage decline from a previously achieved high.\n",
    "    \n",
    "4. **Calmar Ratio** (**CR**)\n",
    "    $$CR = \\frac{r_i - r_f}{MDD},$$\n",
    "Similar to the other ratios, but it uses max drawdown in the denominator as opposed to standard deviation. \n",
    "\n",
    "There are several other version of the Sharpe Ratio which are called \"modified Sharpe Ratio\". They use other risk measures (like the VaR or the modified VaR) to weight the excess of returns of the asset.  \n",
    "\n",
    "\n",
    "\n",
    "**Note**: this measures are often presented annualized (i.e., daily x252, monthly x12,...).\n",
    "\n",
    "Define the function we will use to compute this measures and compute them for the returns of CCJ. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "05c5f8ef-22f7-4eeb-a2f5-0d30ab9da77e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def SharpeRatio(ret,rf=0):\n",
    "    mu = np.mean(ret-rf)*252\n",
    "    std = ret.std()*np.sqrt(252)\n",
    "    return mu/std*100\n",
    "\n",
    "def SortinoRatio(ret,rf=0):\n",
    "    mu = np.mean(ret-rf)*252\n",
    "    std_neg = ret[ret<0].std()*np.sqrt(252)\n",
    "    return mu/std_neg*100\n",
    "\n",
    "def MaxDD(ret,rf=0):\n",
    "    comp_ret = (ret-rf+1).cumprod()\n",
    "    peak = ret.expanding(min_periods=1).max() # increase the window size\n",
    "    dd = (comp_ret/peak)-1\n",
    "    return dd.min()\n",
    "\n",
    "def CalmarRatio(ret,rf=0):\n",
    "    mu = np.mean(ret-rf)*252\n",
    "    mdd = MaxDD(ret,rf)\n",
    "    calmar = mu/(mdd*252)\n",
    "    return calmar*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9af5389d-3fca-4546-961d-797756b53494",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Est</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>SR</th>\n",
       "      <td>91.920998</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SoR</th>\n",
       "      <td>159.492558</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MDD</th>\n",
       "      <td>-0.520111</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CR</th>\n",
       "      <td>-1.253680</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Est\n",
       "SR    91.920998\n",
       "SoR  159.492558\n",
       "MDD   -0.520111\n",
       "CR    -1.253680"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# apply these performance measures to our data\n",
    "e = XX['e_i']\n",
    "mkt = XX['Mkt-RF']\n",
    "rf = XX['RF']\n",
    "\n",
    "SR = SharpeRatio(e,rf=0) # e is already r-rf\n",
    "SO = SortinoRatio(e,rf=0)\n",
    "MDD = MaxDD(e)\n",
    "CR = CalmarRatio(e,rf=0)\n",
    "\n",
    "perf_table = pd.DataFrame({\n",
    "    \"Est\": [SR, SO, MDD, CR]\n",
    "}, index =  [\"SR\", \"SoR\", \"MDD\", \"CR\"],\n",
    ")\n",
    "perf_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "447eaa05-ce8d-440b-89e1-cc93563d9a41",
   "metadata": {},
   "source": [
    "## Simulations\n",
    "\n",
    "Generating synthetic data is crucial for quantitative finance, offering the ability to evaluate model behavior using data with known characteristics. This practice finds diverse applications, including testing the functionality of backtesting simulators for accuracy and exploring hypothetical scenarios like simulated crises.\n",
    "\n",
    "In synthetic data generation, one typically employs two main approaches, which can be combined for increased sophistication. The first approach, often referred to as \"**model-driven**,\" involves specifying a theoretical model for asset behavior and conducting simulations based on this model. Stochastic Differential Equations (SDEs) are a common tool in this approach. The second approach, termed \"**data-driven**,\" entails estimating the distribution of, for example, stock returns from real market data and using this distribution to generate simulated historical data for the same asset.\n",
    "\n",
    "## Bootstrap \n",
    "\n",
    "Bootstrap is a **non-parametric technique** to make inference from our data. \n",
    "It is a resampling technique widely used in statistics to estimate the sampling distribution of a statistic by **repeatedly resampling with replacement from the observed data**. This method allows for the estimation of standard errors, confidence intervals, and other statistical properties without relying on strict assumptions about the underlying distribution. \n",
    "\n",
    "Mathematically, let $X_1, X_2, ..., X_n$ be the observed data sample. The basic steps of the bootstrap algorithm are as follows:\n",
    "\n",
    "1. Randomly select $n$ observations from the original sample with replacement to form a resampled dataset, denoted as $X_1^*, X_2^*, ..., X_n^*$.\n",
    "2. Calculate the statistic of interest (e.g., mean, median, standard deviation, VaR,...) based on the resampled dataset.\n",
    "3. Repeat steps 1 and 2 a large number of times (e.g., $B$ times) to obtain a distribution of the statistic.\n",
    "4. Compute the desired estimate (e.g., standard error, confidence interval) based on the distribution of the statistic obtained from the resampling process.\n",
    "\n",
    "\n",
    "\n",
    "### Backtesting: a financial application of the bootstrap technique\n",
    "\n",
    "Backtesting is a vital process in quantitative finance, involving the evaluation of trading strategies or financial models using historical data. It aims to assess the **effectiveness** and **robustness** of these strategies or models by simulating their performance against **past market conditions**.\n",
    "\n",
    "Here there is an application on how backtest the market risk. \n",
    "\n",
    "When we define a risk measure, we are intersted in determining its accuracy. In particular: \n",
    "\n",
    "1. **Unconditional coverage**: This refers to the ability of a risk prediction model to accurately cover actual returns, regardless of the distribution of returns or the presence of tail events. An unconditional coverage test evaluates whether the **total number of observed violations is consistent with theoretical expectations** under the specified distribution.\n",
    "\n",
    "2. **Conditional coverage**: the ability of a risk model to effectively cover actual returns **given that a loss event** (tail event) **occurs**. In other words, it measures whether the risk measure (e.g., the VaR) is accurate when the asset return exceeds the specified confidence level. A conditional coverage test compares the **actual number of violations** (e.g., when actual returns exceed the estimated VaR) with the **expected number of violations** based on the specified **confidence level** (e.g., 5% for a 95% VaR).\n",
    "\n",
    "**Note**: The conditional coverage is consistent with the volatility clustering which characterize the financial time series. \n",
    "\n",
    "Typically, we backtest a model using the Value-at-Risk (VaR), that we've just defined in Python). Let's move on with a fancier and more appropiate risk measure for the financial returns (because we are cool). \n",
    "\n",
    "\n",
    "In the \"modern\" application, we substitute the computation of the VaR with the **modified VaR**, which accounts for the distribution of the asset. Among the several version, we use the **Cornish-Fisher Modification** which adjusts the VaR measure by incorporating **skewness** and **kurtosis** of financial returns. \n",
    "It enhances VaR calculations, particularly in cases where returns do not follow a normal distribution.\n",
    "\n",
    "\n",
    "$$\\text{VaR}_{\\alpha} = \\mu + \\sigma \\cdot z_c,$$\n",
    "$$ z_c = z_{\\alpha} + \\left( z_{\\alpha}^2 - 1 \\right) \\cdot \\frac{s}{6} + \\left( z_{\\alpha}^3 - 3z_{\\alpha} \\right) \\cdot \\frac{k}{24} - \\frac{s^2}{36}(2z_\\alpha^3-5z_\\alpha),$$\n",
    "where: \n",
    "\n",
    "- $\\mu$ the average returns;\n",
    "- $\\sigma$ the volatility of returns;\n",
    "- $\\alpha$ the confidence level; \n",
    "- $z_\\alpha$ is the $\\alpha$ quantile of the standard normal distribution;\n",
    "- $s$ the skewness;\n",
    "- $k$ the excess of kurtosis.\n",
    "\n",
    "I implemented the function for you in Python: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "4248da60-20b7-4cfe-8dcc-18dd84e7d628",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as stats\n",
    "import numpy as np\n",
    "\n",
    "def mVaR(returns, alpha):\n",
    "    # Calculate the mean and standard deviation of the returns\n",
    "    mu = np.mean(returns)\n",
    "    sigma = np.std(returns)\n",
    "    \n",
    "    # Calculate skewness and kurtosis of the returns\n",
    "    skewness = stats.skew(returns)\n",
    "    kurtosis = stats.kurtosis(returns)\n",
    "    \n",
    "    # Calculate the critical value z_alpha\n",
    "    z_alpha = stats.norm.ppf(alpha)\n",
    "    \n",
    "    # Calculate the corrected value z_c\n",
    "    z_c = z_alpha + ((z_alpha**2 - 1) * skewness / 6) + ((z_alpha**3 - 3*z_alpha) * kurtosis / 24) - (sigma**2 / 36) * (2 * z_alpha**3 - 5 * z_alpha)\n",
    "    \n",
    "    # Calculate mVaR\n",
    "    mVaR = mu + sigma * z_c\n",
    "    \n",
    "    return mVaR*252"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0290998-fab5-4c51-af45-1c6c88b53b79",
   "metadata": {},
   "source": [
    "### Let's find an application of the backtest\n",
    "\n",
    "Search for a (quick) research question. Reading around[$\\rightarrow$](https://www.independent.co.uk/tech/nuclear-battery-betavolt-atomic-china-b2476979.html)\n",
    "\n",
    "\n",
    "*Are nuclear energy stocks becoming less risky and closer to conventional energy companies?*\n",
    "\n",
    "Let's analyze the recent past of the largest nuclear public company and compare its performance with the oil counterpart, i.e., [Exxon](https://corporate.exxonmobil.com/).  \n",
    "\n",
    "\n",
    "Get adjusted closing prices from YahooFinance!;)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "c68ea2ab-9488-4918-a48c-93a53573be13",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[*********************100%***********************]  1 of 1 completed\n",
      "[*********************100%***********************]  1 of 1 completed\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CCJ</th>\n",
       "      <th>XOM</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2025-12-23</th>\n",
       "      <td>93.410004</td>\n",
       "      <td>118.629288</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2025-12-24</th>\n",
       "      <td>93.339996</td>\n",
       "      <td>118.430618</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2025-12-26</th>\n",
       "      <td>92.839996</td>\n",
       "      <td>118.321342</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2025-12-29</th>\n",
       "      <td>92.839996</td>\n",
       "      <td>119.731941</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2025-12-30</th>\n",
       "      <td>91.370003</td>\n",
       "      <td>120.188896</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  CCJ         XOM\n",
       "Date                             \n",
       "2025-12-23  93.410004  118.629288\n",
       "2025-12-24  93.339996  118.430618\n",
       "2025-12-26  92.839996  118.321342\n",
       "2025-12-29  92.839996  119.731941\n",
       "2025-12-30  91.370003  120.188896"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import yfinance as yf\n",
    "\n",
    "start_date = \"2000-01-01\"\n",
    "end_date = \"2025-12-31\"\n",
    "tickers = [\"CCJ\", \"XOM\"]\n",
    "\n",
    "adjCP = []\n",
    "\n",
    "for ticker in tickers:\n",
    "    Stock_data = yf.download(ticker, start=start_date, end=end_date)\n",
    "    adjCP.append(Stock_data[\"Close\"])\n",
    "\n",
    "# Use concatenate to build the dataframe\n",
    "adjCP = pd.concat(adjCP, axis=1)\n",
    "# rename columns\n",
    "adjCP.columns = tickers\n",
    "\n",
    "adjCP.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "7059e94a-eea8-41cc-8852-663fb6d93caa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CCJ</th>\n",
       "      <th>XOM</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2000-01-04</th>\n",
       "      <td>-0.012474</td>\n",
       "      <td>-0.019340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-01-05</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.053082</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-01-06</th>\n",
       "      <td>-0.008404</td>\n",
       "      <td>0.050405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-01-07</th>\n",
       "      <td>-0.004228</td>\n",
       "      <td>-0.002938</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-01-10</th>\n",
       "      <td>-0.004246</td>\n",
       "      <td>-0.014080</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 CCJ       XOM\n",
       "Date                          \n",
       "2000-01-04 -0.012474 -0.019340\n",
       "2000-01-05  0.000000  0.053082\n",
       "2000-01-06 -0.008404  0.050405\n",
       "2000-01-07 -0.004228 -0.002938\n",
       "2000-01-10 -0.004246 -0.014080"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate log-returns for each stock\n",
    "lret = np.log(adjCP / adjCP.shift(1)).dropna() # drop the first line\n",
    "lret.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "114ff5b0-5c29-43f4-9510-89837499747d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CCJ</th>\n",
       "      <th>XOM</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Mean</th>\n",
       "      <td>0.060</td>\n",
       "      <td>0.030</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St.Deviation</th>\n",
       "      <td>2.718</td>\n",
       "      <td>1.651</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Skewness</th>\n",
       "      <td>-0.068</td>\n",
       "      <td>-0.082</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kurtosis</th>\n",
       "      <td>4.719</td>\n",
       "      <td>8.623</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Excess.Kurtosis</th>\n",
       "      <td>1.719</td>\n",
       "      <td>5.623</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Min</th>\n",
       "      <td>-20.134</td>\n",
       "      <td>-15.027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Quant5</th>\n",
       "      <td>-4.131</td>\n",
       "      <td>-2.518</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Quant25</th>\n",
       "      <td>-1.385</td>\n",
       "      <td>-0.777</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Median</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.045</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Quant75</th>\n",
       "      <td>1.569</td>\n",
       "      <td>0.878</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Quant95</th>\n",
       "      <td>4.193</td>\n",
       "      <td>2.437</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Max</th>\n",
       "      <td>21.046</td>\n",
       "      <td>15.863</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>JB.stat</th>\n",
       "      <td>6070.652</td>\n",
       "      <td>20258.199</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>JB.pval</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LB.lag1.stat</th>\n",
       "      <td>2.407</td>\n",
       "      <td>31.984</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LB.lag1.pval</th>\n",
       "      <td>0.121</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BP.lag1.stat</th>\n",
       "      <td>2.405</td>\n",
       "      <td>31.970</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BP.lag1.pval</th>\n",
       "      <td>0.121</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>N.obs</th>\n",
       "      <td>6537.000</td>\n",
       "      <td>6537.000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      CCJ        XOM\n",
       "Mean                0.060      0.030\n",
       "St.Deviation        2.718      1.651\n",
       "Skewness           -0.068     -0.082\n",
       "Kurtosis            4.719      8.623\n",
       "Excess.Kurtosis     1.719      5.623\n",
       "Min               -20.134    -15.027\n",
       "Quant5             -4.131     -2.518\n",
       "Quant25            -1.385     -0.777\n",
       "Median              0.000      0.045\n",
       "Quant75             1.569      0.878\n",
       "Quant95             4.193      2.437\n",
       "Max                21.046     15.863\n",
       "JB.stat          6070.652  20258.199\n",
       "JB.pval             0.000      0.000\n",
       "LB.lag1.stat        2.407     31.984\n",
       "LB.lag1.pval        0.121      0.000\n",
       "BP.lag1.stat        2.405     31.970\n",
       "BP.lag1.pval        0.121      0.000\n",
       "N.obs            6537.000   6537.000"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# descriptive stats of the sample\n",
    "from scipy.stats import skew, kurtosis, jarque_bera\n",
    "import statsmodels.api as sm\n",
    "\n",
    "def multi_fun(x):\n",
    "    stat_tab = {\n",
    "        'Mean': round(np.mean(x) * 100,3),\n",
    "        'St.Deviation': round(np.std(x) * 100,3),\n",
    "        'Skewness': round(skew(x),3),\n",
    "        'Kurtosis': round(kurtosis(x),3),\n",
    "        'Excess.Kurtosis': round(kurtosis(x) - 3,3),\n",
    "        'Min': round(np.min(x) * 100,3),\n",
    "        'Quant5': round(np.quantile(x, 0.05) * 100,3),\n",
    "        'Quant25': round(np.quantile(x, 0.25) * 100,3),\n",
    "        'Median': round(np.quantile(x, 0.50) * 100,3),\n",
    "        'Quant75': round(np.quantile(x, 0.75) * 100,3),\n",
    "        'Quant95': round(np.quantile(x, 0.95) * 100,3),\n",
    "        'Max': round(np.max(x) * 100,3),\n",
    "        'JB.stat': round(jarque_bera(x)[0],3),\n",
    "        'JB.pval': round(jarque_bera(x)[1],3),\n",
    "        'LB.lag1.stat' : round(sm.stats.diagnostic.acorr_ljungbox(x, lags=1, boxpierce=False)['lb_stat'][1],3),\n",
    "        'LB.lag1.pval' : round(sm.stats.diagnostic.acorr_ljungbox(x, lags=1, boxpierce=False)['lb_pvalue'][1],3),\n",
    "        'BP.lag1.stat' : round(sm.stats.diagnostic.acorr_ljungbox(x, lags=1, boxpierce=True)['bp_stat'][1],3),\n",
    "        'BP.lag1.pval' : round(sm.stats.diagnostic.acorr_ljungbox(x, lags=1, boxpierce=True)['bp_pvalue'][1],3),\n",
    "        'N.obs': round(len(x))\n",
    "    }\n",
    "    return stat_tab \n",
    "\n",
    "\n",
    "# apply the function to each asset\n",
    "output = lret.apply(multi_fun, axis=0)\n",
    "descr_stat = pd.DataFrame(output.tolist(), index=output.index).T\n",
    "descr_stat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18b5c68f-a440-4825-b3bb-f3c491acb190",
   "metadata": {},
   "source": [
    "#### Compare the performance of the assets\n",
    "\n",
    "For a given period $\\Pi$ and for each asset $i$, we compute: \n",
    "\n",
    "- Cumulative Returns (100USD initial investment);\n",
    "- Sharpe Ratio (ShR);\n",
    "- Valute-at-Risk;\n",
    "- modified VaR (mVaR)\n",
    "- modified Sharpe Ratio (mShR), which is defined as: \n",
    "$$mShR_i = \\frac{r_i-r_f}{mVaR_i},$$\n",
    "\n",
    "I've defined the function as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "d6988d2e-f294-40d0-966b-0480fd311342",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mVaR function\n",
    "def mShR(ret,rf,alpha):\n",
    "    e = ret-rf\n",
    "    mu = np.mean(e)*252\n",
    "    m_VaR = mVaR(e,alpha)\n",
    "    return mu/m_VaR"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8013bc0c-626d-4f43-8205-be5fec362ea9",
   "metadata": {},
   "source": [
    "### Compute the compounded returns (initial investment = 100 USD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "6c142847-2683-4ff9-950b-3044f38667a2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 76172 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# compute the compounded returns (initial investment = 100 USD)\n",
    "def R100(x):\n",
    "    init_inv = 100\n",
    "    R = (1 + x).cumprod() * init_inv\n",
    "    return R\n",
    "\n",
    "R = lret.apply(R100)\n",
    "\n",
    "# Plot R\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(R.index, R['CCJ'], linestyle='-', label = \"CCJ\")\n",
    "plt.plot(R.index, R['XOM'], linestyle='-', label = \"XOM\")\n",
    "plt.title('Cumulative Rets (100 USD)')\n",
    "plt.ylabel('USD')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "1a5f6f1c-d4c3-4d31-9819-85de6be13db0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 86752 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# repeat the analysis for the last ten years!\n",
    "lret16 = lret.loc['2016-01-01':]\n",
    "\n",
    "R2016 = lret16.apply(R100)\n",
    "\n",
    "# Plot R2022\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(R2016.index, R2016['CCJ'], linestyle='-', label = \"CCJ\")\n",
    "plt.plot(R2016.index, R2016['XOM'], linestyle='-', label = \"XOM\")\n",
    "plt.title('Cumulative Rets (100 USD)')\n",
    "plt.ylabel('USD')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "720f9f54-3ec0-4556-9a39-21a94df1958c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 85816 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# repeat the analysis for the last two years!\n",
    "lret22 = lret.loc['2023-01-01':]\n",
    "\n",
    "R2022 = lret22.apply(R100)\n",
    "\n",
    "# Plot R2022\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(R2022.index, R2022['CCJ'], linestyle='-', label = \"CCJ\")\n",
    "plt.plot(R2022.index, R2022['XOM'], linestyle='-', label = \"XOM\")\n",
    "plt.title('Cumulative Rets (100 USD)')\n",
    "plt.ylabel('USD')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77f9cb48-f63c-4cff-afe0-8bdb677683a9",
   "metadata": {},
   "source": [
    "### What about risk measures? What about risk-weighted performance? But, this time we use a rolling window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "ed6261a0-1eda-44a4-bab4-4cf462fac502",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set a rolling window \n",
    "w = 252\n",
    "\n",
    "T = lret.shape[0] # number of obs (= days)\n",
    "nSubsets = T - w # number of subsets (out of sample days) \n",
    "assets = 2 # number of assets\n",
    "alpha = .05 # define the confidence level\n",
    "\n",
    "# define an empty numpy matrix to collect the output (one for each metric)\n",
    "var = np.empty((len(range(nSubsets)),assets)) \n",
    "es = np.empty((len(range(nSubsets)),assets)) \n",
    "shr = np.empty((len(range(nSubsets)),assets)) \n",
    "mvar = np.empty((len(range(nSubsets)),assets)) \n",
    "mshr = np.empty((len(range(nSubsets)),assets)) \n",
    "\n",
    "\n",
    "for j in range(nSubsets):\n",
    "    # select the window observations\n",
    "    subset_j = lret[j:w + j - 1]\n",
    "    \n",
    "    # compute metrics\n",
    "    var[j,] = subset_j.apply(VaR, alpha=alpha)\n",
    "    es[j,] = subset_j.apply(ES, alpha=alpha)\n",
    "    shr[j,] = subset_j.apply(SharpeRatio,rf=0)\n",
    "    mvar[j,] = subset_j.apply(mVaR, alpha=alpha)\n",
    "    mshr[j,] = subset_j.apply(mShR, rf=0, alpha=alpha)\n",
    "\n",
    "# add time index to the metrics and tranform them in dataframes\n",
    "VaR_t = pd.DataFrame(var*252, index=lret.index[w:], columns=tickers)\n",
    "ES_t = pd.DataFrame(es*252, index=lret.index[w:], columns=tickers)\n",
    "ShR_t = pd.DataFrame(shr, index=lret.index[w:], columns=tickers)\n",
    "mVaR_t = pd.DataFrame(mvar, index=lret.index[w:], columns=tickers)\n",
    "mShR_t = pd.DataFrame(mshr, index=lret.index[w:], columns=tickers)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "cacf2d42-0c1d-40fc-9dda-0b3733fb1f03",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 278484 characters]",
      "text/plain": [
       "<Figure size 1200x1000 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## plot the results\n",
    "\n",
    "# Set up the subplots\n",
    "fig, axs = plt.subplots(2, 2, figsize=(12, 10))\n",
    "\n",
    "# VaR\n",
    "axs[0, 0].plot(VaR_t.index, VaR_t['CCJ'], linestyle='-', label = \"CCJ\")\n",
    "axs[0, 0].plot(VaR_t.index, VaR_t['XOM'], linestyle='-', label = \"XOM\")\n",
    "axs[0, 0].set_title(\"VaR\")\n",
    "axs[0, 0].legend()\n",
    "axs[0, 0].grid(True)\n",
    "\n",
    "# ShR\n",
    "axs[0, 1].plot(ShR_t.index, ShR_t['CCJ'], linestyle='-', label = \"CCJ\")\n",
    "axs[0, 1].plot(ShR_t.index, ShR_t['XOM'], linestyle='-', label = \"XOM\")\n",
    "axs[0, 1].set_title(\"ShR\")\n",
    "axs[0, 1].legend()\n",
    "axs[0, 1].grid(True)\n",
    "\n",
    "# mVaR\n",
    "axs[1, 0].plot(mVaR_t.index, mVaR_t['CCJ'], linestyle='-', label = \"CCJ\")\n",
    "axs[1, 0].plot(mVaR_t.index, mVaR_t['XOM'], linestyle='-', label = \"XOM\")\n",
    "axs[1, 0].set_title(\"mVaR\")\n",
    "axs[1, 0].legend()\n",
    "axs[1, 0].grid(True)\n",
    "\n",
    "# mShR\n",
    "axs[1, 1].plot(mShR_t.index, mShR_t['CCJ'], linestyle='-', label = \"CCJ\")\n",
    "axs[1, 1].plot(mShR_t.index, mShR_t['XOM'], linestyle='-', label = \"XOM\")\n",
    "axs[1, 1].set_title(\"mShR\")\n",
    "axs[1, 1].legend()\n",
    "axs[1, 1].grid(True)\n",
    "\n",
    "# Adjust layout and display the plot\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f849eeaf-5af3-4dc2-9299-4d9ac954b527",
   "metadata": {},
   "source": [
    "### Check which model (VaR or mVaR) fits better the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "842e7553-7798-497e-a5d0-749fde9fd7b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 188680 characters]",
      "text/plain": [
       "<Figure size 1200x800 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## represent the results\n",
    "\n",
    "# exclude the first year of observations\n",
    "lret1 = lret[w:]\n",
    "\n",
    "# Set up the subplots\n",
    "fig, axs = plt.subplots(2, 1, figsize=(12, 8))\n",
    "\n",
    "# CCJ\n",
    "axs[0].plot(lret1.index, lret1['CCJ']*252, linestyle='-', label = \"CCJ\")\n",
    "axs[0].plot(VaR_t.index, VaR_t['CCJ'], linestyle='-', label = \"VaR\")\n",
    "axs[0].plot(mVaR_t.index, mVaR_t['CCJ'], linestyle='-', label = \"mVaR\")\n",
    "axs[0].set_title(\"CCJ\")\n",
    "axs[0].legend()\n",
    "axs[0].grid(True)\n",
    "\n",
    "\n",
    "# XOM\n",
    "axs[1].plot(lret1.index, lret1['XOM']*252, linestyle='-', label = \"XOM\")\n",
    "axs[1].plot(VaR_t.index, VaR_t['XOM'], linestyle='-', label = \"VaR\")\n",
    "axs[1].plot(mVaR_t.index, mVaR_t['XOM'], linestyle='-', label = \"mVaR\")\n",
    "axs[1].set_title(\"XOM\")\n",
    "axs[1].legend()\n",
    "axs[1].grid(True)\n",
    "\n",
    "# Adjust layout and display the plot\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5759e899-1fbd-4f4a-bd92-cacba339d178",
   "metadata": {},
   "source": [
    "### Define and compute unconditional and conditional coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "ac5fae1b-c856-4fff-897a-987683d68a9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def uc(obs, VaR):\n",
    "    # Calculate VaR violations\n",
    "    violations = obs*252 < VaR\n",
    "    # Calculate total number of VaR violations\n",
    "    total_violations = np.sum(violations)\n",
    "    # Calculate coverage ratio as the percentage of violations over total periods where VaR was exceeded\n",
    "    uc_ratio = total_violations / len(VaR)\n",
    "    return round(uc_ratio,4)*100\n",
    "\n",
    "\n",
    "def cc(obs, VaR):\n",
    "    # Calculate VaR violations\n",
    "    violations = obs * 252 < VaR\n",
    "    # Count the number of VaR violations\n",
    "    v_count = np.sum(violations)\n",
    "    # Create a vector representing the position of violations\n",
    "    v_pos = np.where(violations)[0]\n",
    "    # Count the number of exceptions on the next day when there was a violation on the current day\n",
    "    cons_viol = np.zeros(v_count)\n",
    "    for i in range(v_count-1):\n",
    "        if v_pos[i + 1] - v_pos[i] == 1:\n",
    "            cons_viol[i] = 1\n",
    "    # Calculate the conditional coverage as the ratio of exceptions on the next day to the total violations\n",
    "    cc_ratio = sum(cons_viol) / v_count if v_count > 0 else 0\n",
    "    return round(cc_ratio,4)*100\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "708c1a62-aea7-4c87-9f50-1e68ab059043",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>UC (VaR_t)</th>\n",
       "      <th>UC (mVaR_t)</th>\n",
       "      <th>CC (VaR_t)</th>\n",
       "      <th>CC (mVaR_t)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>CCJ</th>\n",
       "      <td>0.51</td>\n",
       "      <td>5.14</td>\n",
       "      <td>3.12</td>\n",
       "      <td>9.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>XOM</th>\n",
       "      <td>0.57</td>\n",
       "      <td>5.22</td>\n",
       "      <td>11.11</td>\n",
       "      <td>12.20</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     UC (VaR_t)  UC (mVaR_t)  CC (VaR_t)  CC (mVaR_t)\n",
       "CCJ        0.51         5.14        3.12         9.91\n",
       "XOM        0.57         5.22       11.11        12.20"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute and comment on the coverage of each asset\n",
    "assets = ['CCJ','XOM']\n",
    "uc_data = {'UC (VaR_t)': [], 'UC (mVaR_t)': []}\n",
    "cc_data = {'CC (VaR_t)': [], 'CC (mVaR_t)': []}\n",
    "\n",
    "for asset in assets:\n",
    "    uc_data['UC (VaR_t)'].append(uc(lret1[asset], VaR_t[asset]))\n",
    "    uc_data['UC (mVaR_t)'].append(uc(lret1[asset], mVaR_t[asset]))\n",
    "    cc_data['CC (VaR_t)'].append(cc(lret1[asset], VaR_t[asset]))\n",
    "    cc_data['CC (mVaR_t)'].append(cc(lret1[asset], mVaR_t[asset]))\n",
    "\n",
    "uc_df = pd.DataFrame(uc_data)\n",
    "cc_df = pd.DataFrame(cc_data)\n",
    "\n",
    "RM = pd.concat([uc_df, cc_df], axis=1)\n",
    "RM.index= assets\n",
    "RM"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae3ac479-e764-4d9a-88cc-bdd554d3a5ee",
   "metadata": {},
   "source": [
    "### The Kupiec test \n",
    "\n",
    "The Kupiec test evaluates the accuracy of a Value at Risk (VaR) model by comparing the actual number of VaR violations with the expected number of violations under a given significance level. It quantifies the discrepancy using a test statistic derived from the **binomial distribution**, providing a measure of the model's reliability in capturing extreme events.\n",
    "\n",
    "$$t_{stat} = -2 \\times \\left( T \\times \\log(1 - \\alpha) + \\text{violations} \\times \\log(\\alpha) \\right)$$\n",
    "\n",
    "Here, the **null hypothesis** states that the number of observed violations in the dataset is consistent with what would be expected under the specified VaR model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "30135685-9cb0-4f30-878d-30bacfa20074",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import chi2\n",
    "\n",
    "def kup_test(p,v):\n",
    "    T = len(v)\n",
    "    sv = sum(v)\n",
    "    stat = -2*(T*np.log(p) + sv*np.log(1-p))\n",
    "    p_value = 1 - chi2.cdf(stat, 1)\n",
    "\n",
    "    return pd.DataFrame([stat, p_value]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "dd085e1c-606f-48fa-aa26-7281e8dc74f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "## prepare the dataset to test\n",
    "# concatenate returns, VaR, and mVaR\n",
    "CCJviol = pd.concat([lret1['CCJ']*252, VaR_t['CCJ'],mVaR_t['CCJ']], axis=1)\n",
    "XOMviol = pd.concat([lret1['XOM']*252, VaR_t['XOM'], mVaR_t['XOM']], axis=1)\n",
    "# rename columns\n",
    "CCJviol.columns = [\"PnL\",\"VaR\",\"mVaR\"]\n",
    "XOMviol.columns = [\"PnL\",\"VaR\",\"mVaR\"]\n",
    "# compute the VaR violation (1 or 0)\n",
    "CCJviol[\"ViolVaR\"] = (CCJviol[\"PnL\"] < CCJviol[\"VaR\"]).astype(int)\n",
    "XOMviol[\"ViolVaR\"] = (XOMviol[\"PnL\"] < XOMviol[\"VaR\"]).astype(int)\n",
    "# compute the mVaR violation (1 or 0)\n",
    "CCJviol[\"ViolmVaR\"] = (CCJviol[\"PnL\"] < CCJviol[\"mVaR\"]).astype(int)\n",
    "XOMviol[\"ViolmVaR\"] = (XOMviol[\"PnL\"] < XOMviol[\"mVaR\"]).astype(int)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "269b50fa-31e0-4de7-a67f-c33c321efff4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Stat_VaR</th>\n",
       "      <th>pval_VaR</th>\n",
       "      <th>Stat_mVaR</th>\n",
       "      <th>pval_mVaR</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>CCJ</th>\n",
       "      <td>836.483576</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2579.999759</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>XOM</th>\n",
       "      <td>860.449434</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2609.957082</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Stat_VaR  pval_VaR    Stat_mVaR  pval_mVaR\n",
       "CCJ  836.483576       0.0  2579.999759        0.0\n",
       "XOM  860.449434       0.0  2609.957082        0.0"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute the kupiec test\n",
    "kupCCJ_VaR = kup_test(0.95, CCJviol[\"ViolVaR\"])\n",
    "kupXOM_VaR = kup_test(0.95, XOMviol[\"ViolVaR\"])\n",
    "\n",
    "kupCCJ_mVaR = kup_test(0.95, CCJviol[\"ViolmVaR\"])\n",
    "kupXOM_mVaR = kup_test(0.95, XOMviol[\"ViolmVaR\"])\n",
    "\n",
    "KupiecTestVaR = pd.concat([kupCCJ_VaR,kupXOM_VaR], axis=0)\n",
    "KupiecTestmVaR = pd.concat([kupCCJ_mVaR,kupXOM_mVaR], axis=0)\n",
    "\n",
    "KupiecTest = pd.concat([KupiecTestVaR, KupiecTestmVaR], axis=1)\n",
    "KupiecTest.index= assets\n",
    "KupiecTest.columns = [\"Stat_VaR\",\"pval_VaR\",\"Stat_mVaR\",\"pval_mVaR\"]\n",
    "KupiecTest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3483d110-97b4-46c7-82a6-29c0647c2edf",
   "metadata": {},
   "source": [
    "So, both measures are not adequate. How can we find a statistical preference for one or the other asset?\n",
    "\n",
    "\n",
    "\n",
    "Also, **a model may pass the Kupiec test** but still fail the Christoffersen test, indicating that while the number of violations is correct, they occur in clusters — a clear sign of model misspecification.\n",
    "\n",
    "\n",
    "\n",
    "### The Christoffersen Test (Conditional Coverage)\n",
    "\n",
    "While the Kupiec test evaluates whether the **frequency of VaR violations** is correct, it does not assess whether violations occur **independently over time**.\n",
    "\n",
    "The Christoffersen (1998) test extends this framework by jointly testing:\n",
    "\n",
    "1. Unconditional coverage (correct number of violations)\n",
    "2. Independence (no clustering of violations)\n",
    "\n",
    "Let $v_t \\in {0,1}$ be the violation indicator.\n",
    "\n",
    "We define the following transition counts:\n",
    "\n",
    "* $n_{00}$: no violation followed by no violation\n",
    "* $n_{01}$: no violation followed by violation\n",
    "* $n_{10}$: violation followed by no violation\n",
    "* $n_{11}$: violation followed by violation\n",
    "\n",
    "Transition probabilities:\n",
    "\n",
    "$$\n",
    "\\pi_0 = \\frac{n_{01}}{n_{00} + n_{01}}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\pi_1 = \\frac{n_{11}}{n_{10} + n_{11}}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\pi = \\frac{n_{01} + n_{11}}{T - 1}\n",
    "$$\n",
    "\n",
    "The likelihood ratio statistic for independence is:\n",
    "\n",
    "$$\n",
    "LR_{ind} = -2 \\log \\left( \\frac{(1 - \\pi)^{n_{00} + n_{10}} \\pi^{n_{01} + n_{11}}}{(1 - \\pi_0)^{n_{00}} \\pi_0^{n_{01}} (1 - \\pi_1)^{n_{10}} \\pi_1^{n_{11}}} \\right)\n",
    "$$\n",
    "\n",
    "The conditional coverage test combines:\n",
    "\n",
    "$$\n",
    "LR_{cc} = LR_{uc} + LR_{ind}\n",
    "$$\n",
    "\n",
    "Under the null hypothesis:\n",
    "\n",
    "$$\n",
    "LR_{cc} \\sim \\chi^2(2)\n",
    "$$\n",
    "\n",
    "A rejection implies that the VaR model either has incorrect coverage or exhibits clustering of violations.\n",
    "\n",
    "\n",
    "Let's implement the function here in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "d087b597-3320-4bac-9f62-ff970a81fe7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import chi2\n",
    "\n",
    "def christoffersen_cc_test(alpha, v, eps=1e-12):\n",
    "    \"\"\"\n",
    "    Christoffersen conditional coverage test.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    alpha : float\n",
    "        VaR significance level.\n",
    "    v : array-like\n",
    "        Violation indicator: 1 = violation, 0 = no violation.\n",
    "    eps : float\n",
    "        Small numerical tolerance to avoid log(0).\n",
    "    \"\"\"\n",
    "\n",
    "    v = np.asarray(v).astype(int)\n",
    "    T = len(v)\n",
    "    x = np.sum(v)\n",
    "\n",
    "    # Empirical violation frequency\n",
    "    phat = x / T\n",
    "    phat = np.clip(phat, eps, 1 - eps)\n",
    "    alpha = np.clip(alpha, eps, 1 - eps)\n",
    "\n",
    "    # Unconditional coverage\n",
    "    LR_uc = -2 * (\n",
    "        (T - x) * np.log(1 - alpha) + x * np.log(alpha)\n",
    "        - (T - x) * np.log(1 - phat) - x * np.log(phat)\n",
    "    )\n",
    "\n",
    "    # Transition counts\n",
    "    n00 = n01 = n10 = n11 = 0\n",
    "\n",
    "    for i in range(1, T):\n",
    "        if v[i-1] == 0 and v[i] == 0:\n",
    "            n00 += 1\n",
    "        elif v[i-1] == 0 and v[i] == 1:\n",
    "            n01 += 1\n",
    "        elif v[i-1] == 1 and v[i] == 0:\n",
    "            n10 += 1\n",
    "        elif v[i-1] == 1 and v[i] == 1:\n",
    "            n11 += 1\n",
    "\n",
    "    # Transition probabilities\n",
    "    pi0 = n01 / (n00 + n01) if (n00 + n01) > 0 else eps\n",
    "    pi1 = n11 / (n10 + n11) if (n10 + n11) > 0 else eps\n",
    "    pi = (n01 + n11) / (n00 + n01 + n10 + n11)\n",
    "\n",
    "    pi0 = np.clip(pi0, eps, 1 - eps)\n",
    "    pi1 = np.clip(pi1, eps, 1 - eps)\n",
    "    pi = np.clip(pi, eps, 1 - eps)\n",
    "\n",
    "    # Independence\n",
    "    LR_ind = -2 * (\n",
    "        (n00 + n10) * np.log(1 - pi) + (n01 + n11) * np.log(pi)\n",
    "        - n00 * np.log(1 - pi0) - n01 * np.log(pi0)\n",
    "        - n10 * np.log(1 - pi1) - n11 * np.log(pi1)\n",
    "    )\n",
    "\n",
    "    # Conditional coverage\n",
    "    LR_cc = LR_uc + LR_ind\n",
    "\n",
    "    # p-values\n",
    "    p_value_uc = 1 - chi2.cdf(LR_uc, 1)\n",
    "    p_value_ind = 1 - chi2.cdf(LR_ind, 1)\n",
    "    p_value_cc = 1 - chi2.cdf(LR_cc, 2)\n",
    "\n",
    "    return pd.DataFrame({\n",
    "        \"violations\": [x],\n",
    "        \"expected_violations\": [T * alpha],\n",
    "        \"LR_uc\": [LR_uc],\n",
    "        \"p_value_uc\": [p_value_uc],\n",
    "        \"LR_ind\": [LR_ind],\n",
    "        \"p_value_ind\": [p_value_ind],\n",
    "        \"LR_cc\": [LR_cc],\n",
    "        \"p_value_cc\": [p_value_cc],\n",
    "        \"pi0\": [pi0],\n",
    "        \"pi1\": [pi1],\n",
    "        \"pi\": [pi],\n",
    "        \"n00\": [n00],\n",
    "        \"n01\": [n01],\n",
    "        \"n10\": [n10],\n",
    "        \"n11\": [n11]\n",
    "    })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "21e2cd00-af55-4abf-9a78-d3f7f09cc8f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.005\n",
    "\n",
    "CCJ_VaR_christoffersen = christoffersen_cc_test(alpha, CCJviol[\"ViolVaR\"])\n",
    "CCJ_mVaR_christoffersen = christoffersen_cc_test(alpha, CCJviol[\"ViolmVaR\"])\n",
    "\n",
    "XOM_VaR_christoffersen = christoffersen_cc_test(alpha, XOMviol[\"ViolVaR\"])\n",
    "XOM_mVaR_christoffersen = christoffersen_cc_test(alpha, XOMviol[\"ViolmVaR\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "427d3e1f-633e-4ed4-baf9-8d5f8a960b8e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Asset</th>\n",
       "      <th>Model</th>\n",
       "      <th>violations</th>\n",
       "      <th>expected_violations</th>\n",
       "      <th>p_value_uc</th>\n",
       "      <th>p_value_ind</th>\n",
       "      <th>p_value_cc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CCJ</td>\n",
       "      <td>VaR</td>\n",
       "      <td>32</td>\n",
       "      <td>31.425</td>\n",
       "      <td>0.918344</td>\n",
       "      <td>1.574024e-01</td>\n",
       "      <td>0.366133</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CCJ</td>\n",
       "      <td>mVaR</td>\n",
       "      <td>323</td>\n",
       "      <td>31.425</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>3.434679e-04</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>XOM</td>\n",
       "      <td>VaR</td>\n",
       "      <td>36</td>\n",
       "      <td>31.425</td>\n",
       "      <td>0.423986</td>\n",
       "      <td>3.799197e-05</td>\n",
       "      <td>0.000150</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>XOM</td>\n",
       "      <td>mVaR</td>\n",
       "      <td>328</td>\n",
       "      <td>31.425</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>3.714395e-07</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Asset Model  violations  expected_violations  p_value_uc   p_value_ind  \\\n",
       "0   CCJ   VaR          32               31.425    0.918344  1.574024e-01   \n",
       "1   CCJ  mVaR         323               31.425    0.000000  3.434679e-04   \n",
       "2   XOM   VaR          36               31.425    0.423986  3.799197e-05   \n",
       "3   XOM  mVaR         328               31.425    0.000000  3.714395e-07   \n",
       "\n",
       "   p_value_cc  \n",
       "0    0.366133  \n",
       "1    0.000000  \n",
       "2    0.000150  \n",
       "3    0.000000  "
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_cc = pd.concat([\n",
    "    CCJ_VaR_christoffersen.assign(Asset=\"CCJ\", Model=\"VaR\"),\n",
    "    CCJ_mVaR_christoffersen.assign(Asset=\"CCJ\", Model=\"mVaR\"),\n",
    "    XOM_VaR_christoffersen.assign(Asset=\"XOM\", Model=\"VaR\"),\n",
    "    XOM_mVaR_christoffersen.assign(Asset=\"XOM\", Model=\"mVaR\")\n",
    "], axis=0)\n",
    "\n",
    "# Select and reorder columns\n",
    "results_cc = results_cc[\n",
    "    [\"Asset\", \"Model\",\n",
    "     \"violations\", \"expected_violations\",\n",
    "     \"p_value_uc\", \"p_value_ind\", \"p_value_cc\"]\n",
    "].reset_index(drop=True)\n",
    "\n",
    "results_cc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2b4e3c9-d505-4888-a7c6-a599b8337a3a",
   "metadata": {},
   "source": [
    "The table reports backtesting results for VaR and modified VaR models. \n",
    "The Kupiec test evaluates whether the number of violations is consistent with the target level, while the Christoffersen test also checks for independence. A failure in the independence test suggests clustering of violations, \n",
    "indicating that the model does not fully capture time-varying risk. \n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "### How can we assess whether two assets show different risk-weighted performance?\n",
    "Let's check the mShR. **Problem**. Given two assets $a$ and $b$ and a risk free. ratee $r_f$, the distribution of this test statistic: \n",
    "\n",
    "$$\\hat{\\Delta}_{a,b} = mShR^\\alpha_{a} - mShR^\\alpha_{b} = \\frac{\\bar{r}_{a}- r_{f}}{mVaR^\\alpha_{a}} -\\frac{\\bar{r}_{b}- r_{f}}{mVaR^\\alpha_{b}}$$ \n",
    "is unknown. \n",
    "\n",
    "\n",
    "What can we do when there is not a statistical test? We can use the **bootstrap**!\n",
    "\n",
    "For instance, this is the only methodology we currenty have to test difference in modified Sharpe Ratio!\n",
    "\n",
    "We will procede this way:\n",
    "1. Compute the difference in modified Sharpe Ratio obtaining a series $\\delta$;\n",
    "2. We sample from $\\delta$ using a bootstrap method to get a larger sample;\n",
    "3. Now, this large sample $\\sim$ (asymptotically) normal distribution;\n",
    "4. We can test whether $\\delta \\neq 0$, statistically.\n",
    "5. Our test statistic will be:  \n",
    "$$Z_{mShR}=\\frac{|\\hat{\\Delta}_{a,b}|}{s(\\hat{\\Delta}_{a,b})}, \\hspace{3em}  H_0: \\hat{\\Delta}_{a,b}=0, $$\n",
    "\n",
    "**Note**: this is a simplified methodology, please refer to *Ledoit and Wolf (2008)* and *Ardia and Boudt (2015)* for further detail.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "7d74e737-5088-4a7e-ba98-e002223bfeca",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "def MShRtest_boot(msr1, msr2, B, alpha):\n",
    "    # Calculate differences between the two MSR series\n",
    "    delta = (mShR_t['CCJ'] - mShR_t['XOM']).values    \n",
    "    \n",
    "    # Generate bootstrap samples from the difference series\n",
    "    T = len(delta)\n",
    "    boot_sample = np.zeros(B)\n",
    "    \n",
    "    for i in range(B):\n",
    "        bootstrap_indices = np.random.choice(T, size=T, replace=True)\n",
    "        bootstrap_difference = delta[bootstrap_indices]\n",
    "        boot_sample[i] = np.mean(bootstrap_difference)\n",
    "    \n",
    "    # Test if this variable is different from zero\n",
    "    bootstrap_mean = np.mean(boot_sample)\n",
    "    bootstrap_std = np.std(boot_sample, ddof=1)\n",
    "    Z = bootstrap_mean / bootstrap_std   \n",
    "    \n",
    "    p_value = 2 * (1 - norm.cdf(abs(Z)))\n",
    "    \n",
    "    return pd.DataFrame([round(Z,2), round(p_value,4)]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "d31a8002-eb40-40e5-a6ea-c86fbf3d3142",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Stat</th>\n",
       "      <th>pval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>CCJvsXOM</th>\n",
       "      <td>-2.04</td>\n",
       "      <td>0.0417</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          Stat    pval\n",
       "CCJvsXOM -2.04  0.0417"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# run the test for the four assets\n",
    "CCJvsXOM = MShRtest_boot(mShR_t['CCJ'],mShR_t['XOM'], 10000,0.05)\n",
    "\n",
    "CCJvsXOM.columns = [\"Stat\",\"pval\"]\n",
    "CCJvsXOM.index = [\"CCJvsXOM\"]\n",
    "CCJvsXOM"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bddf947-60fb-47c9-bee2-851b7a60b20f",
   "metadata": {},
   "source": [
    "## Monte Carlo Simulation\n",
    "\n",
    "Monte Carlo simulation is a computational technique used to estimate the probability distribution of outcomes in complex systems or processes by using random sampling methods. It is named after the famous Monte Carlo Casino, known for its games of chance, to reflect the random nature of the technique.\n",
    "\n",
    "Monte Carlo simulation involves several steps:\n",
    "\n",
    "1. **Model Definition**: We choose what to simulate\n",
    "\n",
    "2. **Model Structure**: Model randomness and interconnenction of the variables that we simulate\n",
    "\n",
    "3. **Random Sampling**: Generate (pseudo) random samples from a certain probability distribution\n",
    "\n",
    "4. **Scenarios** description\n",
    "\n",
    "5. **Simulation**: Iteratively calculate the outcomes of the model, often based on a defined set of parameters.\n",
    "\n",
    "6. **Draw synthesis of simulations**.\n",
    "\n",
    "7. **Validation**: Validate the simulation results by comparing them to real-world data or known analytical solutions, if available. \n",
    "\n",
    "8. **Robustness check** (possible)\n",
    "\n",
    "9. **Sensitivity analysis** (possible). Change the parameters, check thee output\n",
    "\n",
    "In essence, we explore various potential scenarios under the assumption that one of them, **probably**, will come to fruition.\n",
    "\n",
    "We use MC simulation to take decision without having observed data. \n",
    "\n",
    "let me show an example: **Roulette!**. We bet on the **even** numbers!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "c91ad11d-c975-4e84-a73f-7bce3df4705b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 472052 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "\n",
    "def roulette_simulation(num_scenarios=100, initial_balance=1000, bet_amount=20, max_hands=100, seed = None):\n",
    "    np.random.seed(seed)  # Fixing the random seed for reproducibility\n",
    "    all_balances = []\n",
    "    for _ in range(num_scenarios):\n",
    "        balance = initial_balance\n",
    "        balances = [balance]  # Store initial balance\n",
    "        hands_played = 0\n",
    "        while balance > 0 and hands_played < max_hands:\n",
    "            # Generate a random number between 0 and 36 to simulate the roulette outcome\n",
    "            outcome = random.randint(0, 36)\n",
    "            # Update the balance based on the outcome of the bet\n",
    "            if outcome == 0:\n",
    "                balance -= bet_amount\n",
    "            elif outcome % 2 == 0:  # Even\n",
    "                balance += bet_amount\n",
    "            else:  # Odd\n",
    "                balance -= bet_amount\n",
    "            balances.append(balance)  # Store updated balance\n",
    "            hands_played += 1\n",
    "        all_balances.append(balances)\n",
    "    \n",
    "    # Plot the results\n",
    "    plt.figure(figsize=(10, 6))\n",
    "    for balances in all_balances:\n",
    "        plt.plot(range(len(balances)), balances)\n",
    "    plt.xlabel('Hand')\n",
    "    plt.ylabel('Balance')\n",
    "    plt.title('Balance Over Hands in Roulette Simulation')\n",
    "    plt.grid(True)\n",
    "    plt.show()\n",
    "\n",
    "# Run the roulette simulation\n",
    "roulette_simulation(num_scenarios=100, initial_balance=1000, bet_amount=20, max_hands=100, seed=1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7235b4ec-af4b-4b9f-b4d4-0ffd504ae1d5",
   "metadata": {},
   "source": [
    "## Monte Carlo Simulation of an Asset Price Path\n",
    "\n",
    "In the Hull framework, a common way to simulate the evolution of an asset price is to assume that prices follow a **Geometric Brownian Motion (GBM)**.\n",
    "\n",
    "The continuous-time model is:\n",
    "\n",
    "$$\n",
    "dS_t = \\mu S_t dt + \\sigma S_t dW_t\n",
    "$$\n",
    "\n",
    "where:\n",
    "\n",
    "* $S_t$ is the asset price at time $t$\n",
    "* $\\mu$ is the expected return\n",
    "* $\\sigma$ is the volatility\n",
    "* $W_t$ is a standard Brownian motion\n",
    "\n",
    "Using a discrete-time approximation, the simulated price dynamics become:\n",
    "\n",
    "$$\n",
    "S_{t+1} = S_t \\exp\\left[\\left(\\mu - \\frac{1}{2}\\sigma^2\\right)\\Delta t + \\sigma \\sqrt{\\Delta t} Z_t\\right]\n",
    "$$\n",
    "\n",
    "where $Z_t \\sim N(0,1)$.\n",
    "\n",
    "This specification ensures that prices remain strictly positive and provides a standard benchmark for Monte Carlo simulation in finance.\n",
    "\n",
    "In the following example, we simulate a single asset price path over 100 trading days.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "2addc38d-9527-4c1c-b9fa-ab1f25cdd362",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 83320 characters]",
      "text/plain": [
       "<Figure size 1800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "\n",
    "# Set the random seed for reproducibility\n",
    "np.random.seed(2806)\n",
    "\n",
    "# Model parameters\n",
    "S0 = 100\n",
    "mu = 0.08\n",
    "sigma = 0.20\n",
    "T = 100\n",
    "dt = 1 / 252\n",
    "\n",
    "# Generate random shocks\n",
    "Z = np.random.normal(0, 1, T)\n",
    "\n",
    "# Initialize the price path\n",
    "S = np.zeros(T + 1)\n",
    "S[0] = S0\n",
    "\n",
    "# Simulate the path\n",
    "for t in range(1, T + 1):\n",
    "    S[t] = S[t-1] * np.exp((mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z[t-1])\n",
    "\n",
    "# Create business-day dates\n",
    "dates = pd.bdate_range(start=\"2026-01-01\", periods=T + 1)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(figsize=(15, 5), dpi=120, constrained_layout=True)\n",
    "ax.plot(dates, S, linewidth=2)\n",
    "\n",
    "ax.set_title(\"Simulated Asset Price Path under Geometric Brownian Motion\")\n",
    "ax.set_ylabel(\"Price\")\n",
    "ax.grid(True, alpha=0.3)\n",
    "\n",
    "# One tick per year at January\n",
    "ax.xaxis.set_major_locator(mdates.MonthLocator(bymonth=1))\n",
    "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n",
    "\n",
    "plt.setp(ax.get_xticklabels(), rotation=0, ha=\"center\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "eb88decf-1090-4a0f-b1a8-ff85b2f8e689",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 1399228 characters]",
      "text/plain": [
       "<Figure size 1800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set the random seed for reproducibility\n",
    "np.random.seed(2806)\n",
    "\n",
    "# Model parameters\n",
    "S0 = 100\n",
    "mu = 0.08\n",
    "sigma = 0.20\n",
    "T = 100\n",
    "dt = 1 / 252\n",
    "n_paths = 1000\n",
    "\n",
    "# Initialize the matrix of simulated prices\n",
    "S = np.zeros((T + 1, n_paths))\n",
    "S[0, :] = S0\n",
    "\n",
    "# Simulate all paths\n",
    "for t in range(1, T + 1):\n",
    "    Z = np.random.normal(0, 1, n_paths)\n",
    "    S[t, :] = S[t-1, :] * np.exp((mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(figsize=(15, 5), dpi=120, constrained_layout=True)\n",
    "ax.plot(S, linewidth=1.2)\n",
    "\n",
    "ax.set_title(\"Monte Carlo Simulation of Multiple Asset Price Paths\")\n",
    "ax.set_xlabel(\"Trading days\")\n",
    "ax.set_ylabel(\"Price\")\n",
    "ax.grid(True, alpha=0.3)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fccb6a54-71d9-46e3-9b9c-fee86ee262d3",
   "metadata": {},
   "source": [
    "## Monte Carlo Simulation to Price Options\n",
    "\n",
    "A typical financial application of Monte Carlo methods is the pricing of **financial derivatives**. These instruments are often difficult to price analytically, and in many cases closed-form solutions do not exist.\n",
    "\n",
    "To address this problem, we rely on **risk-neutral pricing** and simulate **risk-neutral asset paths**.\n",
    "\n",
    "In continuous time, the price of a derivative satisfies:\n",
    "\n",
    "$$\n",
    "C_t = \\mathbb{E}_{\\mathbb{Q}}\\left[ e^{-\\int_t^T r_s ds} C_T \\mid \\mathcal{F}_t \\right]\n",
    "$$\n",
    "\n",
    "where:\n",
    "\n",
    "* $C_t$ is the value of the derivative at time $t$\n",
    "* $C_T$ is the payoff at maturity\n",
    "* $r_s$ is the instantaneous risk-free rate\n",
    "* $\\mathcal{F}_t$ is the information available at time $t$\n",
    "* $\\mathbb{Q}$ is the risk-neutral probability measure\n",
    "\n",
    "Assuming a constant risk-free rate $r$, this simplifies to:\n",
    "\n",
    "$$\n",
    "C_0 = \\mathbb{E}_{\\mathbb{Q}}\\left[ e^{-rT} C_T \\right]\n",
    "$$\n",
    "\n",
    "In a Monte Carlo framework, we approximate this expectation by simulating $M$ independent realizations of the payoff.\n",
    "\n",
    "For each simulation $i$, we compute the discounted payoff:\n",
    "\n",
    "$$\n",
    "C_{0,i} = e^{-rT} C_{T,i}\n",
    "$$\n",
    "\n",
    "The option price is then estimated as the sample average:\n",
    "\n",
    "$$\n",
    "\\hat{C}_0 = \\frac{1}{M} \\sum_{i=1}^{M} C_{0,i}\n",
    "$$\n",
    "\n",
    "This estimator is subject to sampling error. By the Central Limit Theorem, its standard error is given by:\n",
    "\n",
    "$$\n",
    "SE(\\hat{C}_0) = \\frac{\\sigma(C_{0,i})}{\\sqrt{M}}\n",
    "$$\n",
    "\n",
    "where:\n",
    "\n",
    "$$\n",
    "\\sigma(C_{0,i}) = \\sqrt{\\frac{1}{M-1} \\sum_{i=1}^{M} (C_{0,i} - \\hat{C}_0)^2}\n",
    "$$\n",
    "\n",
    "As the number of simulations $M$ increases, the standard error decreases, improving the accuracy of the estimate.\n",
    "\n",
    "\n",
    "**Note**: Under the risk-neutral measure, the expected return of the underlying asset is the risk-free rate, which allows us to simulate price dynamics without estimating risk premia.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c1cf822-40c7-4141-9c0b-340862dc37dc",
   "metadata": {},
   "source": [
    "## European Call Option in the Black-Scholes World\n",
    "\n",
    "We now consider the pricing of a **European call option** in the Black–Scholes framework.\n",
    "\n",
    "In a risk-neutral world with a constant interest rate $r$, the discount factor in continuous time is:\n",
    "\n",
    "$$\n",
    "d = e^{-rT}\n",
    "$$\n",
    "\n",
    "The stock price dynamics are modeled using a **Geometric Brownian Motion (GBM)**:\n",
    "\n",
    "$$\n",
    "dS_t = r S_t dt + \\sigma S_t dW_t\n",
    "$$\n",
    "\n",
    "To simulate the process efficiently, it is convenient to work with the logarithm of the price:\n",
    "\n",
    "$$\n",
    "x_t = \\ln(S_t)\n",
    "$$\n",
    "\n",
    "Using Itô’s lemma, the dynamics of the log-price process become:\n",
    "\n",
    "$$\n",
    "dx_t = \\nu dt + \\sigma dW_t, \\quad \\nu = r - \\frac{1}{2} \\sigma^2\n",
    "$$\n",
    "\n",
    "where:\n",
    "\n",
    "* $\\nu$ is the drift under the risk-neutral measure\n",
    "* $\\sigma$ is the volatility\n",
    "* $dW_t$ is a standard Brownian motion increment\n",
    "\n",
    "By discretizing the process over a small time step $\\Delta t$, we obtain:\n",
    "\n",
    "$$\n",
    "\\Delta x = \\nu \\Delta t + \\sigma \\sqrt{\\Delta t} , \\epsilon_t\n",
    "$$\n",
    "\n",
    "where $\\epsilon_t \\sim N(0,1)$.\n",
    "\n",
    "This leads to the **exact discretization** of the GBM:\n",
    "\n",
    "$$\n",
    "x_{t+\\Delta t} = x_t + \\nu \\Delta t + \\sigma \\sqrt{\\Delta t} , \\epsilon_t\n",
    "$$\n",
    "\n",
    "and therefore:\n",
    "\n",
    "$$\n",
    "S_{t+\\Delta t} = S_t \\exp\\left( \\nu \\Delta t + \\sigma \\sqrt{\\Delta t} , \\epsilon_t \\right)\n",
    "$$\n",
    "\n",
    "Under this model, log-prices are normally distributed and prices remain strictly positive.\n",
    "\n",
    "For each simulated path, the payoff of a European call option is:\n",
    "\n",
    "$$\n",
    "C_T = \\max(S_T - K, 0)\n",
    "$$\n",
    "\n",
    "\n",
    "### Application\n",
    "\n",
    "We apply this framework to price options on CCJ using Monte Carlo simulation.\n",
    "\n",
    "The list of available options on CCJ is available [here](https://finance.yahoo.com/quote/CCJ/options/).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "c7ea7cdc-f1ad-4333-bc07-d8c5519dc03b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expiration Dates: ('2026-04-24', '2026-05-01', '2026-05-08', '2026-05-15', '2026-05-22', '2026-05-29', '2026-06-18', '2026-09-18', '2026-12-18', '2027-01-15', '2028-01-21')\n"
     ]
    }
   ],
   "source": [
    "import yfinance as yf\n",
    "import pandas as pd\n",
    "\n",
    "# Specify the ticker\n",
    "ticker = \"CCJ\"\n",
    "\n",
    "# Get stock data\n",
    "stock = yf.Ticker(ticker)\n",
    "\n",
    "# Fetch all expiration dates for the options\n",
    "expiration_dates = stock.options\n",
    "print(\"Expiration Dates:\", expiration_dates)\n",
    "\n",
    "# Create empty DataFrames to store all call and put data\n",
    "all_calls = pd.DataFrame()\n",
    "all_puts = pd.DataFrame()\n",
    "\n",
    "# Loop through all expiration dates\n",
    "for date in expiration_dates:    \n",
    "    # Fetch the option chain for the current expiration date\n",
    "    option_chain = stock.option_chain(date)\n",
    "    calls = option_chain.calls  # Call options data\n",
    "    puts = option_chain.puts    # Put options data\n",
    "\n",
    "    # Add a column to indicate the expiration date\n",
    "    calls['expirationDate'] = date\n",
    "    puts['expirationDate'] = date\n",
    "\n",
    "    # Append the data to the cumulative DataFrames\n",
    "    all_calls = pd.concat([all_calls, calls], ignore_index=True)\n",
    "    all_puts = pd.concat([all_puts, puts], ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "9bcaa24e-40c2-4fc7-9404-010d071bca7e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>contractSymbol</th>\n",
       "      <th>lastTradeDate</th>\n",
       "      <th>strike</th>\n",
       "      <th>lastPrice</th>\n",
       "      <th>bid</th>\n",
       "      <th>ask</th>\n",
       "      <th>change</th>\n",
       "      <th>percentChange</th>\n",
       "      <th>volume</th>\n",
       "      <th>openInterest</th>\n",
       "      <th>impliedVolatility</th>\n",
       "      <th>inTheMoney</th>\n",
       "      <th>contractSize</th>\n",
       "      <th>currency</th>\n",
       "      <th>expirationDate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CCJ260424C00075000</td>\n",
       "      <td>2026-03-13 18:34:54+00:00</td>\n",
       "      <td>75.0</td>\n",
       "      <td>35.45</td>\n",
       "      <td>44.75</td>\n",
       "      <td>47.50</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.371098</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CCJ260424C00090000</td>\n",
       "      <td>2026-04-15 14:39:57+00:00</td>\n",
       "      <td>90.0</td>\n",
       "      <td>29.75</td>\n",
       "      <td>29.85</td>\n",
       "      <td>32.45</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.599611</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>CCJ260424C00093000</td>\n",
       "      <td>2026-04-15 13:45:21+00:00</td>\n",
       "      <td>93.0</td>\n",
       "      <td>24.30</td>\n",
       "      <td>26.75</td>\n",
       "      <td>29.55</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.453128</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CCJ260424C00095000</td>\n",
       "      <td>2026-04-09 14:12:24+00:00</td>\n",
       "      <td>95.0</td>\n",
       "      <td>21.24</td>\n",
       "      <td>24.50</td>\n",
       "      <td>27.55</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.272465</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>CCJ260424C00096000</td>\n",
       "      <td>2026-04-13 13:34:23+00:00</td>\n",
       "      <td>96.0</td>\n",
       "      <td>19.00</td>\n",
       "      <td>23.90</td>\n",
       "      <td>26.55</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.352542</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>373</th>\n",
       "      <td>CCJ280121C00175000</td>\n",
       "      <td>2026-04-09 14:30:11+00:00</td>\n",
       "      <td>175.0</td>\n",
       "      <td>18.52</td>\n",
       "      <td>19.50</td>\n",
       "      <td>22.40</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>0.565556</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>374</th>\n",
       "      <td>CCJ280121C00180000</td>\n",
       "      <td>2026-04-08 19:42:53+00:00</td>\n",
       "      <td>180.0</td>\n",
       "      <td>19.30</td>\n",
       "      <td>18.50</td>\n",
       "      <td>22.50</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>0.573246</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>375</th>\n",
       "      <td>CCJ280121C00185000</td>\n",
       "      <td>2026-04-13 13:30:06+00:00</td>\n",
       "      <td>185.0</td>\n",
       "      <td>17.16</td>\n",
       "      <td>17.00</td>\n",
       "      <td>21.25</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>0.565465</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>376</th>\n",
       "      <td>CCJ280121C00190000</td>\n",
       "      <td>2026-04-14 16:55:50+00:00</td>\n",
       "      <td>190.0</td>\n",
       "      <td>16.30</td>\n",
       "      <td>16.50</td>\n",
       "      <td>20.50</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>97.0</td>\n",
       "      <td>0.568913</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>377</th>\n",
       "      <td>CCJ280121C00195000</td>\n",
       "      <td>2026-04-10 18:35:47+00:00</td>\n",
       "      <td>195.0</td>\n",
       "      <td>16.27</td>\n",
       "      <td>15.00</td>\n",
       "      <td>19.50</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>73.0</td>\n",
       "      <td>0.561589</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>378 rows × 15 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         contractSymbol             lastTradeDate  strike  lastPrice    bid  \\\n",
       "0    CCJ260424C00075000 2026-03-13 18:34:54+00:00    75.0      35.45  44.75   \n",
       "1    CCJ260424C00090000 2026-04-15 14:39:57+00:00    90.0      29.75  29.85   \n",
       "2    CCJ260424C00093000 2026-04-15 13:45:21+00:00    93.0      24.30  26.75   \n",
       "3    CCJ260424C00095000 2026-04-09 14:12:24+00:00    95.0      21.24  24.50   \n",
       "4    CCJ260424C00096000 2026-04-13 13:34:23+00:00    96.0      19.00  23.90   \n",
       "..                  ...                       ...     ...        ...    ...   \n",
       "373  CCJ280121C00175000 2026-04-09 14:30:11+00:00   175.0      18.52  19.50   \n",
       "374  CCJ280121C00180000 2026-04-08 19:42:53+00:00   180.0      19.30  18.50   \n",
       "375  CCJ280121C00185000 2026-04-13 13:30:06+00:00   185.0      17.16  17.00   \n",
       "376  CCJ280121C00190000 2026-04-14 16:55:50+00:00   190.0      16.30  16.50   \n",
       "377  CCJ280121C00195000 2026-04-10 18:35:47+00:00   195.0      16.27  15.00   \n",
       "\n",
       "       ask  change  percentChange  volume  openInterest  impliedVolatility  \\\n",
       "0    47.50     0.0            0.0     NaN           1.0           2.371098   \n",
       "1    32.45     0.0            0.0     2.0           3.0           1.599611   \n",
       "2    29.55     0.0            0.0     NaN           1.0           1.453128   \n",
       "3    27.55     0.0            0.0     2.0           5.0           1.272465   \n",
       "4    26.55     0.0            0.0     1.0           3.0           1.352542   \n",
       "..     ...     ...            ...     ...           ...                ...   \n",
       "373  22.40     0.0            0.0     1.0          23.0           0.565556   \n",
       "374  22.50     0.0            0.0     1.0          37.0           0.573246   \n",
       "375  21.25     0.0            0.0     1.0          27.0           0.565465   \n",
       "376  20.50     0.0            0.0    20.0          97.0           0.568913   \n",
       "377  19.50     0.0            0.0    11.0          73.0           0.561589   \n",
       "\n",
       "     inTheMoney contractSize currency expirationDate  \n",
       "0          True      REGULAR      USD     2026-04-24  \n",
       "1          True      REGULAR      USD     2026-04-24  \n",
       "2          True      REGULAR      USD     2026-04-24  \n",
       "3          True      REGULAR      USD     2026-04-24  \n",
       "4          True      REGULAR      USD     2026-04-24  \n",
       "..          ...          ...      ...            ...  \n",
       "373       False      REGULAR      USD     2028-01-21  \n",
       "374       False      REGULAR      USD     2028-01-21  \n",
       "375       False      REGULAR      USD     2028-01-21  \n",
       "376       False      REGULAR      USD     2028-01-21  \n",
       "377       False      REGULAR      USD     2028-01-21  \n",
       "\n",
       "[378 rows x 15 columns]"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_calls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "43ceff7b-579a-4edd-ac94-cdd0f22cee55",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>contractSymbol</th>\n",
       "      <th>lastTradeDate</th>\n",
       "      <th>strike</th>\n",
       "      <th>lastPrice</th>\n",
       "      <th>bid</th>\n",
       "      <th>ask</th>\n",
       "      <th>change</th>\n",
       "      <th>percentChange</th>\n",
       "      <th>volume</th>\n",
       "      <th>openInterest</th>\n",
       "      <th>impliedVolatility</th>\n",
       "      <th>inTheMoney</th>\n",
       "      <th>contractSize</th>\n",
       "      <th>currency</th>\n",
       "      <th>expirationDate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CCJ260424P00070000</td>\n",
       "      <td>2026-03-24 13:41:30+00:00</td>\n",
       "      <td>70.0</td>\n",
       "      <td>0.18</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.67</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>7.0</td>\n",
       "      <td>3.001956</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CCJ260424P00075000</td>\n",
       "      <td>2026-04-07 19:55:41+00:00</td>\n",
       "      <td>75.0</td>\n",
       "      <td>0.25</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.67</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>5.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>2.679691</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>CCJ260424P00080000</td>\n",
       "      <td>2026-04-17 17:19:59+00:00</td>\n",
       "      <td>80.0</td>\n",
       "      <td>0.28</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.26</td>\n",
       "      <td>0.14</td>\n",
       "      <td>100.000000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>6954.0</td>\n",
       "      <td>1.683595</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CCJ260424P00085000</td>\n",
       "      <td>2026-04-16 13:30:14+00:00</td>\n",
       "      <td>85.0</td>\n",
       "      <td>0.13</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.75</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>58.0</td>\n",
       "      <td>1.755861</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>CCJ260424P00090000</td>\n",
       "      <td>2026-04-17 14:12:52+00:00</td>\n",
       "      <td>90.0</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.08</td>\n",
       "      <td>-0.01</td>\n",
       "      <td>-16.666664</td>\n",
       "      <td>90.0</td>\n",
       "      <td>104.0</td>\n",
       "      <td>1.078130</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-04-24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>345</th>\n",
       "      <td>CCJ280121P00140000</td>\n",
       "      <td>2026-01-20 17:24:43+00:00</td>\n",
       "      <td>140.0</td>\n",
       "      <td>42.90</td>\n",
       "      <td>45.90</td>\n",
       "      <td>48.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.541478</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>346</th>\n",
       "      <td>CCJ280121P00145000</td>\n",
       "      <td>2026-03-12 18:49:04+00:00</td>\n",
       "      <td>145.0</td>\n",
       "      <td>49.50</td>\n",
       "      <td>43.00</td>\n",
       "      <td>47.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>0.485540</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>347</th>\n",
       "      <td>CCJ280121P00150000</td>\n",
       "      <td>2026-03-05 17:36:32+00:00</td>\n",
       "      <td>150.0</td>\n",
       "      <td>55.40</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000010</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>348</th>\n",
       "      <td>CCJ280121P00155000</td>\n",
       "      <td>2026-01-22 17:08:29+00:00</td>\n",
       "      <td>155.0</td>\n",
       "      <td>52.00</td>\n",
       "      <td>55.90</td>\n",
       "      <td>58.50</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.528661</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>349</th>\n",
       "      <td>CCJ280121P00195000</td>\n",
       "      <td>2026-03-20 19:56:00+00:00</td>\n",
       "      <td>195.0</td>\n",
       "      <td>97.90</td>\n",
       "      <td>80.50</td>\n",
       "      <td>84.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0.432928</td>\n",
       "      <td>True</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2028-01-21</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>350 rows × 15 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         contractSymbol             lastTradeDate  strike  lastPrice    bid  \\\n",
       "0    CCJ260424P00070000 2026-03-24 13:41:30+00:00    70.0       0.18   0.00   \n",
       "1    CCJ260424P00075000 2026-04-07 19:55:41+00:00    75.0       0.25   0.00   \n",
       "2    CCJ260424P00080000 2026-04-17 17:19:59+00:00    80.0       0.28   0.00   \n",
       "3    CCJ260424P00085000 2026-04-16 13:30:14+00:00    85.0       0.13   0.00   \n",
       "4    CCJ260424P00090000 2026-04-17 14:12:52+00:00    90.0       0.05   0.01   \n",
       "..                  ...                       ...     ...        ...    ...   \n",
       "345  CCJ280121P00140000 2026-01-20 17:24:43+00:00   140.0      42.90  45.90   \n",
       "346  CCJ280121P00145000 2026-03-12 18:49:04+00:00   145.0      49.50  43.00   \n",
       "347  CCJ280121P00150000 2026-03-05 17:36:32+00:00   150.0      55.40   0.00   \n",
       "348  CCJ280121P00155000 2026-01-22 17:08:29+00:00   155.0      52.00  55.90   \n",
       "349  CCJ280121P00195000 2026-03-20 19:56:00+00:00   195.0      97.90  80.50   \n",
       "\n",
       "       ask  change  percentChange  volume  openInterest  impliedVolatility  \\\n",
       "0     1.67    0.00       0.000000     NaN           7.0           3.001956   \n",
       "1     1.67    0.00       0.000000     5.0           8.0           2.679691   \n",
       "2     0.26    0.14     100.000000     1.0        6954.0           1.683595   \n",
       "3     0.75    0.00       0.000000     1.0          58.0           1.755861   \n",
       "4     0.08   -0.01     -16.666664    90.0         104.0           1.078130   \n",
       "..     ...     ...            ...     ...           ...                ...   \n",
       "345  48.00    0.00       0.000000     2.0           1.0           0.541478   \n",
       "346  47.00    0.00       0.000000     1.0           4.0           0.485540   \n",
       "347   0.00    0.00       0.000000     1.0           0.0           0.000010   \n",
       "348  58.50    0.00       0.000000     NaN           2.0           0.528661   \n",
       "349  84.00    0.00       0.000000     3.0           3.0           0.432928   \n",
       "\n",
       "     inTheMoney contractSize currency expirationDate  \n",
       "0         False      REGULAR      USD     2026-04-24  \n",
       "1         False      REGULAR      USD     2026-04-24  \n",
       "2         False      REGULAR      USD     2026-04-24  \n",
       "3         False      REGULAR      USD     2026-04-24  \n",
       "4         False      REGULAR      USD     2026-04-24  \n",
       "..          ...          ...      ...            ...  \n",
       "345        True      REGULAR      USD     2028-01-21  \n",
       "346        True      REGULAR      USD     2028-01-21  \n",
       "347        True      REGULAR      USD     2028-01-21  \n",
       "348        True      REGULAR      USD     2028-01-21  \n",
       "349        True      REGULAR      USD     2028-01-21  \n",
       "\n",
       "[350 rows x 15 columns]"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_puts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "910fcb7d-ac1a-4773-a19d-260e89f2b2a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>contractSymbol</th>\n",
       "      <th>lastTradeDate</th>\n",
       "      <th>strike</th>\n",
       "      <th>lastPrice</th>\n",
       "      <th>bid</th>\n",
       "      <th>ask</th>\n",
       "      <th>change</th>\n",
       "      <th>percentChange</th>\n",
       "      <th>volume</th>\n",
       "      <th>openInterest</th>\n",
       "      <th>impliedVolatility</th>\n",
       "      <th>inTheMoney</th>\n",
       "      <th>contractSize</th>\n",
       "      <th>currency</th>\n",
       "      <th>expirationDate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>100</th>\n",
       "      <td>CCJ260501C00160000</td>\n",
       "      <td>2026-03-17 16:08:18+00:00</td>\n",
       "      <td>160.0</td>\n",
       "      <td>0.46</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.96</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.939454</td>\n",
       "      <td>False</td>\n",
       "      <td>REGULAR</td>\n",
       "      <td>USD</td>\n",
       "      <td>2026-05-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         contractSymbol             lastTradeDate  strike  lastPrice  bid  \\\n",
       "100  CCJ260501C00160000 2026-03-17 16:08:18+00:00   160.0       0.46  0.0   \n",
       "\n",
       "      ask  change  percentChange  volume  openInterest  impliedVolatility  \\\n",
       "100  0.96     0.0            0.0     NaN           1.0           0.939454   \n",
       "\n",
       "     inTheMoney contractSize currency expirationDate  \n",
       "100       False      REGULAR      USD     2026-05-01  "
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# take a random call and compute the price\n",
    "call_random = all_calls.iloc[[100]]\n",
    "call_random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "d64685c9-7f55-4417-ad67-80ab69c4e786",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Call value is $0.61 with SE +/- 0.12\n",
      "Observed market value: $0.46\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Compute time to maturity\n",
    "expiration = pd.to_datetime(call_random[\"expirationDate\"])\n",
    "today = pd.to_datetime(CCJ_d.index[-1])\n",
    "\n",
    "TtM_days = (expiration - today).dt.days.item()\n",
    "T = TtM_days / 365\n",
    "\n",
    "# Initial option parameters\n",
    "S0 = CCJ_d.iloc[-1].item()                         # current stock price\n",
    "K = call_random[\"strike\"].item()                   # strike price\n",
    "sigma = call_random[\"impliedVolatility\"].item()    # implied volatility\n",
    "r = 0.04                                           # risk-free rate\n",
    "M = 1000                                           # number of simulations\n",
    "\n",
    "market_value = call_random[\"lastPrice\"].item()     # observed market price\n",
    "\n",
    "# Simulate terminal stock prices under the risk-neutral measure\n",
    "Z = np.random.normal(0, 1, M)\n",
    "ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)\n",
    "\n",
    "# Compute call payoffs\n",
    "CT = np.maximum(ST - K, 0)\n",
    "\n",
    "# Discount the average payoff\n",
    "C0 = np.exp(-r * T) * np.mean(CT)\n",
    "\n",
    "# Compute Monte Carlo standard error\n",
    "SE = np.exp(-r * T) * np.std(CT, ddof=1) / np.sqrt(M)\n",
    "\n",
    "print(\"Call value is ${0} with SE +/- {1}\".format(np.round(C0, 2), np.round(SE, 2)))\n",
    "print(\"Observed market value: ${0}\".format(np.round(market_value, 2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "122e8285-8548-474b-8961-adfce29f6c66",
   "metadata": {},
   "source": [
    "### BONUS SECTION: Option Trading with Python: \n",
    "\n",
    "We could immplement these strategies together: \n",
    "\n",
    "https://zodiactrading.medium.com/options-trading-with-python-5-easy-strategies-with-codes-a4f560e60f70\n",
    "\n",
    "Here is a function to compute the option premium using B-S formulae:\n",
    "\n",
    "https://gist.github.com/jongan69/23e8cfee231819ebe78741a9a7508ec1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c67072f7-26c8-4440-acf6-82c68b9a88b0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
