SPX.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import cvxpy as cp\n",
    "from scipy.stats import norm\n",
    "np.set_printoptions(precision=5, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IvyDB contains a complete historical record of end-of-day data on all US exchange-traded equity and index options (including options on ETFs and ADRs) from January 1996 onward. The data includes both daily option pricing information (symbol, date, closing bid and ask quote, volume, and open interest) as well as high, low, and closing prices for the underlying equity or index. IvyDB also provides all interest rate, dividend, and corporate action information for each security, so you can correlate your own option pricing models with calculations.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 5999 5999 2744.45\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(108, 49)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load SPX options data\n",
    "price_df = pd.read_csv('211c7adb142fa896.csv', parse_dates=['date'])\n",
    "df = pd.read_csv('0b9c74e2057156a8.csv', parse_dates=['date', 'exdate', 'last_date'])\n",
    "df = pd.merge(df, price_df[['date', 'close']])\n",
    "\n",
    "derivs = df[(df['date'] == '2019-06-03') & (df['exdate'] == '2019-06-07')]\n",
    "price = derivs.close.iloc[0]\n",
    "derivs = derivs[derivs.volume > 10]\n",
    "derivs = derivs[np.abs(derivs.strike_price / 1000 - price) / price < .3]\n",
    "\n",
    "calls = derivs[derivs.cp_flag == 'C']\n",
    "puts = derivs[derivs.cp_flag == 'P']\n",
    "p = np.arange(1, 6000, 1)\n",
    "print (p.min(), p.max(), p.size, price)\n",
    "\n",
    "calls_bid = np.array(calls.best_bid)\n",
    "calls_ask = np.array(calls.best_offer)\n",
    "calls_strike = np.array(calls.strike_price / 1000)\n",
    "\n",
    "puts_bid = np.array(puts.best_bid)\n",
    "puts_ask = np.array(puts.best_offer)\n",
    "puts_strike = np.array(puts.strike_price / 1000)\n",
    "\n",
    "puts_strike.size, calls_strike.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = [\"buy put %.3f for %.3f\" % (s, p) for p, s in zip(puts_ask, puts_strike)] + \\\n",
    "        [\"write put %.3f for %.3f\" % (s, p) for p, s in zip(puts_bid, puts_strike)] + \\\n",
    "        [\"buy call %.3f for %.3f\" % (s, p) for p, s in zip(calls_ask, calls_strike)] + \\\n",
    "        [\"write call %.3f for %.3f\" % (s, p) for p, s in zip(calls_bid, calls_strike)] + \\\n",
    "        [\"buy underlying\", \"short underlying\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "option_fee = .65/100\n",
    "short_fee = .0003\n",
    "long_fee = .0003\n",
    "P = np.hstack([\n",
    "    np.maximum(puts_strike[None, :] - p[:, None], 0) - option_fee,\n",
    "    -np.maximum(puts_strike[None, :] - p[:, None], 0) - option_fee,\n",
    "    np.maximum(p[:, None] - calls_strike[None, :], 0) - option_fee,\n",
    "    -np.maximum(p[:, None] - calls_strike[None, :], 0) - option_fee,\n",
    "    p[:, None],\n",
    "    -p[:, None]\n",
    "])\n",
    "d = np.concatenate([\n",
    "    puts_ask,\n",
    "    -puts_bid,\n",
    "    calls_ask,\n",
    "    -calls_bid,\n",
    "    np.array([price + long_fee * price]),\n",
    "    np.array([-price + short_fee * price])\n",
    "])\n",
    "m, n = P.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if arbitrage exists\n",
    "R = P - d[None,:]\n",
    "x = cp.Variable(R.shape[1])\n",
    "prob = cp.Problem(cp.Minimize(0), [x >= 0, cp.min(R @ x) >= 1])\n",
    "assert prob.solve() == np.inf, \"Arbitrage exists\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum entropy / min dirichlet energy, total variation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi = cp.Variable(m)\n",
    "\n",
    "dirichlet = cp.sum_squares(cp.diff(pi))\n",
    "entr = cp.sum(cp.entr(pi))\n",
    "tv = cp.norm1(cp.diff(pi, k=1))\n",
    "\n",
    "cons = [P.T @ pi <= d, cp.sum(pi) == 1, pi >= 0]\n",
    "\n",
    "plt.close()\n",
    "cp.Problem(cp.Minimize(dirichlet), cons).solve(solver=cp.MOSEK)\n",
    "plt.plot(p, pi.value, label='entr', c='k')\n",
    "plt.xlabel(\"$p$\")\n",
    "plt.savefig(\"../price_distr_derivs/figs/dirichlet.pdf\")\n",
    "\n",
    "plt.close()\n",
    "cp.Problem(cp.Maximize(entr), cons).solve(solver=cp.MOSEK)\n",
    "pi_entr = np.copy(pi.value)\n",
    "plt.plot(p, pi.value, label='entr', c='k')\n",
    "plt.xlabel(\"$p$\")\n",
    "plt.savefig(\"../price_distr_derivs/figs/entr.pdf\")\n",
    "\n",
    "plt.close()\n",
    "cp.Problem(cp.Minimize(tv), cons).solve(solver=cp.MOSEK)\n",
    "plt.plot(p, pi.value, label='entr', c='k')\n",
    "plt.xlabel(\"$p$\")\n",
    "plt.savefig(\"../price_distr_derivs/figs/tv.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Closest log-normal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pilognormal = np.ones(m) / m\n",
    "for k in range(5):\n",
    "    mu = pilognormal @ np.log(p)\n",
    "    sigma = np.sqrt(pi_entr @ np.square(np.log(p) - mu))\n",
    "    nu = np.append(0, np.diff(norm.cdf(np.log(p), loc=mu, scale=sigma)))\n",
    "    nu /= nu.sum()\n",
    "    result = cp.Problem(cp.Minimize(cp.sum(cp.kl_div(pi, nu))), cons).solve(solver=cp.MOSEK)\n",
    "    pilognormal = pi.value\n",
    "    print (k, result, mu, sigma)\n",
    "mu = pilognormal @ np.log(p)\n",
    "sigma = np.sqrt(pi_entr @ np.square(np.log(p) - mu))\n",
    "print (\"mu: %.3f, sigma: %.3f\" % (mu, sigma))\n",
    "plt.close()\n",
    "plt.plot(p, pilognormal, c='black')\n",
    "plt.savefig(\"../price_distr_derivs/figs/lognormal.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "annualized_volatility = np.std(np.prod(np.exp(np.random.normal(loc=mu, scale=sigma, size=(1_000_000, int(365/36)))) / price, axis=1))\n",
    "annualized_volatility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lower and upper bounds on the mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "result = cp.Problem(cp.Maximize(p @ pi), cons).solve(solver=cp.MOSEK)\n",
    "print (\"lower bound:\", result)\n",
    "result = cp.Problem(cp.Minimize(p @ pi), cons).solve(solver=cp.MOSEK)\n",
    "print (\"upper bound:\", result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Upper/lower bounds on CDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "upper_bounds = []\n",
    "lower_bounds = []\n",
    "for i in np.arange(0, m, 50):\n",
    "    try:\n",
    "        upper_bound = cp.Problem(cp.Maximize(cp.sum(pi[:i])), cons).solve(solver=cp.MOSEK)\n",
    "    except:\n",
    "        upper_bound = np.nan\n",
    "    upper_bounds.append(upper_bound)\n",
    "    try:\n",
    "        lower_bound = cp.Problem(cp.Minimize(cp.sum(pi[:i])), cons).solve(solver=cp.MOSEK)\n",
    "    except:\n",
    "        lower_bound = np.nan\n",
    "    lower_bounds.append(lower_bound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(p[::50], upper_bounds, c='black', label='upper')\n",
    "plt.plot(p[::50], lower_bounds, '--', c='black', label='lower')\n",
    "plt.axvline(price, c='gray', label='price')\n",
    "plt.xlim(0,20000)\n",
    "plt.ylim(.0001)\n",
    "plt.ylim(5e-3)\n",
    "plt.legend()\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('$\\mathbf{Prob}(p \\leq x)$')\n",
    "plt.savefig('../price_distr_derivs/figs/cdf_bounds.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Upper/lower bounds on VaR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "plt.semilogx(lower_bounds, (p[::10] - price) / price, c='black', label='upper')\n",
    "plt.semilogx(upper_bounds, (p[::10] - price) / price, '--', c='black', label='lower')\n",
    "plt.legend()\n",
    "plt.xlim(.00001, 1)\n",
    "plt.xlabel('$\\epsilon$')\n",
    "plt.ylabel('$\\mathbf{VaR}(p, \\epsilon)$')\n",
    "plt.savefig('../price_distr_derivs/figs/var_bounds.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Upper bound on CVaR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvar_upper = []\n",
    "epss = np.logspace(-4,-.1,25)\n",
    "for eps in epss:\n",
    "    t = cp.Variable(1)\n",
    "    result = cp.Problem(cp.Maximize(t), cons + [p[i] + pi @ np.maximum(p - p[i], 0) / (1 - eps) >= t for i in range(m)]).solve(solver=cp.MOSEK)\n",
    "    cvar_upper.append(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogx(np.logspace(-4,-.1,25), (np.array(cvar_upper) - price) / price)\n",
    "plt.semilogx(upper_bounds, (p[::10] - price) / price, '--', c='black', label='lower')\n",
    "plt.xlim(1e-4,10**(-.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Upper/lower bounds on complementary CDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_bounds = []\n",
    "lower_bounds = []\n",
    "for i in np.arange(0, m, 10):\n",
    "    try:\n",
    "        upper_bound = cp.Problem(cp.Maximize(cp.sum(pi[i:])), cons).solve(solver=cp.MOSEK)\n",
    "    except:\n",
    "        upper_bound = np.nan\n",
    "    upper_bounds.append(upper_bound)\n",
    "    try:\n",
    "        lower_bound = cp.Problem(cp.Minimize(cp.sum(pi[i:])), cons).solve(solver=cp.MOSEK)\n",
    "    except:\n",
    "        lower_bound = np.nan\n",
    "    lower_bounds.append(lower_bound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "plt.semilogy((p[::10] - price) / price, upper_bounds, c='black', label='upper')\n",
    "plt.semilogy((p[::10] - price) / price, lower_bounds, '--', c='black', label='lower')\n",
    "plt.axvline(0, c='gray', label='price')\n",
    "plt.xlim(-.2, .2)\n",
    "plt.ylim(.0001)\n",
    "plt.legend()\n",
    "plt.xlabel('a')\n",
    "plt.ylabel('prob(return >= a)')\n",
    "plt.savefig('../price_distr_derivs/figs/ccdf_bounds.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogx(upper_bounds, (p[::10] - price) / price, c='black', label='upper')\n",
    "plt.semilogx(lower_bounds, (p[::10] - price) / price, '--', c='black', label='lower')\n",
    "plt.axvline(0, c='blue', label='price')\n",
    "plt.xlim(.0001, 1)\n",
    "plt.legend()\n",
    "plt.xlabel('$\\epsilon$')\n",
    "plt.ylabel('$\\mathbf{VaR}(-p,\\epsilon)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bounds on costs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ask_lower_bounds = []\n",
    "bid_upper_bounds = []\n",
    "for i in range(puts_ask.size):\n",
    "    idx = np.ones(n, dtype=bool)\n",
    "    idx[i] = False\n",
    "    idx[puts_ask.size + i] = False\n",
    "    pi = cp.Variable(m)\n",
    "    cbuy = cp.Variable()\n",
    "    csell = cp.Variable()\n",
    "    \n",
    "    cons = [pi >= 0, cp.sum(pi) == 1, P[:,idx].T @ pi <= d[idx], csell <= cbuy,\n",
    "            P[:,~idx].T @ pi <= cp.vstack([cbuy, -csell])[:,0]]\n",
    "    ask_lower_bound = cp.Problem(cp.Minimize(cbuy), cons).solve(solver=cp.MOSEK)\n",
    "    bid_upper_bound = cp.Problem(cp.Maximize(csell), cons).solve(solver=cp.MOSEK)\n",
    "    \n",
    "    ask_lower_bounds.append(ask_lower_bound)\n",
    "    bid_upper_bounds.append(bid_upper_bound)\n",
    "    assert puts_ask[i] >= ask_lower_bound\n",
    "    assert puts_bid[i] <= bid_upper_bound\n",
    "    print (\"Put\", i)\n",
    "    print (\"ask %3.3f bid %3.3f\" % (puts_ask[i], puts_bid[i]))\n",
    "    print (\"ask lower bound %3.3f\" % (ask_lower_bound))\n",
    "    print (\"bid upper bound %3.3f\" % (bid_upper_bound))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "plt.plot(puts_strike[idx], (puts_ask - np.array(ask_lower_bounds))[idx], c='k')\n",
    "plt.xlabel('strike price')\n",
    "plt.ylabel('distance to lower bound')\n",
    "plt.savefig('../price_distr_derivs/figs/ask_bounds.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = np.argsort(puts_strike)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(puts_strike[idx], (np.array(bid_upper_bounds) - puts_bid)[idx], c='k')\n",
    "plt.xlabel('strike price')\n",
    "plt.ylabel('distance to upper bound')\n",
    "plt.savefig('../price_distr_derivs/figs/bid_bounds.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bounds on the price of a new investment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a binary option that pays 1 if price stays within +/- 2%, and 0 otherwise\n",
    "pi = cp.Variable(m)\n",
    "c = cp.Variable()\n",
    "\n",
    "cons = [pi >= 0, cp.sum(pi) == 1, P[:,idx].T @ pi <= d[idx], (np.abs(p - price) / price < .02).astype(np.float) @ pi == c]\n",
    "lower = cp.Problem(cp.Minimize(c), cons).solve(solver=cp.MOSEK)\n",
    "upper = cp.Problem(cp.Maximize(c), cons).solve(solver=cp.MOSEK)\n",
    "lower, upper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}