Market Making in Short-Term Alpha (Chapter 10.4.2 of Algorithmic and High-Frequency Trading by Cartea, Jaimungal, Penalva, published by Cambridge University Press)

short_term_alpha.ipynb

shot_term_alpha_helpers.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import necessary libraries\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib as mpl\n",
    "from matplotlib import cm\n",
    "from matplotlib.ticker import LinearLocator\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.pylab as pylab\n",
    "params = {'legend.fontsize': 10,\n",
    "          'figure.figsize': (8, 4),\n",
    "         'axes.labelsize': 20,\n",
    "         'axes.titlesize': 20,\n",
    "         'xtick.labelsize': 15,\n",
    "         'ytick.labelsize': 15}\n",
    "pylab.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# these are helper functions for the main solver\n",
    "\n",
    "# for the various constants of interest from the alpha process\n",
    "class ShortTermAlpha:\n",
    "\n",
    "    def __init__(self, zeta, epsilon, eta):\n",
    "        self.zeta = zeta\n",
    "        self.epsilon = epsilon\n",
    "        self.eta = eta\n",
    "        self.dalpha = None\n",
    "        self.value = None\n",
    "\n",
    "    def generate_dalpha(self, dt, lambda_p, lambda_m):\n",
    "        self.dalpha = math.sqrt((self.eta ** 2 + self.epsilon ** 2 * (lambda_p + lambda_m)) * 3 * dt)\n",
    "\n",
    "    def generate_alpha(self, alpha_threshold, decimal_places=4):\n",
    "        if self.dalpha is not None:\n",
    "            max_alpha = math.ceil(alpha_threshold / self.dalpha) * self.dalpha\n",
    "            min_alpha = -max_alpha\n",
    "            self.value = np.round(np.arange(min_alpha, max_alpha + self.dalpha, self.dalpha), decimal_places)\n",
    "            \n",
    "# for displaying progress in solving the DPE\n",
    "def counter(current_idx, total_idx):\n",
    "    M = round(total_idx/10)\n",
    "    percent = int(round(current_idx/total_idx*100,0))\n",
    "    if np.mod(current_idx,M) == 0 :\n",
    "        print(\"Processing: \" + str(percent) + \" percent completed.\")\n",
    "\n",
    "\n",
    "def nan_to_num(x, nan):\n",
    "    \"\"\" Change the NaNs in a numpy array to the desired values.\n",
    "    :param x: a numpy array\n",
    "    :param nan: desired value\n",
    "    :return: a deep copy of x with changed values\n",
    "    \"\"\"\n",
    "    y = np.copy(x)\n",
    "    y[np.isnan(y)] = nan\n",
    "    return y\n",
    "\n",
    "\n",
    "def linear_extrapolation(x, y, e):\n",
    "    \"\"\" Extrapolation and interpolation\n",
    "    \n",
    "    :param x: a numpy array\n",
    "    :param y: a numpy array\n",
    "    :param e: a numpy array, equivalent of x\n",
    "    :return: a numpy array\n",
    "    \"\"\"\n",
    "    new_x = np.sort(x)\n",
    "    new_y = y[np.argsort(x)]\n",
    "\n",
    "    def point_wise(ep):\n",
    "        if ep < new_x[0]:\n",
    "            return new_y[0] + (ep - new_x[0]) * (new_y[1] - new_y[0]) / (new_x[1] - new_x[0])\n",
    "        elif ep > new_x[-1]:\n",
    "            return new_y[-1] + (ep - new_x[-1]) * (new_y[-1] - new_y[-2]) / (new_x[-1] - new_x[-2])\n",
    "        else:\n",
    "            return np.interp([ep], x, y)[0]\n",
    "\n",
    "    return np.array([point_wise(i) for i in e])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_h_matrix(short_term_alpha, q, t, Delta, varphi, dt, phi, lambda_p, lambda_m):\n",
    "    \"\"\" Solves the DPE on page 269.\n",
    "    \n",
    "    :param short_term_alpha: a short_term_alpha object\n",
    "    :param q: an inventory vector\n",
    "    :param t: a time \n",
    "    :param Delta: variable\n",
    "    :param varphi: variable\n",
    "    :param dt: variable\n",
    "    :param phi: variable\n",
    "    :param lambda_p: a numpy vector\n",
    "    :param lambda_m: a numpy vector\n",
    "    :return: h value matrix, l+ and l- of Equation 10.32\n",
    "    \"\"\"\n",
    "    \n",
    "    # set up the time, inventory, and space grids\n",
    "    dalpha = short_term_alpha.dalpha\n",
    "    alpha = short_term_alpha.value\n",
    "    eta = short_term_alpha.eta\n",
    "    zeta = short_term_alpha.zeta\n",
    "    Nalpha = alpha.shape[0]\n",
    "    epsilon = short_term_alpha.epsilon\n",
    "    alpha_up = alpha + epsilon\n",
    "    alpha_down = alpha - epsilon\n",
    "\n",
    "    Nt = t.shape[0]\n",
    "    Nq = q.shape[0]\n",
    "\n",
    "    # stores the h as a function of q, alpha, t\n",
    "    h = np.full((Nq, Nalpha, Nt), fill_value=np.nan)\n",
    "    \n",
    "    # stores the optimal posting strategies as a function of q, alpha, t\n",
    "    lp = np.zeros((Nq, Nalpha, Nt), dtype=bool)\n",
    "    lm = np.zeros((Nq, Nalpha, Nt), dtype=bool)\n",
    "\n",
    "    # Terminal conditions for all q\n",
    "    h[:, :, h.shape[-1] - 1] = np.around((-(0.5 * Delta * np.sign(q) + varphi * q) * q), decimals=4)[:, np.newaxis]\n",
    "\n",
    "    d_alpha_h = np.zeros(Nalpha)\n",
    "\n",
    "    # Index of alpha smaller than 0\n",
    "    idx_m = np.where(alpha < 0)[0]\n",
    "    # Index of alpha greater than 0\n",
    "    idx_p = np.where(alpha > 0)[0]\n",
    "    # Index of alpha equals to 0\n",
    "    idx_0 = np.where(alpha == 0)[0]\n",
    "\n",
    "    for i in range(h.shape[2] - 2, -1, -1):\n",
    "\n",
    "        # print counter\n",
    "        counter(Nt - i, Nt)\n",
    "\n",
    "        # Solve the HJB in the continuation region  backward from t + dt to t\n",
    "        for k in range(0, h.shape[0], 1):\n",
    "            \n",
    "            # compute the optimal strategy over the time interval [t,t+dt) Eq. (10.32)\n",
    "            # interpolate/extrapolate to alpha + eta\n",
    "            h_p2 = linear_extrapolation(alpha, h[k, :, i + 1], alpha_up)\n",
    "            if q[k] > -(Nq - 1) / 2:\n",
    "                h_p1 = linear_extrapolation(alpha, h[k - 1, :, i + 1], alpha_up)\n",
    "                post_h_p = 0.5 * Delta + h_p1 - h_p2\n",
    "            else:\n",
    "                post_h_p = np.zeros(Nalpha)\n",
    "\n",
    "            # interpolate/extrapolate to alpha - eta\n",
    "            h_m2 = linear_extrapolation(alpha, h[k, :, i + 1], alpha_down)\n",
    "            if q[k] < (Nq - 1) / 2:\n",
    "                h_m1 = linear_extrapolation(alpha, h[k + 1, :, i + 1], alpha_down)\n",
    "                post_h_m = 0.5 * Delta + h_m1 - h_m2\n",
    "            else:\n",
    "                post_h_m = np.zeros(Nalpha)\n",
    "\n",
    "            lp[k, :, i + 1] = post_h_p > 0\n",
    "            lm[k, :, i + 1] = post_h_m > 0\n",
    "            \n",
    "            \n",
    "            # solve DPE for h function one time-step backwards\n",
    "            d_alpha_h[idx_m] = (h[k, idx_m + 1, i + 1] - h[k, idx_m, i + 1]) / dalpha\n",
    "            d_alpha_h[idx_p] = (h[k, idx_p, i + 1] - h[k, idx_p - 1, i + 1]) / dalpha\n",
    "            d_alpha_h[idx_0] = (h[k, idx_0 + 1, i + 1] - h[k, idx_0 - 1, i + 1]) / (2 * dalpha)\n",
    "\n",
    "            d2_alpha_h = (h[k, 2:h.shape[1], i + 1] - 2 * h[k, 1:(h.shape[1] - 1), i + 1] + h[k, 0:(h.shape[1] - 2),\n",
    "                                                                                            i + 1]) / (dalpha ** 2)\n",
    "\n",
    "            h[k, 1:(h.shape[1] - 1), i] = h[k, 1:(h.shape[1] - 1), i + 1] \\\n",
    "                                          + dt * (-zeta * alpha[1:(alpha.shape[0] - 1)] * d_alpha_h[\n",
    "                                                                                          1:(d_alpha_h.shape[0] - 1)]\n",
    "                                                  + 0.5 * (eta ** 2) * d2_alpha_h\n",
    "                                                  + zeta * alpha[1:(alpha.shape[0] - 1)] * q[k] - phi * (q[k] ** 2)\n",
    "                                                  + lambda_p * np.maximum(post_h_p[1:(post_h_p.shape[0] - 1)], 0)\n",
    "                                                  + lambda_m * np.maximum(post_h_m[1:(post_h_m.shape[0] - 1)], 0)\n",
    "                                                  + lambda_p * (h_p2[1:(h_p2.shape[0] - 1)] - h[k, 1:(h.shape[1] - 1), i + 1])\n",
    "                                                  + lambda_m * (h_m2[1:(h_m2.shape[0] - 1)] - h[k, 1:(h.shape[1] - 1), i + 1]))\n",
    "\n",
    "            # impose second derivative vanishes along maximum and minimum values of alpha grid\n",
    "            h[k, h.shape[1] - 1, i] = 2 * h[k, h.shape[1] - 2, i] - h[k, h.shape[1] - 3, i]\n",
    "            h[k, 0, i] = 2 * h[k, 1, i] - h[k, 2, i]\n",
    "            \n",
    "    print(\"done!\")\n",
    "\n",
    "    return h, lp, lm\n",
    "\n",
    "def generate_sell_buy_posts(short_term_alpha, t, q, lp, lm):\n",
    "    \"\"\" Calculate the inventory level above/below which it is optimal to post buy/sell limit orders\n",
    "    based on lp amd lm solved from the DPE previously.\n",
    "    \n",
    "    :param short_term_alpha: a short_term_alpha object\n",
    "    :param t: a time vector\n",
    "    :param q: an inventory vector\n",
    "    :param lp: a 3d matrix\n",
    "    :param lm: a 3d mrtix\n",
    "    :return: sell side posts and buy side posts\n",
    "    \"\"\"\n",
    "    Nt = t.shape[0]\n",
    "    Nalpha = short_term_alpha.value.shape[0]\n",
    "    lp_b = np.full((Nalpha, Nt), np.NaN)\n",
    "    lm_b = np.full((Nalpha, Nt), np.NaN)\n",
    "\n",
    "    for k in range(0, Nt):\n",
    "        for j in range(0, Nalpha):\n",
    "            idx = np.where(lp[:, j, k] == 1)[0]\n",
    "            if idx.size != 0:\n",
    "                # First Index\n",
    "                idx = idx[0]\n",
    "                lp_b[j, k] = q[idx]\n",
    "            idx = np.where(lm[:, j, k] == 1)[0]\n",
    "            if idx.size != 0:\n",
    "                # Last Index\n",
    "                idx = idx[-1]\n",
    "                lm_b[j, k] = q[idx]\n",
    "    return lp_b, lm_b\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting functions for visualizing the solutions\n",
    "def plot_post_surface(t, alpha, inventory, q, title):\n",
    "    \"\"\" Generate the plots (b) or (c) of Figure 10.10 on page 270. The surface plot with Inventory as a function\n",
    "    vs alpha and t. You can change the cmap parameter to see more beatiful suface plots, for example, let cmap='coolwarm'.\n",
    "    \"\"\"\n",
    "    x, y = np.meshgrid(alpha, t)\n",
    "    z = np.transpose(inventory)\n",
    "    z[np.isnan(z)] = np.nan\n",
    "    fig = plt.figure()\n",
    "    ax = fig.gca(projection='3d')\n",
    "    surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap='rainbow', vmin=np.nanmin(z),vmax=np.nanmax(z), linewidth=0, antialiased=False)\n",
    "    \n",
    "    #ax.xaxis._axinfo['label']['space_factor'] = 20.0\n",
    "    #ax.yaxis._axinfo['label']['space_factor'] = 20.0\n",
    "    #ax.zaxis._axinfo['label']['space_factor'] = 20.0\n",
    "    \n",
    "    ax.set_xlabel(r'$\\alpha$')\n",
    "    ax.set_ylabel(r'$t$')\n",
    "    ax.set_zlabel(r'$Inventory$')\n",
    "    ax.set_zlim(np.min(q), np.max(q))\n",
    "    plt.title(label=title)\n",
    "    ax.dist=10\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_optimal_posting(alpha, lp_b, lm_b, title):\n",
    "    \"\"\" Plot the figure 10.10 (a) on Page 270. Using the buy and sell postings calculated using the function\n",
    "    generate_sell_buy_posts.\n",
    "    \"\"\"\n",
    "    x = alpha\n",
    "    y1 = lp_b[:, 1]\n",
    "    y2 = lm_b[:, 1]\n",
    "    plt.plot(x, y1, color=\"black\", linewidth=1, linestyle=\"solid\")\n",
    "    plt.plot(x, y2, color=\"black\", linewidth=1, linestyle=\"solid\")\n",
    "    plt.ylim(min(np.nanmin(y1), np.nanmin(y2)), max(np.nanmax(y1), np.nanmax(y2)))\n",
    "    plt.xlim((np.nanmin(alpha), np.nanmax(alpha)))\n",
    "\n",
    "    ax = plt.gca()\n",
    "    area1 = ax.fill_between(x, np.maximum(y1, y2), max(np.nanmax(y1), np.nanmax(y2)), facecolor=\"red\", interpolate=True\n",
    "                            , alpha=0.5)\n",
    "    area2 = ax.fill_between(x, y1, y2, where=(y2 >= y1), facecolor=\"green\", interpolate=True, alpha=0.5)\n",
    "    area3 = ax.fill_between(x, np.minimum(y1, y2), min(np.nanmin(y1), np.nanmin(y2)), facecolor=\"blue\", interpolate=True\n",
    "                            , alpha=0.5)\n",
    "    ax.set_xlabel(r'$\\alpha$')\n",
    "    ax.set_ylabel(r'$Inventory$')\n",
    "    plt.legend([area1, area2, area3], [r'$sell$', r'$buy+sell$', r'$buy$'])\n",
    "    plt.title(label=title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_simulation(Nsims, t, dt, Delta, sigma, lambda_p, lambda_m, short_term_alpha, lm, lp, Nq, varphi):\n",
    "    \"\"\" Generate simulation.\n",
    "    \n",
    "    :param Nsims: number of simulations needed\n",
    "    :param t: a time vector\n",
    "    :param dt: granularity of time vector\n",
    "    :param sigma: a variable\n",
    "    :param lambda_p: a numpy vector\n",
    "    :param lambda_m: a numpy vector\n",
    "    :param short_term_alpha: a short-term-alpha object\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    np.random.seed(20)\n",
    "    \n",
    "    # set up time grid\n",
    "    Nt = t.shape[0]\n",
    "    alpha = short_term_alpha.value\n",
    "    dalpha = short_term_alpha.dalpha\n",
    "\n",
    "    # store the paths of postings, fills, inventory, cash, short-term alpha, and asset price\n",
    "    lp_path = np.full((Nsims, Nt), np.nan)\n",
    "    lm_path = np.full((Nsims, Nt), np.nan)\n",
    "    \n",
    "    isfilled_p = np.full((Nsims, Nt), np.nan)\n",
    "    isfilled_m = np.full((Nsims, Nt), np.nan)\n",
    "    \n",
    "    q_path = np.full((Nsims, Nt), np.nan)\n",
    "    q_path[:, 0] = np.zeros(Nsims)\n",
    "\n",
    "    X = np.full((Nsims, Nt), np.nan)\n",
    "    X[:, 0] = 0\n",
    "\n",
    "    alpha_path = np.full((Nsims, Nt), np.nan)\n",
    "    alpha_path[:, 0] = 0    \n",
    "    \n",
    "    s_path = np.full((Nsims, Nt), np.nan)\n",
    "    s_path[:, 0] = 30\n",
    "\n",
    "    # used for noting if a  market order arrives and if it is buy/sell\n",
    "    isMO = np.zeros((Nsims, Nt), dtype=bool)\n",
    "    buySellMO = np.zeros((Nsims, Nt))\n",
    "\n",
    "\n",
    "    zeta = short_term_alpha.zeta\n",
    "    eta = short_term_alpha.eta\n",
    "    epsilon = short_term_alpha.epsilon\n",
    "\n",
    "    \n",
    "    for i in range(0, Nt - 1, 1):\n",
    "        \n",
    "        # Print Counter\n",
    "        counter(i, Nt - 1)\n",
    "\n",
    "        # Update S by simulating (10.30) forward\n",
    "        s_path[:, i+1] = s_path[:, i] + alpha_path[:, i] * dt + sigma * np.random.randn(Nsims)\n",
    "\n",
    "        # Decide if a Market order arrives\n",
    "        isMO[:, i] = np.random.rand(Nsims) < np.around((1 - np.exp(-dt * (lambda_p + lambda_m))), decimals=4)\n",
    "\n",
    "        # Decide if it is a buy or sell market order\n",
    "        buySellMO[:, i] = 2 * (np.random.rand(Nsims) < lambda_p / (lambda_p + lambda_m)) - 1\n",
    "\n",
    "        # Update alpha by simulating (10.31) forward\n",
    "        alpha_path[:, i+1] = np.exp(-zeta * dt) * alpha_path[:, i] + eta * np.sqrt(dt) * np.random.randn(Nsims) + \\\n",
    "                             epsilon * np.multiply(isMO[:, i], buySellMO[:, i])\n",
    "\n",
    "        tOld = dt * i\n",
    "        tNow = dt * (i + 1)\n",
    "\n",
    "        # find the current limit order volume\n",
    "        idx_t = np.where(t <= tOld)[0][-1]\n",
    "\n",
    "        # perhaps interpolation is better?\n",
    "        alpha_path_idx = np.minimum(alpha.shape[0], np.maximum(1, np.floor((alpha_path[:, i] - alpha[0]) / dalpha) + 1))\n",
    "\n",
    "        # Make decisions of whether to buy or sell\n",
    "        lp_path[:, i] = [lp[int(q_path[j, i] + Nq), int(alpha_path_idx[j]), int(idx_t)] for j in range(0, Nsims)]\n",
    "        lm_path[:, i] = [lm[int(q_path[j, i] + Nq), int(alpha_path_idx[j]), int(idx_t)] for j in range(0, Nsims)]\n",
    "\n",
    "        # Update Inventory\n",
    "        isfilled_p[:, i] = np.multiply((lp_path[:, i]).astype(int), (isMO[:, i]).astype(int), (buySellMO[:, i] == 1).astype(int))\n",
    "        isfilled_m[:, i] = np.multiply((lm_path[:, i]).astype(int), (isMO[:, i]).astype(int), (buySellMO[:, i] == -1).astype(int))\n",
    "\n",
    "        q_path[:, i+1] = q_path[:, i] + isfilled_m[:, i] - isfilled_p[:, i]\n",
    "\n",
    "        X[:, i + 1] = X[:, i] - np.multiply(s_path[:, i] - 0.5 * Delta, isfilled_m[:, i]) + \\\n",
    "                      np.multiply(s_path[:, i] + 0.5 * Delta, isfilled_p[:, i])\n",
    "\n",
    "    # account for terminal liquidating trading actions\n",
    "    X[:, X.shape[1]-1] += np.multiply(s_path[:, s_path.shape[1]-1] +\n",
    "                                      (0.5 * Delta) * np.sign(q_path[:, q_path.shape[1]-1]) +\n",
    "                                      varphi * q_path[:, q_path.shape[1]-1],\n",
    "                                      q_path[:, q_path.shape[1]-1])\n",
    "    \n",
    "    q_path[:, q_path.shape[1]-1] = 0\n",
    "    \n",
    "    print(\"done!\")\n",
    "    \n",
    "    return X, q_path, alpha_path, s_path, lm_path, lp_path, isMO, buySellMO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting functionality for the simulations\n",
    "def plot_simulation(x, y, xlabel, ylabel):\n",
    "    \"\"\" This function can be used to generate the right hand side picture of Figure 10.11 on page 271.\n",
    "    The short-term alpha generated by the simulator vs time. It can also be used to generate inventory size simulated vs time.\n",
    "    \"\"\"\n",
    "    plt.plot(x, y, color=\"blue\", linewidth=1, linestyle=\"solid\")\n",
    "\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.ylabel(ylabel)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_simulation_decision(Delta, t, s_path, lm_path, lp_path, isMO, buySellMO, zoom=None):\n",
    "    \"\"\" This function generates the left hand side picture of Figure 10.11 at page 271, which is the sample path of the \n",
    "    optimal strategy. Green lines show when and at what price the agent is posted. Solid red circles indicates MOs that arrive \n",
    "    and hit/lift the posted bid/offers. Open red circles indicate MOs that arrive and hit/lift the posted bid/offers. \n",
    "    Open red circles indicate MOs that arrive but do not hit/lift the agent's posts. Shaded region is the bid-ask spread.\n",
    "    \"\"\"\n",
    "    s_post_buy = s_path - (0.5 * Delta)\n",
    "    s_post_buy[lm_path == 0] = np.NaN\n",
    "    \n",
    "    is_market_buy_order = np.logical_and(isMO == 1, buySellMO == -1)\n",
    "    # Index for closed circles\n",
    "    idx_closed_buy = np.where(np.logical_and(is_market_buy_order, nan_to_num(lm_path, nan=0) == 1))\n",
    "    # Index for open circles\n",
    "    idx_open_buy = np.where(np.logical_and(is_market_buy_order, nan_to_num(lm_path, nan=1) == 0))\n",
    "\n",
    "    s_post_sell = s_path + (0.5 * Delta)\n",
    "    s_post_sell[lp_path == 0] = np.NaN\n",
    "\n",
    "    is_market_sell_order = np.logical_and(isMO == 1, buySellMO == 1)\n",
    "    # Index for closed circles\n",
    "    idx_closed_sell = np.where(np.logical_and(is_market_sell_order, nan_to_num(lp_path, nan=0) == 1))\n",
    "    # Index for open circles\n",
    "    idx_open_sell = np.where(np.logical_and(is_market_sell_order, nan_to_num(lp_path, nan=1) == 0))\n",
    "\n",
    "    \n",
    "    ax = plt.gca()\n",
    "    ax.plot(t, s_path, color=\"black\", linewidth=2, linestyle=\"solid\")\n",
    "    # Plot the buy side of the LOB\n",
    "    ax.plot(t, s_path - 0.5 * Delta, color=\"black\", linewidth=1, linestyle=\"solid\", zorder=1)\n",
    "    ax.plot(t, s_post_buy, color=\"green\", linewidth=3, linestyle=\"solid\", zorder=2)\n",
    "    # Plot the sell side of LOB\n",
    "    ax.plot(t, s_path + 0.5 * Delta, color=\"black\", linewidth=1, linestyle=\"solid\", zorder=3)\n",
    "    ax.plot(t, s_post_sell, color=\"blue\", linewidth=3, linestyle=\"solid\", zorder=4)\n",
    "    # Plot shaded regions\n",
    "    ax.fill_between(t, s_path + 0.5 * Delta, s_path - 0.5 * Delta, facecolor=\"blue\", alpha=0.2)\n",
    "    # Add closed circles for buy side\n",
    "    ax.scatter(t[idx_closed_buy], s_path[idx_closed_buy] - 0.5 * Delta, color='red' , zorder=5)\n",
    "    # Add open circles for buy side\n",
    "    ax.scatter(t[idx_open_buy], s_path[idx_open_buy] - 0.5 * Delta, facecolors='none', edgecolors='r', zorder=6)\n",
    "    # Add closed circles for sell side\n",
    "    ax.scatter(t[idx_closed_sell], s_path[idx_closed_sell] + 0.5 * Delta, color='red', zorder=7)\n",
    "    # Add open circles for sell side\n",
    "    ax.scatter(t[idx_open_sell], s_path[idx_open_sell] + 0.5 * Delta, facecolors='none', edgecolors='r', zorder=8)\n",
    "    \n",
    "    if zoom is not None:\n",
    "        plt.xlim(zoom)\n",
    "        ylim_min = np.nanmin(s_post_buy[np.logical_and(zoom[0] <= t, t < zoom[1])])\n",
    "        ylim_max = np.nanmax(s_post_sell[np.logical_and(zoom[0] <= t, t < zoom[1])])\n",
    "        plt.ylim((ylim_min, ylim_max))\n",
    "\n",
    "    plt.xlabel(r'$Time$')\n",
    "    plt.ylabel(r'$Midprice$')\n",
    "    plt.show()\n",
    "    \n",
    "def plot_wealth_hist(t, X, q, S, which_path):\n",
    "    \n",
    "    # plot marked to market of wealth sample path\n",
    "    wealth = X[which_path,:] + q[which_path,:] * S[which_path,:]\n",
    "    plt.plot(t, wealth)\n",
    "    plt.xlabel(r'Time ($t$)')\n",
    "    plt.ylabel(r'Marked-to-Market Value')\n",
    "    plt.title(r'Wealth Path')\n",
    "    plt.show()\n",
    "    \n",
    "    # plot histogram of terminal wealth\n",
    "    p = 0.005 # truncate plot to contain 0.005 of data in tails\n",
    "    q= np.quantile(X[:,-1], np.array([p, 1-p]))\n",
    "    h_out = plt.hist(X[:,-1],bins=np.linspace(np.round(q[0]),np.round(q[1]),25))\n",
    "    max_freq = 1.1*np.max(h_out[0])\n",
    "\n",
    "    q = np.quantile(X[:,-1], np.array([0.01, 0.1, 0.5, 0.9, 0.99]))\n",
    "    for i in range(len(q)):\n",
    "        plt.plot(q[i]*np.array([1,1]),np.array([0,max_freq]),'--')\n",
    "    plt.xlabel(r'Terminal Profit & Loss')\n",
    "    plt.ylabel(r'Frequency')\n",
    "    plt.title(r'Terminal Cash Histogram and Quantiles')\n",
    "    plt.legend(('qtl=0.01','qtl=0.1','qtl=0.5','qtl=0.9','qtl=0.99'))\n",
    "    plt.ylim((0,max_freq ))\n",
    "    plt.show()\n",
    "    \n",
    "def plot_turnover_hist(q):\n",
    "    \n",
    "    Y = np.sum(np.abs(np.diff(q,axis=1)), axis=1)\n",
    "    \n",
    "    # plot histogram of terminal wealth\n",
    "    p = 0.005 # truncate plot to contain 0.005 of data in tails\n",
    "    q= np.round(np.quantile(Y, np.array([p, 1-p])))\n",
    "    h_out = plt.hist(Y,bins=np.linspace(0,np.round(q[1]),np.round(q[1])+1))\n",
    "    max_freq = 1.1*np.max(h_out[0])\n",
    "\n",
    "    q = np.quantile(Y, np.array([0.01, 0.1, 0.5, 0.9, 0.99]))\n",
    "    for i in range(len(q)):\n",
    "        plt.plot(q[i]*np.array([1,1]),np.array([0,max_freq]),'--')\n",
    "    plt.xlabel(r'$\\sum_{t=1}^T|q_t-q_{t-1}|$')\n",
    "    plt.ylabel(r'Frequency')\n",
    "    plt.title(r'Total inventory turnover')\n",
    "    plt.legend(('qtl=0.01','qtl=0.1','qtl=0.5','qtl=0.9','qtl=0.99'))\n",
    "    plt.ylim((0,max_freq ))\n",
    "    plt.show()"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}