quickguide_code.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e297bf72",
   "metadata": {},
   "source": [
    "*This file, part of the* [SdePy](https://github.com/sdepy/sdepy) *package*,\n",
    "*was automatically generated from* `./doc/quickguide.rst`\n",
    "\n",
    "-----------------------------------------------\n",
    "Quick Guide\n",
    "===========\n",
    "\n",
    "\n",
    "------------------\n",
    "Install and import\n",
    "------------------\n",
    "\n",
    "Install using ``pip install sdepy``, or copy the package source code\n",
    "in a directory in your Python path.\n",
    "\n",
    "Import as\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a09a96c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'png'\n",
    "import sdepy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt  # optional, if plots are needed\n",
    "plt.rcParams['figure.figsize'] = (11., 5.5)\n",
    "plt.rcParams['lines.linewidth'] = 1.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d062dded",
   "metadata": {},
   "source": [
    "-------------------\n",
    "How to state an SDE\n",
    "-------------------\n",
    "\n",
    "Here follows a bare-bone definition of a Stochastic Differential\n",
    "Equation (SDE), in this case a Ornstein-Uhlenbeck process:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b2ed027",
   "metadata": {},
   "outputs": [],
   "source": [
    "@sdepy.integrate\n",
    "def my_process(t, x, theta=1., k=1., sigma=1.):\n",
    "    return {'dt': k*(theta - x), 'dw': sigma}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e33b3253",
   "metadata": {},
   "source": [
    "This represents the SDE ``dX = k*(theta - X)*dt + sigma*dW(t)``,\n",
    "where ``theta``, ``k`` and ``sigma`` are parameters and ``dW(t)`` are Wiener\n",
    "process increments. A further ``'dn'`` or ``'dj'`` entry in the returned\n",
    "dictionary would allow for Poisson or compound Poisson jumps.\n",
    "\n",
    "A number of preset processes are provided, including lognormal processes,\n",
    "Hull-White n-factor processes, Heston processes, and jump-diffusion processes.\n",
    "\n",
    "\n",
    "-----------------------\n",
    "How to integrate an SDE\n",
    "-----------------------\n",
    "\n",
    "Now ``my_process`` is a class, a subclass of the cooperating\n",
    "``sdepy.SDE`` and ``sdepy.integrator`` classes:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fdf48bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "issubclass(my_process, sdepy.integrator), issubclass(my_process, sdepy.SDE)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20879777",
   "metadata": {},
   "source": [
    "It is to be instantiated with a number\n",
    "of parameters, including the SDE parameters ``theta``, ``k`` and ``sigma``;\n",
    "its instances are callable, given a timeline they will integrate and\n",
    "return the process along it. Decorating ``my_process`` with ``sdepy.kunfc``\n",
    "allows for more concise handling of parameters:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0be3322",
   "metadata": {},
   "outputs": [],
   "source": [
    "myp = sdepy.kfunc(my_process)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a28a05c",
   "metadata": {},
   "source": [
    "It is best explained by examples, involving ``my_process``, ``myp`` and\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a9ed16",
   "metadata": {},
   "outputs": [],
   "source": [
    "coarse_timeline = (0., 0.25, 0.5, 0.75, 1.0)\n",
    "timeline = np.linspace(0., 1., 500)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b31544af",
   "metadata": {},
   "source": [
    "1. **Scalar process** in 100000 paths, with default parameters, computed\n",
    "   at 5 time points (``coarse_timeline``), using 100 steps in between:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40e0faf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = my_process(x0=1, paths=100*1000,\n",
    "               steps=100)(coarse_timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cbe9e1c",
   "metadata": {},
   "source": [
    "2. The same scalar process computed on a **fine-grained timeline**\n",
    "   (``timeline``) and 1000 paths, using **one integration step for each point\n",
    "   in the timeline** (no ``steps`` parameter):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63e3cf2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = my_process(x0=1, paths=1000,\n",
    "               steps=100)(timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "510cef6e",
   "metadata": {},
   "source": [
    "   A plot of a few paths may be used to **inspect the integration result**:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36734b35",
   "metadata": {},
   "outputs": [],
   "source": [
    "gr = plt.plot(timeline, x[:, :30])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1360f06",
   "metadata": {},
   "source": [
    "2. **Vector process** with three components and\n",
    "   **correlated Wiener increments**\n",
    "   (same other parameters as above):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afde71cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = ((1, .2, -.3), (.2, 1, .1), (-.3, .1, 1))\n",
    "x = my_process(x0=1, vshape=3, corr=corr,\n",
    "               paths=1000)(timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e56fd03",
   "metadata": {},
   "source": [
    "3. Vector process as above, with 10000 paths and\n",
    "   **time-dependent parameters and correlations**:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a289328",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = lambda t: 0.1 + t\n",
    "theta = lambda t: 2-t\n",
    "k = lambda t: 2/(t+1)\n",
    "c02 = lambda t: -0.1*np.cos(3*t)\n",
    "c12 = lambda t: 0.1*np.sign(0.5-t)\n",
    "corr = lambda t: ((     1,    -.2, c02(t)),\n",
    "                  (   -.2,      1, c12(t)),\n",
    "                  (c02(t), c12(t),      1))\n",
    "x = my_process(x0=1, vshape=3, corr=corr,\n",
    "               theta=theta, k=k, sigma=sigma, paths=10*1000)(timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d19c81d",
   "metadata": {},
   "source": [
    "   This plot illustrates the correlations among the components of\n",
    "   ``x`` increments, as a function of time and as compared to\n",
    "   ``corr(t)``:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b9a5dbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = np.diff(x, axis=0)\n",
    "for i in range(3):\n",
    "    for j in range(i + 1, 3):\n",
    "        gr = plt.plot(\n",
    "            timeline, corr(timeline)[i][j] + 0*timeline,\n",
    "            timeline[1:], [np.cov(z)[i, j]/(z[i].std()*z[j].std())\n",
    "                           for z in dx]\n",
    "            )\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b05368db",
   "metadata": {},
   "source": [
    "4. A 1000 paths scalar process with **path-dependent initial conditions\n",
    "   and parameters**, integrated **backwards** (``i0=-1``):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15bdd008",
   "metadata": {},
   "outputs": [],
   "source": [
    "x0, sigma = np.zeros(1000), np.zeros(1000)\n",
    "x0[::2], x0[1::2] = 0., 2.\n",
    "sigma[::2], sigma[1::2] = 0.5, 0.1\n",
    "x = my_process(x0=x0, sigma=sigma, paths=1000,\n",
    "               theta=1, k=-2,\n",
    "               i0=-1)(timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fa48acb",
   "metadata": {},
   "source": [
    "   When integrating backwards, the inital conditions are applied\n",
    "   at the final point in the given timeline:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41ff53f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert (x[-1, :] == x0).all()\n",
    "gr = plt.plot(timeline, x[:, :30])\n",
    "gr = plt.plot(timeline, np.full_like(timeline, 1), 'k--')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20a4e7bd",
   "metadata": {},
   "source": [
    "   Note the negative value of ``k``, with mean reversion\n",
    "   towards ``theta=1`` occurring backwards in time.\n",
    "\n",
    "5. A scalar process computed on a **10 x 15 grid of parameters** ``sigma`` and\n",
    "   ``k`` (note that the shape of the initial conditions and of each\n",
    "   parameter should be broadcastable to the values of the process across\n",
    "   paths, i.e. to shape ``vshape + (paths,)``):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a4bd3b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = np.linspace(0., 1., 10).reshape(10, 1, 1)\n",
    "k = np.linspace(1., 2., 15).reshape(1, 15, 1)\n",
    "x = my_process(x0=1, theta=2, k=k, sigma=sigma, vshape=(10, 15),\n",
    "               paths=10*1000)(coarse_timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7aa4a5b",
   "metadata": {},
   "source": [
    "   A plot of the final average process values against ``k``\n",
    "   illustrates a faster reversion to ``theta=2`` as\n",
    "   ``k`` increases, as well as the independence of the\n",
    "   process mean from ``sigma``:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b48dd969",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(10):\n",
    "    gr = plt.plot(k[0, :, 0], x[-1, i, :, :].mean(axis=-1))\n",
    "lb = plt.xlabel('k'), plt.ylabel('x(t=2).mean()')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5eb8805",
   "metadata": {},
   "source": [
    "   In the example above, set ``steps>=100`` to go from inaccurate and fast,\n",
    "   to meaningful and slow.\n",
    "\n",
    "6. **Interactive modification** of process and integration parameters\n",
    "   using the ``sdepy.kfunc`` decorator ``myp = sdepy.kfunc(my_process)``.\n",
    "\n",
    "   The ``sdepy.kfunc`` decorated version of ``my_process`` is a subclass\n",
    "   of ``sdepy.integrator`` and ``sdepy.SDE``, as ``my_process`` is,\n",
    "   and fully replicates its functionality and interface:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9573e722",
   "metadata": {},
   "outputs": [],
   "source": [
    "issubclass(myp, sdepy.integrator), issubclass(myp, sdepy.SDE)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a89deeb",
   "metadata": {},
   "source": [
    "   In addition, and in contrast to ``my_process``,\n",
    "   ``myp`` instances accept either an integration timeline,\n",
    "   or a modified value of some integration or SDE parameters,\n",
    "   or both, as illustrated below:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a555a02",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = myp(x0=1, sigma=1, paths=1000)\n",
    "\n",
    "x = p(timeline)\n",
    "x1, x2 = p(timeline, sigma=0.5), p(timeline, sigma=1.5)\n",
    "q = p(paths=100, vshape=(3,), k=2)\n",
    "y = q(timeline, sigma=0.5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e782b1e7",
   "metadata": {},
   "source": [
    "   ``x`` is the result of integrating ``p`` along ``timeline``\n",
    "   (no difference here from a ``my_process`` instance);\n",
    "   ``x1, x2`` are obtained by integration along ``timeline``\n",
    "   by setting ``sigma`` to the given values,\n",
    "   and keeping other parameters as stated when ``p`` was instantiated;\n",
    "   ``q`` is another ``myp`` instance with updated default values\n",
    "   for ``paths``, ``vshape`` and ``k``, and all else set as in ``p``;\n",
    "   and finally, ``y`` was obtained by integrating ``q`` along ``timeline``,\n",
    "   with its own parameters, save for ``sigma`` that was modified to ``0.5``.\n",
    "\n",
    "   Moreover, for ``sdepy.kfunc`` classes, instantiation and computation\n",
    "   may happen contextually:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdc10ba7",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = myp(timeline, x0=1, sigma=1, paths=1000)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4558e090",
   "metadata": {},
   "source": [
    "   is equivalent to:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e819c6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = my_process(x0=1, sigma=1, paths=1000)(timeline)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46e631f4",
   "metadata": {},
   "source": [
    "   ``sdepy.kfunc``-decorated classes allow to **state some central values\n",
    "   of parameters for a given problem**, and to explore the effects\n",
    "   of variatons in some of them via a concise interface,\n",
    "   that **keeps the modified parameters in focus**\n",
    "   and all the rest in the background.\n",
    "\n",
    "   To **inspect the parameters stored** in a ``sdepy.kfunc`` instance,\n",
    "   use the read-only ``params`` attribute:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a97cd7c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "q.params\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e78686e5",
   "metadata": {},
   "source": [
    "   To test if an object is a kfunc, use ``sdepy.iskfunc()``:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2317d65",
   "metadata": {},
   "outputs": [],
   "source": [
    "sdepy.iskfunc(myp), sdepy.iskfunc(p), sdepy.iskfunc(my_process)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e8f605f",
   "metadata": {},
   "source": [
    "   The examples that follow illustrate, among other things, the\n",
    "   use of ``myp`` as a ``sdepy.kfunc`` class.\n",
    "\n",
    "7. Processes generated using **integration results as stochasticity sources**\n",
    "   (mind using consistent ``vshape`` and ``paths``, and synchronizing timelines):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8736e5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_dw = sdepy.integrate(lambda t, x: {'dw': 1})(vshape=1, paths=1000)(timeline)\n",
    "p = myp(dw=my_dw, vshape=3, paths=1000,\n",
    "        x0=1, sigma=((1,), (2,), (3,)))\n",
    "x = p(timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcb3f138",
   "metadata": {},
   "source": [
    "   Now, ``x1, x2, x3 = = x[:, 0], x[:, 1], x[:, 2]`` have different ``sigma``,\n",
    "   but share the same ``dw`` increments, as can be seen plotting a path:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f56ae385",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 0  # path to be plotted\n",
    "gr = plt.plot(timeline, x[:, :, k])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68bac7fa",
   "metadata": {},
   "source": [
    "   If more integrations steps are needed between points in the output timeline,\n",
    "   use ``steps`` to keep the integration timeline consistent with the one\n",
    "   of ``my_dw``:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0857895",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = p(coarse_timeline, steps=timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dddd2a31",
   "metadata": {},
   "source": [
    "8. Using **stochasticity sources with memory**\n",
    "   (mind using consistent ``vshape`` and ``paths``):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ed665b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_dw = sdepy.true_wiener_source(paths=1000)\n",
    "p = myp(x0=1, theta=1, k=1, sigma=1, dw=my_dw, paths=1000)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6909a572",
   "metadata": {},
   "source": [
    "   ``my_dw``, as a ``sdepy.true_wiener_source`` instance\n",
    "   has memory of, and generates new Wiener process increments\n",
    "   consistent with, its formerly realized values.\n",
    "   As a consequence, processes defined invoking ``p``\n",
    "   share the same underlying Wiener process increments:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1efb8714",
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = np.linspace(0., 1.,  30)\n",
    "t2 = np.linspace(0., 1., 100)\n",
    "t3 = t = np.linspace(0., 1., 300)\n",
    "x1, x2, x3 = p(t1), p(t2), p(t3)\n",
    "y1, y2, y3 = p(t, theta=1.5), p(t, theta=1.75), p(t, theta=2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e98052f",
   "metadata": {},
   "source": [
    "   ``x1, x2, x3`` illustrate SDE integration convergence as time steps\n",
    "   become smaller, and ``y1, y2, y3`` illustrate\n",
    "   how ``theta`` affects paths, all else being equal:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d104367d",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0 # path to be plotted\n",
    "gr = plt.plot(t, x1(t)[:, i], t, x2(t)[:, i], t, x3(t)[:, i])\n",
    "plt.show() \n",
    "gr = plt.plot(t, y1[:, i], t, y2[:, i], t, y3[:, i])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1280c49c",
   "metadata": {},
   "source": [
    "------------------------------------\n",
    "How to handle the integration output\n",
    "------------------------------------\n",
    "\n",
    "SDE integrators return instances of ``sdepy.process``,\n",
    "a subclass of ``np.ndarray``\n",
    "with a timeline stored in the ``t`` attribute\n",
    "(note the shape of ``x``,\n",
    "repeatedly used in the examples below):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24d9161d",
   "metadata": {},
   "outputs": [],
   "source": [
    "coarse_timeline = (0., 0.25, 0.5, 0.75, 1.0)\n",
    "timeline = np.linspace(0., 1., 101)\n",
    "x = my_process(x0=1, vshape=3, paths=1000)(timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0f21e41",
   "metadata": {},
   "source": [
    "``x`` is a ``sdepy.process`` instance:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "416f5074",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c7a553b",
   "metadata": {},
   "source": [
    "and is based on the given timeline:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1a2a4b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isclose(timeline, x.t).all()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d3e437b",
   "metadata": {},
   "source": [
    "Whenever possible, a process will store references, not copies, of timeline\n",
    "and values. In fact:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17a5e733",
   "metadata": {},
   "outputs": [],
   "source": [
    "timeline is x.t\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aed76f2f",
   "metadata": {},
   "source": [
    "The first axis is reserved for the timeline, the last for paths, and axes\n",
    "in the middle match the shape of process values:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2939976",
   "metadata": {},
   "outputs": [],
   "source": [
    "x.shape == x.t.shape + x.vshape + (x.paths,)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e668b2eb",
   "metadata": {},
   "source": [
    "Calling processes interpolates in time:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b5fbfaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x(coarse_timeline)\n",
    "y.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7e85842",
   "metadata": {},
   "source": [
    "The result is always an array, not a process:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0510e8b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(y)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8d03993",
   "metadata": {},
   "source": [
    "Indexing works as usual, and returns NumPy arrays:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0de3ef9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(x[0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab87d1a8",
   "metadata": {},
   "source": [
    "All array methods are unchanged  (no overriding),\n",
    "and return NumPy arrays as well:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeef3a18",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(x.mean(axis=0))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d357337",
   "metadata": {},
   "source": [
    "You can slice processes along time, values and paths with special indexing.\n",
    "\n",
    "- Time indexing:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ce67ae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x['t', ::2]\n",
    "y.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5285162",
   "metadata": {},
   "source": [
    "- Values indexing:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29928688",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x['v', 0]\n",
    "y.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49a9b333",
   "metadata": {},
   "source": [
    "- Paths indexing:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66beee32",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x['p', :10]\n",
    "y.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e57dac65",
   "metadata": {},
   "source": [
    "The output of a special indexing operation is a process:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c50e3709",
   "metadata": {},
   "outputs": [],
   "source": [
    "isinstance(y, sdepy.process)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "870f3648",
   "metadata": {},
   "source": [
    "Smart indexing is allowed. To select paths that cross ``x=0``\n",
    "at some point and for some component, use:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbbe921d",
   "metadata": {},
   "outputs": [],
   "source": [
    "i_negative = x.min(axis=(0, 1)) < 0\n",
    "y = x['p', i_negative]\n",
    "y.shape == (101, 3, i_negative.sum())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05f5befb",
   "metadata": {},
   "source": [
    "You can do algebra with processes that either share the same timeline, or are constant\n",
    "(a process with a one-point timeline is assumed to be constant), and either have the\n",
    "same number of paths, or are deterministic (with one path):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad3593ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_const = x['t', 0]  # a constant process\n",
    "x_one_path = x['p', 0]  # a process with one path\n",
    "\n",
    "y = np.exp(x) - x_const\n",
    "z = np.maximum(x, x_one_path)\n",
    "\n",
    "isinstance(y, sdepy.process), isinstance(z, sdepy.process)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08eafcf2",
   "metadata": {},
   "source": [
    "When integrating SDEs, the SDE parameters and/or stochasticity sources\n",
    "accept processes as valid values (mind using deterministic processes, or\n",
    "synchronizing the number of paths, and make sure that the shape of values\n",
    "do broadcast together). To use a realization of ``my_process``\n",
    "as the volatility of a 3-component lognormal process, do as follows:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2d3d8b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "stochastic_vol = my_process(x0=1, paths=10*1000)(timeline)\n",
    "stochastic_vol_x = sdepy.lognorm_process(x0=1, vshape=3, paths=10*1000,\n",
    "    mu=0, sigma=stochastic_vol)(timeline)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6a4931a",
   "metadata": {},
   "source": [
    "Processes have specialized methods, and may be analyzed, and their statistics\n",
    "cumulated across multiple runs, using the ``sdepy.montecarlo`` class. Some examples follow:\n",
    "\n",
    "1. Cumulative probability distribution function at t=0.5\n",
    "   of the process values of ``x`` across paths:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bec817da",
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf = x.cdf(0.5, x=np.linspace(-2, 2, 100))  # an array\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b737957",
   "metadata": {},
   "source": [
    "2. Characteristic function at t=0.5 of the same distribution:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48988241",
   "metadata": {},
   "outputs": [],
   "source": [
    "chf = x.chf(0.5, u=np.linspace(-2, 2, 100))  # an array\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3f0113b",
   "metadata": {},
   "source": [
    "3. Standard deviation across paths:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31e930b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "std = x.pstd()  # a one-path process\n",
    "std.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e1937cc",
   "metadata": {},
   "source": [
    "4. Maximum value reached along the timeline:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79a6ddf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "xmax = x.tmax()  # a constant process\n",
    "xmax.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3590c20",
   "metadata": {},
   "source": [
    "5. A linearly interpolated, or Gaussian kernel estimate (default)\n",
    "   of the probability distribution function (pdf) and its cumulated\n",
    "   values (cdf) across paths, at a given time point,\n",
    "   may be obtained using the ``montecarlo`` class:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a32a3bf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x(1)[0]  # 0-th component of x at time t=1\n",
    "a = sdepy.montecarlo(y, bins=30)\n",
    "ygrid = np.linspace(y.min(), y.max(), 200)\n",
    "gr = plt.plot(ygrid, a.pdf(ygrid), ygrid, a.cdf(ygrid))\n",
    "gr = plt.plot(ygrid, a.pdf(ygrid, method='interp', kind='nearest'))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c58cd6a",
   "metadata": {},
   "source": [
    "6. A ``sdepy.montecarlo`` instance can be used to cumulate the results\n",
    "   of multiple simulations, across multiple components of process values:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa42fdd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = my_process(x0=1, vshape=3, paths=10*1000)\n",
    "a = sdepy.montecarlo(bins=100)  # empty montecarlo instance\n",
    "for _ in range(10):\n",
    "    x = p(timeline)  # run simulation\n",
    "    a.update(x(1))  # cumulate x values at t=1\n",
    "a.paths\n",
    "gr = plt.plot(ygrid, a[0].pdf(ygrid), ygrid, a[0].cdf(ygrid))\n",
    "gr = plt.plot(ygrid, a[0].pdf(ygrid, method='interp', kind='nearest'))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79855fdb",
   "metadata": {},
   "source": [
    "--------------------------------\n",
    "Example - Stochastic Runge-Kutta\n",
    "--------------------------------\n",
    "\n",
    "Minimal implementation of a basic stochastic Runge-Kutta integration\n",
    "scheme, as a subclass of ``sdepy.integrator``\n",
    "(the ``A`` and ``dZ`` methods below are the standardized way\n",
    "in which equations are exposed to integrators):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "248bb75a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import sqrt\n",
    "class my_integrator(sdepy.integrator):\n",
    "    def next(self):\n",
    "        t, new_t = self.itervars['sw']\n",
    "        x, new_x = self.itervars['xw']\n",
    "        dt = new_t - t\n",
    "        A, dZ = self.A(t, x), self.dZ(t, dt)\n",
    "        a, b, dw = A['dt'], A['dw'], dZ['dw']\n",
    "        b1 = self.A(t, x + a*dt + b*sqrt(dt))['dw']\n",
    "        new_x[...] = x + a*dt + b*dw + (b1 - b)/2 * (dw**2 - dt)/sqrt(dt)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62c7ae7d",
   "metadata": {},
   "source": [
    "SDE of a lognormal process, as a subclass of ``sdepy.SDE``,\n",
    "and classes that integrate it with the default integration method\n",
    "(``euler``) and via ``my_integrator`` (``rk``):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36cb3727",
   "metadata": {},
   "outputs": [],
   "source": [
    "class my_SDE(sdepy.SDE):\n",
    "    def sde(self, t, x):\n",
    "        return {'dt': 0, 'dw': x}\n",
    "\n",
    "class euler(my_SDE, sdepy.integrator):\n",
    "    pass\n",
    "\n",
    "class rk(my_SDE, my_integrator):\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6731b27f",
   "metadata": {},
   "source": [
    "Comparison of integration errors, as the integration from ``t=0`` to\n",
    "``t=1`` is carried out with an increasing number of steps, against\n",
    "the integration result of ``sdepy.lognorm_process``, which returns\n",
    "an exact result irrespective of the number and size\n",
    "of the integration steps (this happens since, by implementation,\n",
    "it integrates the linear SDE for ``log(x)``):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc64244f",
   "metadata": {},
   "outputs": [],
   "source": [
    "args = dict(dw=sdepy.true_wiener_source(paths=100),\n",
    "            paths=100, x0=10)\n",
    "timeline = (0, 1)\n",
    "steps = np.array((2, 3, 5, 10, 20, 30, 50, 100,\n",
    "                  200, 300, 500, 1000, 2000, 3000))\n",
    "\n",
    "# exact integration results at t=1\n",
    "exact = sdepy.lognorm_process(mu=0, sigma=1, **args)(timeline)[-1].mean()\n",
    "\n",
    "# errors of approximate integration results at t=1\n",
    "errors = np.abs(np.array([\n",
    "    [euler(**args, steps=s)(timeline)[-1].mean()/exact - 1,\n",
    "     rk(**args, steps=s)(timeline)[-1].mean()/exact - 1]\n",
    "    for s in steps]))\n",
    "\n",
    "# plots\n",
    "ax = plt.axes(label=0); ax.set_xscale('log'); ax.set_yscale('log')\n",
    "gr = ax.plot(steps, errors)\n",
    "plt.show()  \n",
    "print(f'euler error: {errors[-1,0]:.2e}\\n'\n",
    "      f'   rk error: {errors[-1,1]:.2e}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ab8a6d5",
   "metadata": {},
   "source": [
    "--------------------------------\n",
    "Example - Fokker-Planck Equation\n",
    "--------------------------------\n",
    "\n",
    "Monte Carlo integration of partial differential equations, illustrated\n",
    "in the simplest example of the heat equation ``diff(u, t) - k*diff(u, x, 2) == 0``,\n",
    "for the function ``u(x, t)``, i.e. the Fokker-Planck equation for the SDE\n",
    "``dX(t) = sqrt(2*k)*dW(t)``. Initial conditions at ``t=t0``, two examples:\n",
    "\n",
    "1.  ``u(x, t0) = 1`` for ``lb < x < hb`` and ``0`` otherwise,\n",
    "2.  ``u(x, t0) = sin(x)``.\n",
    "\n",
    "Setup:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7958dbfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import exp, sin\n",
    "from scipy.special import erf\n",
    "from scipy.integrate import quad\n",
    "\n",
    "k = .5\n",
    "x0, x1 = 0, 10;\n",
    "t0, t1 = 0, 1\n",
    "lb, hb = 4, 6\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d65019e4",
   "metadata": {},
   "source": [
    "Exact green function and solutions for initial conditions\n",
    "1. and 2., to be checked against results:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87b037e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def green_exact(y, s, x, t):\n",
    "    return exp(-(x - y)**2/(4*k*(t - s)))/sqrt(4*np.pi*k*(t - s))\n",
    "\n",
    "def u1_exact(x, t):\n",
    "    return (erf((x - lb)/2/sqrt(k*(t - t0))) - erf((x - hb)/2/sqrt(k*(t - t0))))/2\n",
    "\n",
    "def u2_exact(x, t):\n",
    "    return exp(-k*(t - t0))*sin(x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c1d8b55",
   "metadata": {},
   "source": [
    "Realization of the needed stochastic process, by backward integration from\n",
    "a grid of final values of ``x`` at ``t=t1``, using the preset\n",
    "``wiener_process`` class (the ``steps`` keyword is added as a reminder\n",
    "of the setup needed for less-than-trivial equations, it does not actually\n",
    "make a difference here):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5eee7967",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgrid = np.linspace(x0, x1, 51)\n",
    "tgrid = np.linspace(t0, t1, 5)\n",
    "xp = sdepy.wiener_process(\n",
    "    paths=10000, steps=100,\n",
    "    sigma=sqrt(2*k),\n",
    "    vshape=xgrid.shape, x0=xgrid[..., np.newaxis],\n",
    "    i0=-1,\n",
    "    )(timeline=tgrid)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c89bf13",
   "metadata": {},
   "source": [
    "Computation of the green function and of the solutions ``u(x, t1)``\n",
    "via Monte Carlo integration\n",
    "(note the liberal use of ``scipy.integrate.quad`` below, enabled by\n",
    "the smoothness of the Gaussian kernel estimate ``a[i, j].pdf``):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fb3901f",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = sdepy.montecarlo(xp, bins=100)\n",
    "\n",
    "def green(y, i, j):\n",
    "    \"\"\"green function from (y=y, s=tgrid[i]) to (x=xgrid[j], t=t1)\"\"\"\n",
    "    return a[i, j].pdf(y)\n",
    "\n",
    "u1, u2 = np.empty(51), np.empty(51)\n",
    "for j in range(51):\n",
    "    u1[j] = quad(lambda y: green(y, 0, j), lb, hb)[0]\n",
    "    u2[j] = quad(lambda y: sin(y)*green(y, 0, j), -np.inf, np.inf)[0]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe6bfa78",
   "metadata": {},
   "source": [
    "Comparison against exact values:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52dee8b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.linspace(x0, x1, 500)\n",
    "for i, j in ((1, 20), (2, 30), (3, 40)):\n",
    "    gr = plt.plot(y, green(y, i, j),\n",
    "                  y, green_exact(y, tgrid[i], xgrid[j], t1), ':')\n",
    "plt.show()  \n",
    "\n",
    "gr = plt.plot(xgrid, u1, y, u1_exact(y, t1), ':')\n",
    "gr = plt.plot(xgrid, u2, y, u2_exact(y, t1), ':')\n",
    "plt.show()  \n",
    "\n",
    "print(f'u1 error: {np.abs(u1 - u1_exact(xgrid, t1)).mean():.2e}\\n'\n",
    "      f'u2 error: {np.abs(u2 - u2_exact(xgrid, t1)).mean():.2e}'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18dc864c",
   "metadata": {},
   "source": [
    "--------------------------------\n",
    "Example - Basket Lookback Option\n",
    "--------------------------------\n",
    "\n",
    "Take a basket of 4 financial securities, with risk-neutral probabilities following\n",
    "lognormal processes in the Black-Scholes framework. Correlations, dividend yields\n",
    "and term structure of volatility (will be linearly interpolated) are given below:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c87e370",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = [\n",
    "    [1,    0.50, 0.37, 0.35],\n",
    "    [0.50,    1, 0.47, 0.46],\n",
    "    [0.37, 0.47,    1, 0.19],\n",
    "    [0.35, 0.46,  0.19,   1]]\n",
    "\n",
    "dividend_yield = sdepy.process(c=(0.20, 4.40, 0., 4.80))/100\n",
    "riskfree = 0  # to keep it simple\n",
    "\n",
    "vol_timepoints = (0.1, 0.2, 0.5, 1, 2, 3)\n",
    "vol = np.array([\n",
    "    [0.40, 0.38, 0.30, 0.28, 0.27, 0.27],\n",
    "    [0.31, 0.29, 0.22, 0.16, 0.18, 0.21],\n",
    "    [0.24, 0.22, 0.19, 0.19, 0.21, 0.22],\n",
    "    [0.35, 0.31, 0.21, 0.18, 0.19, 0.19]])\n",
    "sigma = sdepy.process(t=vol_timepoints, v=vol.T)\n",
    "sigma.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63e7f599",
   "metadata": {},
   "source": [
    "The prices of the securities at the end of each quarter for the next 2 years,\n",
    "simulated across 50000 independent paths and their antithetics\n",
    "(``sdepy.odd_wiener_source`` is used), are:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78a1b8c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "maturity = 2\n",
    "timeline = np.linspace(0, maturity, 4*maturity + 1)\n",
    "p = sdepy.lognorm_process(\n",
    "    x0=100, corr=corr, dw=sdepy.odd_wiener_source,\n",
    "    mu=(riskfree - dividend_yield),\n",
    "    sigma=sigma,\n",
    "    vshape=4, paths=100*1000, steps=maturity*250)\n",
    "x = p(timeline)\n",
    "x.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2742ef32",
   "metadata": {},
   "source": [
    "A call option knocks in if any of the securities reaches a price below 80\n",
    "at any quarter (starting from 100), and pays the lookback maximum attained\n",
    "by the basket (equally weighted), minus 105, if positive.\n",
    "Its price is:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb57edc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_worst = x.min(axis=1)\n",
    "x_basket = x.mean(axis=1)\n",
    "down_and_in_paths = (x_worst.min(axis=0) < 80)\n",
    "lookback_x_basket = x_basket.max(axis=0)\n",
    "payoff = np.maximum(0, lookback_x_basket - 105)\n",
    "payoff[np.logical_not(down_and_in_paths)] = 0\n",
    "a = sdepy.montecarlo(payoff, use='even')\n",
    "print(a)\n"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}