Add new tutorial for multi-currencies arbitrage using Bellman-Ford Algorithm

Documentation/Python/Bellman-Ford Currency Arbitrage.ipynb

@@ -0,0 +1,311 @@
+{
+  "nbformat": 4,
+  "nbformat_minor": 0,
+  "metadata": {
+    "colab": {
+      "provenance": [],
+      "toc_visible": true
+    },
+    "kernelspec": {
+      "name": "python3",
+      "display_name": "Python 3"
+    },
+    "language_info": {
+      "name": "python"
+    }
+  },
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "source": [
+        "![QuantConnect Logo](https://cdn.quantconnect.com/web/i/icon.png)\n",
+        "<hr>"
+      ],
+      "metadata": {
+        "id": "j0a1UHtbMwGX"
+      }
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "# Crypto-Arbitrage with Bellman-Ford Algorithm written in Python\n",
+        "## Definitions\n",
+        "> Triangular arbitrage is the result of a discrepancy between three foreign currencies that occurs when the currency's exchange rates do not exactly match up. These opportunities are rare and traders who take advantage of them usually have advanced computer equipment and/or programs to automate the process.[^1]\n",
+        "\n",
+        "###### Graph:\n",
+        "> A graph is a combinatorial object composed of a set of vertices V (also known as nodes) and a set of edges E. The edges correspond to pairs of vertices, which are generally distinct, and without a notion of order in the sense where (u,v) and (v,u) denote the same edge.\n",
+        "At times, we consider a variant, the directed graph, where the edges have an ori- entation. In this case, the edges are usually known as arcs. The arc (u,v) has origin u and destination v. Most of the algorithms described in this book operate on directed graphs but can be applied to non-directed graphs by replacing each edge (u,v) by two arcs (u,v) and (v,u).\n",
+        "Graphs can contain additional information, such as weights or letters, in the form of labels on the vertices or the edges.\n",
+        "\n",
+        "###### Bellman-Ford algorithm:\n",
+        "> The Bellman-Ford algorithm finds the minimum weight path from a single source vertex to all other vertices on a weighted directed graph.\n",
+        "\n",
+        "Our goal is to develop a systematic method for detecting arbitrage opportunities by framing the problem in the language of graphs. \n",
+        "\n",
+        "## Approach\n",
+        "We will assign currencies to different vertices, and let the edge weight represent the exchange rate.\n",
+        "Find a cycle in the graph such that multiplying all the edge weights along that cycle results in a value greater than 1. In fact we have already described an algorithm that can find such a path – the problem is equivalent to finding a negative-weight cycle, provided we do some preprocessing on the edges.\n",
+        "\n",
+        "We note that Bellman-Ford computes the path weight by adding the individual edge weights. To make this work for exchange rates, which are multiplicative, a fix is to first take the logs of all the edge weights. Thus when we sum edge weights along a path we are actually multiplying exchange rates – we can recover the multiplied quantity by exponentiating the sum. Secondly, Bellman-Ford attempts to find minimum weight paths and negative edge cycles, whereas our arbitrage problem is about maximising the amoun t of currency received. Thus as a simple modification, we must also make our log weights negative.\n",
+        "Now we are able to apply Bellman-Ford. The minimal weight between two currency vertices corresponds to the optimal exchange rate, the value of which can be found by by exponentiating the negative sum of weights along the path. A negative-weight cycle on the negative-log graph corresponds to an arbitrage opportunity.\n",
+        "\n",
+        "## Code\n",
+        "> List of functions:\n",
+        "```\n",
+        "get_price()\n",
+        "```\n",
+        "- get last price from *QuantConnect[^3] API* and put together into a dataframe by using pandas\n",
+        "```\n",
+        "Trading.strategy()\n",
+        "```\n",
+        "- recall the `Graph` class[^2] and the `Graph.bellman_ford()` to perform the strategy and print the  *boolean* variable `bol` (`True` if negative cycles were detected, `False` otherwise) and the **profit**  expressed as %\n",
+        "\n",
+        "\n",
+        "\n",
+        "[^1]:[Investopedia: Triangular Arbitrage](https://www.investopedia.com/terms/t/triangulararbitrage.asp)\\\n",
+        "[^2]: The code to implement it has been taken from this [book](https://amzn.to/3bBI8tP)\\\n",
+        "[^3]: [Datasets from QuantConnect]https://www.quantconnect.com/docs/v2/research-environment/datasets/crypto\n"
+      ],
+      "metadata": {
+        "id": "A2I-7YxuM2y6"
+      }
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "##### Importing modules\t"
+      ],
+      "metadata": {
+        "id": "FVW7GG_IM98g"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "import datetime as dt\n",
+        "import numpy as np\n",
+        "import pandas as pd\n",
+        "import warnings\n",
+        "warnings.filterwarnings(\"ignore\")\n",
+        "\n",
+        "\n",
+        "# QuantBook Analysis Tool \n",
+        "# For more information see [https://www.quantconnect.com/docs/v2/our-platform/research/getting-started]\n",
+        "qb = QuantBook()\n"
+      ],
+      "metadata": {
+        "id": "39qjX9M2M3Mk"
+      },
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "##### Defining function for getting data\t"
+      ],
+      "metadata": {
+        "id": "vduiwBO_NEfW"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "\n",
+        "    \n",
+        "def get_data(sym1:str,sym2:str, sym3:str , start_time: datetime, end_time: datetime):\n",
+        "    \"\"\"\n",
+        "    @sym1: ticker name\n",
+        "    @sym2: ticker name\n",
+        "    @sym3: ticker name\n",
+        "    @start_time: date to start getting data\n",
+        "    @end_time: last day of getting data\n",
+        "    \"\"\"\n",
+        "\n",
+        "    ticker1 = qb.AddCrypto(sym1).Symbol\n",
+        "    ticker2 = qb.AddCrypto(sym2).Symbol\n",
+        "    ticker3 = qb.AddCrypto(sym3).Symbol\n",
+        "\n",
+        "    price1 = qb.History(ticker1, start_time, end_time, Resolution.Tick)['lastprice']\n",
+        "    price2 = qb.History(ticker2, start_time, end_time, Resolution.Tick)['lastprice']\n",
+        "    price3 = qb.History(ticker3, start_time, end_time, Resolution.Tick)['lastprice']\n",
+        "\n",
+        "    df = pd.DataFrame({'symbol':list(),'price':list()})\n",
+        "    price1 = -np.log(price1.values.tolist())\n",
+        "    price2 = -np.log(price2.values.tolist())\n",
+        "    price3 = -np.log(price3.values.tolist())\n",
+        "\n",
+        "    df=df.append([{'symbol': sym1 + '_' + str(price1.index()[-1]) , 'price': price1[-1]}], ignore_index=True)\n",
+        "    df=df.append([{'symbol': sym2 + '_' + str(price2['Date'][-1]), 'price': price2[-1]}], ignore_index=True)\n",
+        "    df=df.append([{'symbol': sym3 + '_' + str(price3['Date'][-1]) , 'price': price3[-1]}], ignore_index=True)\n",
+        "    return df\n"
+      ],
+      "metadata": {
+        "id": "3qWmjQHjNGsu"
+      },
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "##### Create the Graph class\t"
+      ],
+      "metadata": {
+        "id": "yGDNr1NhNUvB"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "\n",
+        "class Graph:\n",
+        "    def __init__(self):\n",
+        "        self.neighbors = []\n",
+        "        self.name2node = {}\n",
+        "        self.node2name = []\n",
+        "        self.weight = []\n",
+        "    \n",
+        "    def __len__(self):\n",
+        "        return len(self.node2name)\n",
+        "    def __getitem__(self,v):\n",
+        "        return self.neighbors[v]\n",
+        "    \n",
+        "    def add_node(self,name):\n",
+        "        assert name not in self.name2node\n",
+        "        self.name2node[name] = len(self.name2node)\n",
+        "        self.node2name.append(name)\n",
+        "        self.neighbors.append([]) \n",
+        "        self.weight.append({})\n",
+        "        return self.name2node[name]\n",
+        "    \n",
+        "    def add_edge(self,name_u,name_v,weight_uv=None):\n",
+        "        self.add_arc(name_u, name_v, weight_uv) \n",
+        "        self.add_arc(name_v, name_u, weight_uv)\n",
+        "\n",
+        "    def add_arc(self,name_u,name_v,weight_uv=None):\n",
+        "        u = self.name2node[name_u]\n",
+        "        v = self.name2node[name_v] \n",
+        "        self.neighbors[u].append(v)\n",
+        "        self.weight[u][v] = weight_uv\n",
+        "\n",
+        "    def bellman_ford(self, weight, source=0):\n",
+        "        graph = self\n",
+        "        n = len(graph)\n",
+        "        dist = [float('inf')] * n\n",
+        "        prec = [None]*n\n",
+        "        dist[source] = 0\n",
+        "        for nb_iterations in range(n):\n",
+        "            changed = False\n",
+        "            for node in range(n):\n",
+        "                for neighbor in graph[node]:\n",
+        "                    alt = dist[node] + weight[node][neighbor]\n",
+        "                    if alt < dist[neighbor]:\n",
+        "                        dist[neighbor] = alt\n",
+        "                        prec[neighbor] = node\n",
+        "                        changed = True\n",
+        "                if not changed:\n",
+        "                    return dist,prec,False\n",
+        "        return dist, prec, True\n",
+        "\n",
+        "\n"
+      ],
+      "metadata": {
+        "id": "DSPN5zyNNX4j"
+      },
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "##### Defining the strategy\t"
+      ],
+      "metadata": {
+        "id": "HHHCYHp4Nb1D"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "def strategy(df: pd.DataFrame):\n",
+        "        \"\"\"\n",
+        "        df: DataFrame of prices\n",
+        "        \"\"\"\n",
+        "        global g\n",
+        "        g = Graph()\n",
+        "        for i in df.symbol:\n",
+        "            g.add_node(i)\n",
+        "        for j in df.price:\n",
+        "            g.weight.append((j))\n",
+        "        \n",
+        "        for m in range(len(df)-1):\n",
+        "            g.add_arc(df.symbol[m],df.symbol[m+1],df.price[m])\n",
+        "        for n in reversed(range(len(df)-1)):\n",
+        "            g.add_arc(df.symbol[n],df.symbol[n+1],df.price[n])\n",
+        "        \n",
+        "        dist, prec, bol = g.bellman_ford(g.weight,source=0)\n",
+        "        #####\n",
+        "        tot = 0\n",
+        "        for i in dist:\n",
+        "            tot *= i\n",
+        "        profit = np.exp(-tot)-1\n",
+        "        if bol:\n",
+        "            print(f\"Profit from the strategy is: {profit*100:.2g}%\\n\")\n",
+        "        return bol, profit\n",
+        "\n"
+      ],
+      "metadata": {
+        "id": "GvU_aew-NcWu"
+      },
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "##### implementing the strategy\t"
+      ],
+      "metadata": {
+        "id": "2eoGzIJONgPX"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "import time \n",
+        "start_time = datetime.datetime(2023, 5, 1)\n",
+        "end_time = datetime.datetime(2023, 5, 3)\n",
+        "for i in range(100):\n",
+        "    df = get_data(\"BTCUSD\",\"ETHUSD\",\"LTCUSD\",start_time,end_time)\n",
+        "    print(df)\n",
+        "    print(strategy(df))\n",
+        "    time.sleep(5)"
+      ],
+      "metadata": {
+        "id": "0tymEMjwNgpw"
+      },
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "##### Final Consideration:\n",
+        "This algorith was built for only three cryptocurrencies, but you can modify the code to apply with real time data and most important: **multi currencies**. You can actually add as many tickers as you want to test this algorithm (you need a proper subscription) and if you want to trade in a seconds environment, remember to change the code accordingly \t"
+      ],
+      "metadata": {
+        "id": "XJThPMMRNmia"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [],
+      "metadata": {
+        "id": "Rt0ADfeKNm4y"
+      },
+      "execution_count": null,
+      "outputs": []
+    }
+  ]
+}

Documentation/Python/capital_structure_arbitrage.ipynb

@@ -0,0 +1,496 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Capital Structure Arbitrage"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Capital structure arbitrage is the term used to describe the fashion for arbitraging equity claims against fixed income and convertible claims. At its most sophisticated, practitioners build elaborate models of the capital structure of a company to determine the relative values of the various claims—in particular, stock, bonds, and convertible bonds. At its simplest, the trader looks to see if equity puts are cheaper than credit derivatives and if so buys the one and sells the other"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Denoting the value of a risk-free put, call and bond by $P_0$, $C_0$, and $B_0$ and the value of risky claims on the issuer (I) of the stock by $P_I$, $C_I$, and $B_I$, we obtain"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "the risk-free put is worth more than the risky put. The excess value is equal to the difference in risky and risk-free bond prices (times the strike price). With maturity-independent rates and credit spreads for clarity and setting t = 0, we obtain\n",
+    "\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$B_0-B_I=e^{-r T}\\left(1-e^{-\\lambda T}\\right)$\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sympy as smp"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle B_{0} - B_{1}$"
+      ],
+      "text/plain": [
+       "B_0 - B_1"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "B_0,B_1,r,T,lambda_ = smp.symbols('B_0 B_1 r T lambda_', real=True) \n",
+    "expr = B_0 - B_1 \n",
+    "expr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle B_{0} - B_{1} = \\left(1 - e^{- T \\lambda_{}}\\right) e^{- T r}$"
+      ],
+      "text/plain": [
+       "Eq(B_0 - B_1, (1 - exp(-T*lambda_))*exp(-T*r))"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    " #for writing equation using smp.Ex()\n",
+    "smp.Eq(expr,smp.exp(-r*T)*(1-smp.exp(-lambda_*T)))\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "the extra value is the strike price times the (pseudo-) probability that default occurs. This payoff is also more or less exactly the payoff of a default put in the credit derivatives market.\\\n",
+    "the trader buys an equity option on the exchange at a ‘‘very high’’ implied volatility and sells a default put on the same stock in the credit derivatives market locking in a risk-free return.\\\n",
+    " buy one at-the-money put and sell two puts struck at half the current stock price"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\\begin{aligned}\n",
+    "&\\begin{aligned}\n",
+    "\\sigma_{l o c}^2(K, T, S) & =\\sigma^2-\\lambda \\frac{K \\frac{\\partial C}{\\partial K}}{\\frac{1}{2} K^2 \\frac{\\partial^2 C}{\\partial K^2}} \\\\\n",
+    "& =\\sigma^2+2 \\lambda \\sigma \\sqrt{T} \\frac{N\\left(d_2\\right)}{N^{\\prime}\\left(d_2\\right)}\n",
+    "\\end{aligned}\\\\\n",
+    "&d_2=\\frac{\\log S / K+\\lambda T}{\\sigma \\sqrt{T}}-\\frac{\\sigma \\sqrt{T}}{2}\n",
+    "\\end{aligned}"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Then, for very low strikes,"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$\\sigma_{l o c}^2(K, T, S) \\approx \\sigma^2+2 \\lambda \\sigma \\sqrt{T} \\sqrt{2 \\pi} e^{+d_2^2 / 2}$"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It follows that implied volatility in the jump-to-ruin model increases very fast as the strike decreases from at-the-money and tends to the constant σ for high strikes"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sigma_loc_squared = smp.symbols('sigma_loc_squared', cls=smp.Function) \n",
+    "sigma,K,d_2, S, f_d_2, mu_stock, sigma_stock = smp.symbols('sigma K d_2 S f_d_2 mu_stock sigma_stock', real=True) "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{2 \\sqrt{T} \\lambda_{} \\sigma \\left(\\begin{cases} \\frac{\\operatorname{erf}{\\left(\\frac{\\sqrt{2} \\left(T - \\mu_{stock}\\right)}{2 \\sigma_{stock}} \\right)}}{2} + \\frac{1}{2} & \\text{for}\\: \\left(\\left|{\\arg{\\left(\\sigma_{stock} \\right)}}\\right| \\leq \\frac{\\pi}{4} \\wedge \\left|{2 \\arg{\\left(T \\right)} - 4 \\arg{\\left(\\sigma_{stock} \\right)} + 2 \\arg{\\left(\\frac{T - \\mu_{stock}}{T} \\right)} + 2 \\pi}\\right| < \\pi\\right) \\vee \\left|{\\arg{\\left(\\sigma_{stock} \\right)}}\\right| < \\frac{\\pi}{4} \\\\\\frac{\\sqrt{2} \\int\\limits_{-\\infty}^{0} e^{- \\frac{\\left(z + T - \\mu_{stock}\\right)^{2}}{2 \\sigma_{stock}^{2}}}\\, dz}{2 \\sqrt{\\pi} \\sigma_{stock}} & \\text{otherwise} \\end{cases}\\right)}{d_{2}} + \\sigma^{2}$"
+      ],
+      "text/plain": [
+       "2*sqrt(T)*lambda_*sigma*Piecewise((erf(sqrt(2)*(T - mu_stock)/(2*sigma_stock))/2 + 1/2, (Abs(arg(sigma_stock)) < pi/4) | ((Abs(arg(sigma_stock)) <= pi/4) & (Abs(2*arg(T) - 4*arg(sigma_stock) + 2*arg((T - mu_stock)/T) + 2*pi) < pi))), (sqrt(2)*Integral(exp(-(_z + T - mu_stock)**2/(2*sigma_stock**2)), (_z, -oo, 0))/(2*sqrt(pi)*sigma_stock), True))/d_2 + sigma**2"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "from sympy.stats import P, E, variance, Die, Normal\n",
+    "from sympy import simplify\n",
+    "f_d_2 = Normal('d_2',mu_stock, sigma_stock) \n",
+    "F_d_2 = simplify(P(f_d_2 <= T))\n",
+    "d_2 = (smp.log(S/K) + lambda_ * T) / (sigma * smp.sqrt(T)) - ((sigma * smp.sqrt(T)) / 2)\n",
+    "sigma_loc_squared = sigma**2 + 2 * lambda_ * sigma * smp.sqrt(T) * (F_d_2 / f_d_2)\n",
+    "sigma_loc_squared"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "fair price of a zero coupon bond of GT (assuming zero rates) should be given by"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$P_t=e^{-\\lambda t} R+\\left(1-e^{-\\lambda t}\\right)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "P_t = smp.symbols('P_t', cls=smp.Function) \n",
+    "t,R = smp.symbols('t R', real=True) \n",
+    "P_t = smp.exp(-lambda_*t)*R+(1-smp.exp(-lambda_*t))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where R is the recovery rate\\\n",
+    "most of the volatility skew for stocks with high credit spreads can be ascribed to default risk."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Model setup"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The level of V at which the company defaults is given by L D where D is\n",
+    "today’s value of its debt (per share) and L is the recovery rate. As discussed above, it is further assumed that the recovery rate L is a lognormally distributed random variable with mean $\\bar{L}$ and standard deviation λ so that"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$L D=\\bar{L} D e^{\\lambda Z-\\lambda^2 / 2}$\\\n",
+    "where Z ∼ N(0, 1). The random variable Z is assumed to be independent of $W_t$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle D L = D barL e^{- \\frac{\\lambda_{dev}^{2}}{2} + \\frac{\\lambda_{dev} Z}{2}}$"
+      ],
+      "text/plain": [
+       "Eq(D*L, D*barL*exp(-lambda_dev**2/2 + lambda_dev*Z/2))"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "L,D,barL,lambda_dev = smp.symbols('L D barL lambda_dev', real=True) \n",
+    "Z = Normal('Z',0, 1) \n",
+    "exprr = L*D\n",
+    "smp.Eq(exprr,barL*D*smp.exp((lambda_dev*Z-lambda_dev**2)/2))\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Survival Probability"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$X_t:=\\sigma W_t-\\lambda Z-\\frac{\\sigma^2 t}{2}-\\frac{\\lambda^2}{2}$\\\n",
+    "\\\n",
+    "Then $X_t$ is normally distributed with\\\n",
+    "\\\n",
+    "$\\begin{aligned} \\mathbb{E}\\left[X_t\\right] & =-\\frac{\\sigma^2}{2}\\left(t+\\frac{\\lambda^2}{\\sigma^2}\\right) \\\\ \\operatorname{Var}\\left[X_t\\right] & =\\sigma^2\\left(t+\\frac{\\lambda^2}{\\sigma^2}\\right)\\end{aligned}$\\\n",
+    "Default occurs when\\\n",
+    "\\\n",
+    "$V=V_0 e^{\\sigma W_t-\\sigma^2 t / 2}=L D=\\bar{L} D e^{\\lambda Z-\\lambda^2 / 2}$\\\n",
+    "\\\n",
+    "or equivalently when\\\n",
+    "\\\n",
+    "$X_t=\\log \\left(\\frac{\\bar{L} D}{V_0}\\right)-\\lambda^2$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle V_{0} e^{W_{t} \\sigma - \\frac{\\sigma^{2} t}{2}}$"
+      ],
+      "text/plain": [
+       "V_0*exp(W_t*sigma - sigma**2*t/2)"
+      ]
+     },
+     "execution_count": 19,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "X_t,sigma,W_t,E_X_t,Var_X_t,V,V_0 = smp.symbols('X_t sigma W_t E_X_t Var_X_t V V_0', real=True) \n",
+    "X_t = sigma*W_t-lambda_dev*Z-((sigma**2*t)/2) - (lambda_dev**2/2)\n",
+    "E_X_t = -(sigma**2/2)*(t+(lambda_dev**2)/2)\n",
+    "Var_X_t = sigma**2*(t+(lambda_dev**2)/2)\n",
+    "V = V_0 * smp.exp(sigma*W_t-sigma**2*t/2)\n",
+    "\n",
+    "#Default case\n",
+    "if V == exprr or X_t == smp.log((barL*D)/V_0)-lambda_dev**2:\n",
+    "    print(\"Default\")\n",
+    "V\n",
+    "\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Since $\\hat{X}$ is a Brownian motion with drift, the probability of survival (or the probability of not hitting the default barrier) is given by the Black-Scholes- like formula\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$P_t=N\\left(-\\frac{A_t}{2}+\\frac{\\log d}{A_t}\\right)-d N\\left(-\\frac{A_t}{2}-\\frac{\\log d}{A_t}\\right)$\\\n",
+    "\\\n",
+    "with\\\n",
+    "\\\n",
+    "$d=\\frac{V_0 e^{\\lambda^2}}{\\bar{L} D} ; A_t^2=\\sigma^2 t+\\lambda^2$\\\n",
+    "\\\n",
+    "Since $P_t$ is the probability of survival up to time t, it may be estimated directly from the prices of risky instruments such as bonds and credit default swaps (CDS).\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 37,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle expr_{1} - \\frac{V_{0} e^{\\lambda_{dev}^{2}} expr_{2}}{D L}$"
+      ],
+      "text/plain": [
+       "expr1 - V_0*exp(lambda_dev**2)*expr2/(D*L)"
+      ]
+     },
+     "execution_count": 37,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A_t_squared,d,A_t = smp.symbols('A_t_squared d A_t', real=True) \n",
+    "P_t = smp.symbols('P_t', cls=smp.Function)\n",
+    "exprs1 = Normal('expr1',0,1 )\n",
+    "exprs2 = Normal('expr2',0,1 )\n",
+    "expr1 = (-A_t/2) + (smp.log(d)/A_t) \n",
+    "expr2 = (-A_t/2) - (smp.log(d)/A_t) \n",
+    "d = V_0 * smp.exp(lambda_dev**2) / (L*D) \n",
+    "A_t_squared = sigma**2 * t + lambda_dev**2\n",
+    "P_t = exprs1 - d * exprs2\n",
+    "P_t"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Equity volatility"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The stock price S is approximately related (neglecting the time value of the option) to the firm value V via\n",
+    "$V \\approx L D+S$ then\\\n",
+    "\\\n",
+    "$\\sigma \\sim \\frac{\\delta V}{V} \\approx \\frac{\\delta S}{S+L D}=\\frac{\\delta S}{S} \\frac{S}{S+L D} \\sim \\sigma_S \\frac{S}{S+L D}$\\\n",
+    "\\\n",
+    "where $σ_S$ is the stock volatility. We see that as the stock price rises, keeping σ fixed, the volatility $σ_S$ of the stock declines. Conversely, as the stock price S approaches zero, the stock volatility increases as 1/S."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Model calibration"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We end up with the following model in terms of market observables\\\n",
+    "\\\n",
+    "$P_t=N\\left(-\\frac{A_t}{2}+\\frac{\\log d}{A_t}\\right)-d N\\left(-\\frac{A_t}{2}-\\frac{\\log d}{A_t}\\right)$\\\n",
+    "\\\n",
+    "$d=\\frac{S_0+\\bar{L} D}{\\bar{L} D} e^{\\lambda^2} ; A_t^2=\\left(\\sigma_S^* \\frac{S^*}{S^*+\\bar{L} D}\\right)^2 t+\\lambda^2$\\\n",
+    "\\\n",
+    "where $S^*$ is some reference stock price and $σ_{S}^*$ the stock volatility at that price."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "S_0,S_t = smp.symbols('S_0 S_t', real=True)\n",
+    "expr1 = (-A_t/2) + (smp.log(d)/A_t) \n",
+    "expr2 = (-A_t/2) - (smp.log(d)/A_t) \n",
+    "d = ((S_0 + barL*D ) / (L*D)) *smp.exp(lambda_dev**2)\n",
+    "A_t_squared = (sigma_stock**2 * (S_t/(S_t + L*D))**2)*t + lambda_dev**2\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Getting $\\bar{L}$, λ and D from company and industry data rather than from the term structure of credit spreads would theoretically enable us to identify rich and cheap claims."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# To evaluate the equation use the expression: \n",
+    "A_t_squared.subs((sigma_stock,sigma_stock),(S_t,S_t),(L,L),(D,D),(t,t),(lambda_dev,lambda_dev))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> all the text and formulas are taken from \"The Volatility Surface: A Practitioner's Guide\" by Jim Gatheral"
+   ]
+  }
+ ],
+ "metadata": {
+  "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.10.2"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

Local_Stochastic_Volatility_Model.ipynb

@@ -0,0 +1,266 @@
+{
+  "nbformat": 4,
+  "nbformat_minor": 0,
+  "metadata": {
+    "colab": {
+      "provenance": []
+    },
+    "kernelspec": {
+      "name": "python3",
+      "display_name": "Python 3"
+    },
+    "language_info": {
+      "name": "python"
+    }
+  },
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "source": [
+        "![QuantConnect Logo](https://cdn.quantconnect.com/web/i/icon.png)\n",
+        "<hr>"
+      ],
+      "metadata": {
+        "id": "6iZIAgxW8Ari"
+      }
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "# Dupire model\n",
+        "$C\\left(S_0, K, T\\right)=\\int_K^{\\infty} d S_T \\varphi\\left(S_T, T ; S_0\\right)\\left(S_T-K\\right)$\n",
+        "### risk-neutral density function φ of the final spot $S_T$\n",
+        "$\\varphi\\left(K, T ; S_0\\right)=\\frac{\\partial^2 C}{\\partial K^2}$\n",
+        "## Drift\n",
+        "$\\mu = r_t − D_t$\n",
+        "## local volatility\n",
+        "$\\sigma^2\\left(K, T, S_0\\right)=\\frac{\\frac{\\partial C}{\\partial T}}{\\frac{1}{2} K^2 \\frac{\\partial^2 C}{\\partial K^2}}$\n",
+        "## undiscounted option price C in terms of the strike price K:\n",
+        "### which is the Dupire equation when the underlying stock has risk-neutral drift μ.\n",
+        "$\\frac{\\partial C}{\\partial T}=\\frac{\\sigma^2 K^2}{2} \\frac{\\partial^2 C}{\\partial K^2}+\\left(r_t-D_t\\right)\\left(C-K \\frac{\\partial C}{\\partial K}\\right)$\n",
+        "# BSM model\n",
+        "$Call=S  N\\left(d_1\\right)-K e^{-r \\tau} N\\left(d_2\\right)$ \n",
+        "\\\n",
+        "$d_1=\\frac{\\ln (S / K)+\\left(r-y+\\sigma^2 / 2\\right) \\tau}{\\sigma \\sqrt{\\tau}}$\n",
+        "\\\n",
+        "$d_2=\\frac{\\ln (S / K)+\\left(r-y-\\sigma^2 / 2\\right) \\tau}{\\sigma \\sqrt{\\tau}}=d_1-\\sigma \\sqrt{\\tau}$\n",
+        "\n",
+        "## PDF of Normal  \n",
+        "$f(x)=\\frac{1}{\\sigma \\sqrt{2 \\pi}} e^{-\\frac{1}{2}\\left(\\frac{x-\\mu}{\\sigma}\\right)^2}$\n",
+        "\n",
+        "## CDF of Normal\n",
+        "$N(x)=\\frac{1}{\\sqrt{2 \\pi}} \\int_{-\\infty}^x e^{-z^2 / 2} d z$\n"
+      ],
+      "metadata": {
+        "id": "fqocn6D58JDW"
+      }
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "# Importing Libraries"
+      ],
+      "metadata": {
+        "id": "LeIimwLE8FOu"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "import sympy as smp \n",
+        "import numpy as np\n",
+        "from sympy.stats import P, E, variance, Die, Normal\n",
+        "from sympy import simplify\n"
+      ],
+      "metadata": {
+        "id": "dk37ioq7wpBX"
+      },
+      "execution_count": 2,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "\n",
+        "# Risk-neutral density function φ of the final spot $S_T$\n",
+        "$\\varphi\\left(K, T ; S_0\\right)=\\frac{\\partial^2 C}{\\partial K^2}$"
+      ],
+      "metadata": {
+        "id": "BOs4Ha4_9Go_"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "#define the function symbols\n",
+        "Call,varphi,N = smp.symbols('Call varphi N', cls=smp.Function) \n",
+        "S_0,K,T,S_T,D,r,sigma,tau ,d_1,d_2 = smp.symbols('S_0 K T S_T D r sigma tau d_1 d_2',real=True) \n",
+        "\n",
+        "\n",
+        "Call = smp.symbols('Call', cls=smp.Function) \n",
+        "f_d_1,f_d_2,F_d_1,F_d_2 = smp.symbols('f_d_1 f_d_2 F_d_1 F_d_2' , cls=smp.Function) \n",
+        "mu_google,sigma_google= smp.symbols('mu_google sigma_google',real=True) \n",
+        "f_d_1 = Normal('d_1', mu_google, sigma_google)\n",
+        "f_d_2 = Normal('d_2', mu_google, sigma_google)\n",
+        "\n",
+        "\n",
+        "F_d_1 = simplify(P(f_d_1>tau))\n",
+        "F_d_2 = simplify(P(f_d_2>tau))\n",
+        "Call = S_T  * F_d_1 - K * smp.exp(-r*tau) * F_d_2\n",
+        "\n",
+        "\n",
+        "varphi = smp.diff(Call,K,2)\n",
+        "\n",
+        "#varphi.subs([(Call,Call),(K,int(K))]).evalf()"
+      ],
+      "metadata": {
+        "id": "0G1VQsJF8M0M"
+      },
+      "execution_count": 4,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "\n",
+        "## Undiscounted option price C in terms of the strike price K:\n",
+        "### which is the Dupire equation when the underlying stock has risk-neutral drift μ.\n",
+        "$\\frac{\\partial C}{\\partial T}=\\frac{\\sigma^2 K^2}{2} \\frac{\\partial^2 C}{\\partial K^2}+\\left(r_t-D_t\\right)\\left(C-K \\frac{\\partial C}{\\partial K}\\right)$"
+      ],
+      "metadata": {
+        "id": "ORnbFO8A8XAJ"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "dCdT = sigma**2 * K / 2 * varphi + (r - D)* (Call - K * smp.diff(Call,K))\n",
+        "dCdT"
+      ],
+      "metadata": {
+        "colab": {
+          "base_uri": "https://localhost:8080/",
+          "height": 137
+        },
+        "id": "Xf5AbHNS8XT6",
+        "outputId": "fec23eba-66b6-4c57-8588-38d7d839f9ff"
+      },
+      "execution_count": 5,
+      "outputs": [
+        {
+          "output_type": "execute_result",
+          "data": {
+            "text/plain": [
+              "S_T*(-D + r)*Piecewise((erf(sqrt(2)*(mu_google - tau)/(2*sigma_google))/2 + 1/2, (Abs(arg(sigma_google)) < pi/4) | ((Abs(arg(sigma_google)) <= pi/4) & (Abs(2*arg(mu_google) - 4*arg(sigma_google) + 2*arg((mu_google - tau)/mu_google) + 2*pi) < pi))), (sqrt(2)*Integral(exp(-(_z - mu_google + tau)**2/(2*sigma_google**2)), (_z, 0, oo))/(2*sqrt(pi)*sigma_google), True))"
+            ],
+            "text/latex": "$\\displaystyle S_{T} \\left(- D + r\\right) \\left(\\begin{cases} \\frac{\\operatorname{erf}{\\left(\\frac{\\sqrt{2} \\left(\\mu_{google} - \\tau\\right)}{2 \\sigma_{google}} \\right)}}{2} + \\frac{1}{2} & \\text{for}\\: \\left(\\left|{\\arg{\\left(\\sigma_{google} \\right)}}\\right| \\leq \\frac{\\pi}{4} \\wedge \\left|{2 \\arg{\\left(\\mu_{google} \\right)} - 4 \\arg{\\left(\\sigma_{google} \\right)} + 2 \\arg{\\left(\\frac{\\mu_{google} - \\tau}{\\mu_{google}} \\right)} + 2 \\pi}\\right| < \\pi\\right) \\vee \\left|{\\arg{\\left(\\sigma_{google} \\right)}}\\right| < \\frac{\\pi}{4} \\\\\\frac{\\sqrt{2} \\int\\limits_{0}^{\\infty} e^{- \\frac{\\left(z - \\mu_{google} + \\tau\\right)^{2}}{2 \\sigma_{google}^{2}}}\\, dz}{2 \\sqrt{\\pi} \\sigma_{google}} & \\text{otherwise} \\end{cases}\\right)$"
+          },
+          "metadata": {},
+          "execution_count": 5
+        }
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "## Local volatility\n",
+        "$\\sigma^2\\left(K, T, S_0\\right)=\\frac{\\frac{\\partial C}{\\partial T}}{\\frac{1}{2} K^2 \\frac{\\partial^2 C}{\\partial K^2}}$\n",
+        "\n",
+        "The right-hand side of this equation can be computed from known European option prices. So, given a complete set of European option prices for all strikes and expirations, local volatilities are given uniquely by equation above.\n",
+        "We can view this equation as a definition of the local volatility function regardless of what kind of process (stochastic volatility for example) actually governs the evolution of volatility."
+      ],
+      "metadata": {
+        "id": "tLMrOdpH8ci9"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "local_volatility = smp.sqrt( dCdT / (1/2*K**2 * varphi))\n",
+        "local_volatility"
+      ],
+      "metadata": {
+        "colab": {
+          "base_uri": "https://localhost:8080/",
+          "height": 141
+        },
+        "id": "6-ut-zDY8aMU",
+        "outputId": "7696d056-f450-4e39-e466-d3f270052379"
+      },
+      "execution_count": 6,
+      "outputs": [
+        {
+          "output_type": "execute_result",
+          "data": {
+            "text/plain": [
+              "sqrt(zoo*S_T*(-D + r)*Piecewise((erf(sqrt(2)*(mu_google - tau)/(2*sigma_google))/2 + 1/2, (Abs(arg(sigma_google)) < pi/4) | ((Abs(arg(sigma_google)) <= pi/4) & (Abs(2*arg(mu_google) - 4*arg(sigma_google) + 2*arg((mu_google - tau)/mu_google) + 2*pi) < pi))), (sqrt(2)*Integral(exp(-(_z - mu_google + tau)**2/(2*sigma_google**2)), (_z, 0, oo))/(2*sqrt(pi)*sigma_google), True)))"
+            ],
+            "text/latex": "$\\displaystyle \\sqrt{\\tilde{\\infty} S_{T} \\left(- D + r\\right) \\left(\\begin{cases} \\frac{\\operatorname{erf}{\\left(\\frac{\\sqrt{2} \\left(\\mu_{google} - \\tau\\right)}{2 \\sigma_{google}} \\right)}}{2} + \\frac{1}{2} & \\text{for}\\: \\left(\\left|{\\arg{\\left(\\sigma_{google} \\right)}}\\right| \\leq \\frac{\\pi}{4} \\wedge \\left|{2 \\arg{\\left(\\mu_{google} \\right)} - 4 \\arg{\\left(\\sigma_{google} \\right)} + 2 \\arg{\\left(\\frac{\\mu_{google} - \\tau}{\\mu_{google}} \\right)} + 2 \\pi}\\right| < \\pi\\right) \\vee \\left|{\\arg{\\left(\\sigma_{google} \\right)}}\\right| < \\frac{\\pi}{4} \\\\\\frac{\\sqrt{2} \\int\\limits_{0}^{\\infty} e^{- \\frac{\\left(z - \\mu_{google} + \\tau\\right)^{2}}{2 \\sigma_{google}^{2}}}\\, dz}{2 \\sqrt{\\pi} \\sigma_{google}} & \\text{otherwise} \\end{cases}\\right)}$"
+          },
+          "metadata": {},
+          "execution_count": 6
+        }
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "source": [
+        "# Using real values from QuantBook to calculate Local Volatility"
+      ],
+      "metadata": {
+        "id": "JSHTTZeU9ONk"
+      }
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "# QuantBook Analysis Tool \n",
+        "# For more information see [https://www.quantconnect.com/docs/v2/our-platform/research/getting-started]\n",
+        "import datetime\n",
+        "qb = QuantBook()\n",
+        "option = qb.AddOption(\"GOOG\") \n",
+        "option.SetFilter(-2, 2, 0, 90)\n",
+        "history = qb.GetOptionHistory(option.Symbol, datetime(2023, 5, 7), datetime(2023, 5, 26))\n",
+        "\n",
+        "all_history = history.GetAllData()\n",
+        "expiries = history.GetExpiryDates() \n",
+        "strikes = history.GetStrikes()\n",
+        "\n",
+        "goog = qb.AddEquity(\"GOOG\")\n",
+        "goog_df = qb.History(qb.Securities.Keys, datetime(2023, 5, 7), datetime(2023, 5, 26), Resolution.Daily)\n",
+        "goog_df = goog_df['close'].unstack(level=0)\n",
+        "S_T = goog_df.values[-1]\n",
+        "all_history.head()"
+      ],
+      "metadata": {
+        "id": "pYG7bmQv88tI"
+      },
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "source": [
+        "#BSM model of Call option\n",
+        "import math\n",
+        "import sympy.stats\n",
+        "\n",
+        "####################### parameters #########################\n",
+        "tau = 3/220 #expiry the 5/19\n",
+        "r = 0.0341 \n",
+        "D = 0.0\n",
+        "K = strikes[-1]\n",
+        "sigma = goog_df.std()\n",
+        "S_T = float(goog_df.values[-1])\n",
+        "mu_goog = goog_df.mean()\n",
+        "#Calculate numerically the local volatility\n",
+        "local_volatility.subs([(S_T,S_T),(K,int(K)),(tau,tau),(D,D),(sigma_google,sigma),(mu_google,mu_google)]).evalf()\n"
+      ],
+      "metadata": {
+        "id": "lcZxAl8K8msU"
+      },
+      "execution_count": null,
+      "outputs": []
+    }
+  ]
+}