4_Mispricing_Index_using_KDE_Research.ipynb

{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["![QuantConnect Logo](https://cdn.quantconnect.com/web/i/icon.png)\n", "<hr>"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### KDE Model MI Index Research\n", "\n", "Research the best method for calculating the Mispricing index given a KDE model. \n", "Factors to consider are: \n", "* Computational speed\n", "* "]}, {"cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": ["#----------------------------------------------------------------------------------------------------\n", "# UTIL FUNCTIONS\n", "#----------------------------------------------------------------------------------------------------\n", "import numpy as np\n", "class Parameter():\n", "    def __init__(self,\n", "                resolution,\n", "                start_year, \n", "                start_month=1, \n", "                start_day=1,\n", "                n_years=0,\n", "                n_months=0, \n", "                n_days=0): \n", "        \n", "        self.Resolution = resolution\n", "        \n", "        end_year = start_year + n_years\n", "        end_month = start_month + n_months\n", "        end_day = start_day + n_days \n", "        self.Start = datetime(start_year,start_month,start_day,9,30,0)\n", "        self.End = datetime(end_year, end_month, end_day,16,30,0)\n", "        \n", "class Data():\n", "    def __init__(self, df):\n", "        self._df = df.dropna()\n", "    \n", "    @property \n", "    def prices(self):\n", "        return self._df\n", "    \n", "    @property\n", "    def returns(self):\n", "        return self._df.diff() \n", "    \n", "    @property\n", "    def columns(self):\n", "        return self._df.columns\n", "    \n", "    def plot(self, *args, **kwargs): \n", "        return self._df.plot(*args, **kwargs)\n", "\n", "import random\n", "\n", "def sample_from_bivariate(x_domain, y_domain, weights, n_samples):\n", "    x_domain = x_domain.ravel()\n", "    y_domain = y_domain.ravel()\n", "    weights = np.nan_to_num(weights.ravel()) # account for nans\n", "    \n", "    population = np.array([x_domain, y_domain]).T\n", "    print(population)\n", "    return np.array(random.choices(population, weights, k=n_samples))\n", "\n", "# Utils\n", "def np_remove_nan(x): \n", "    return x[np.logical_not(np.isnan(x).any(axis=1))]\n", "\n", "def np_remove_inf(x):\n", "    return x[np.logical_not(np.isinf(x).any(axis=1))]\n", "\n", "def np_remove_inf_1D(x):\n", "    return x[~np.isinf(x)]\n", "\n", "def np_remove_nan_1D(x):\n", "    return x[~np.isinf(x)]\n", "\n", "# test remove nan \n", "data = np.array([[np.inf,np.inf],\n", "                 [1,2],\n", "                 [1,3],\n", "                 [np.nan, 4],\n", "                 [5,6]])\n", "\n", "data_no_nan = np_remove_nan(data)\n", "data_no_inf = np_remove_inf(data)\n", "data_no_non_num = np_remove_nan(np_remove_inf(data))\n", "\n", "check_for_nan = lambda x: True not in np.isnan(x)\n", "check_for_inf = lambda x: True not in np.isinf(x)\n", "\n", "assert check_for_nan(data_no_nan)\n", "assert check_for_inf(data_no_inf)\n", "assert check_for_nan(data_no_non_num) and check_for_inf(data_no_non_num)\n", "assert data_no_non_num.shape == (3,2)\n"]}, {"cell_type": "code", "execution_count": 39, "metadata": {"tags": []}, "outputs": [], "source": ["#----------------------------------------------------------------------------------------------------\n", "# LOAD DATA\n", "#----------------------------------------------------------------------------------------------------\n", "\n", "print(\"Setting parameters...\")\n", "# parameters\n", "START_YEAR = 2011\n", "START_MONTH = 1\n", "START_DAY = 1\n", "N_YEARS = 5\n", "N_MONTHS = 0\n", "N_DAYS = 0\n", "RESOLUTION = Resolution.Daily\n", "PARAM = Parameter(RESOLUTION, START_YEAR, START_MONTH, START_DAY,\n", "                  N_YEARS, N_MONTHS, N_DAYS)\n", "print(\"Loading QuantBook...\")\n", "qb = QuantBook()\n", "\n", "# Specify list of correlated tickers for S&P 500 \n", "tickers = [\"SPY\",\"XLK\", \"VGT\", \"IYW\", \"IGV\"]\n", "\n", "# register tickers to quantbook\n", "for ticker in tickers: \n", "    qb.AddEquity(ticker)\n", "    \n", "print(\"Loading historical data...\")\n", "# get historical prices\n", "history = qb.History(tickers, PARAM.Start, PARAM.End, PARAM.Resolution)\n", "\n", "print(\"Preparing historical data...\")\n", "# Unpack dataframe\n", "Open = history['open'].unstack(level=0)\n", "Close = history['close'].unstack(level=0)\n", "High = history['high'].unstack(level=0)\n", "Low = history['low'].unstack(level=0)\n", "\n", "# create data object\n", "close = Data(Close)\n", "print(\"Done\")"]}, {"cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": ["close.prices.plot()"]}, {"cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": ["#----------------------------------------------------------------------------------------------------\n", "# ESTIMATE EMPIRICAL MARGINAL DISTRIBUTIONS AND CALCULATE MARGINAL VALUES\n", "#----------------------------------------------------------------------------------------------------\n", "\n", "from statsmodels.distributions.empirical_distribution import ECDF\n", "\n", "# Empirically estimate marginal distributions\n", "marginal_dist = {}\n", "for asset in close.columns:\n", "    marginal_dist[asset] = ECDF(close.returns[asset])\n", "\n", "marginal_values = pd.DataFrame()\n", "for asset in close.columns: \n", "    marginal_values[asset] = close.returns[asset].apply(marginal_dist[asset])\n", "    "]}, {"cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": ["np.array(close.returns).shape"]}, {"cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": ["\n", "def get_max_corr_pair(corr_matrix): \n", "    sol = (corr.where(np.triu(np.ones(corr.shape), k=1).astype(np.bool))\n", "                 .stack()\n", "                 .sort_values(ascending=False))\n", "    \n", "    col1 = close.columns[sol.index[sol == max(sol)].labels[0]].values[0]\n", "    col2 = close.columns[sol.index[sol == max(sol)].labels[1]].values[0]\n", "    return(col1, col2)\n", "\n", "#----------------------------------------------------------------------------------------------------\n", "# CALCULATE KENDALL TAU CORRELATION\n", "#----------------------------------------------------------------------------------------------------\n", "corr = close.returns.corr(method='kendall')\n", "corr.style.background_gradient(cmap='viridis')\n", "\n", "pair_override = (\"SPY R735QTJ8XC9X\", \"XLK RGRPZX100F39\") # tuple pair\n", "pair = get_max_corr_pair(corr)\n", "pair = pair_override if  pair_override is not None else pair \n", "print(\"Pair with maximum correlation: {}\".format(pair))\n", "\n", "marginals = np.array(marginal_values[list(pair)].values)"]}, {"cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [], "source": ["corr"]}, {"cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": ["\n", "\n", "from time import time\n", "from scipy.stats import norm\n", "import statsmodels.api as sm\n", "\n", "#----------------------------------------------------------------------------------------------------\n", "# FIT KDE MODEL TO DATA MAPPED USING GAUSSIAN TRANSFORMATION\n", "#----------------------------------------------------------------------------------------------------\n", "# transform samples to R^2 domain\n", "marginals_R_raw = norm.ppf(marginals)\n", "marginals_R = np_remove_inf(marginals_R_raw)\n", "\n", "# fit a kernel density estimate \n", "_BANDWIDTH = 'cv_ml' #'normal_reference' # 'cv_ml', 'cv_ls'\n", "kernel = sm.nonparametric.KDEMultivariate\n", "t = time()\n", "model = kernel(data=marginals_R, var_type='cc', bw=_BANDWIDTH)\n", "elapsed = time() - t \n", "\n", "print(\"Model fit in {} seconds\".format(elapsed))\n", "print(\"Model bandwidth {}\".format(model.bw))\n"]}, {"cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": ["plt.figure(figsize=(15,15))\n", "plt.scatter(marginals[:,0], marginals[:,1], s=0.3)\n", "\n", "plt.figure(figsize=(15,15))\n", "plt.scatter(marginals_R[:,0], marginals_R[:,1])"]}, {"cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [], "source": ["# define domain in [0,1]^2 \n", "_equidistant_grid = False\n", "if _equidistant_grid: \n", "    _min = 0.01 # offset to avoid inf\n", "    _max = 1\n", "    w = 0.001 \n", "    x_valuesU = np.arange(_min, _max, w)\n", "    y_valuesU = np.arange(_min, _max, w)\n", "else: \n", "    # by transforming the uniformly distributed grid using the normal cdf, \n", "    # the grid density is thicker near the extremes. \n", "    _min = -3; _max = 3; w = 0.01\n", "    x_valuesU = norm.cdf(np.arange(_min, _max, w))\n", "    y_valuesU = norm.cdf(np.arange(_min, _max, w))\n", "    \n", "# transform to R^2 domain and make meshgrid\n", "x_hat_valuesR = norm.ppf(x_valuesU)\n", "y_hat_valuesR = norm.ppf(y_valuesU)\n", "x_hat_meshR, y_hat_meshR = np.meshgrid(x_hat_valuesR, y_hat_valuesR)\n", "\n", "# make predictions in R^2 domain and make meshgrid\n", "data_predictR = pd.DataFrame({\"Xhat\":x_hat_meshR.ravel(), \"Yhat\":y_hat_meshR.ravel()})\n"]}, {"cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": ["import time\n", "# sample from KDE model to fit interpolation model\n", "t = time.time()\n", "USE_CDF = False\n", "if USE_CDF: \n", "    z_hat_meshR = model.cdf(data_predict=data_predictR.to_numpy())\n", "else: \n", "    z_hat_meshR = model.pdf(data_predict=data_predictR.to_numpy())\n", "    \n", "elapsed = time.time() - t \n", "print(\"Sampling took {} seconds\".format(elapsed))\n"]}, {"cell_type": "code", "execution_count": 95, "metadata": {}, "outputs": [], "source": ["# Calculating pdf for a chunk of data with n_worker = 4 takes around 14 seconds\n", "t = time.time() \n", "z_hat_meshR = model.pdf(data[0:int(len(data)/n_worker)])\n", "elapsed = time.time() - t \n", "print(f\"Sampling took {elapsed} seconds\")\n", "print(z_hat_meshR)"]}, {"cell_type": "code", "execution_count": 121, "metadata": {}, "outputs": [], "source": ["'''\n", "The overhead with creating a pool is too great to be advantageous apparently - we get at least \n", "double the computation time. \n", "'''\n", "i = 1000\n", "n_worker = 4\n", "testdata = data[0:i,:]\n", "_chunksize = int(len(testdata)/n_worker)\n", "_chunksize = _chunksize if _chunksize > 1 else 1\n", "print(f\"num workers: {n_worker}\")\n", "print(f\"Chunksize: {_chunksize}\")\n", "\n", "t = time.time()\n", "res = model.pdf(testdata)\n", "elapsed = time.time() - t\n", "print(f\"basic func took {elapsed} seconds\")\n", "print(res)\n", "\n", "with mp.Pool(n_worker) as pool:\n", "    t = time.time()\n", "    res = pool.map(model.pdf, testdata, chunksize=500)\n", "    elapsed = time.time() - t \n", "    print(f\"pooled basic_func {elapsed} seconds\")\n", "    print(res)"]}, {"cell_type": "code", "execution_count": 145, "metadata": {}, "outputs": [], "source": ["'''\n", "Pretty much the same result with this method. Might be due to overhead with serialisation of the KDE object \n", "itself + the data. In any ways the best way to speed up this process seems to be by writing an optimised KDE \n", "library myself. \n", "'''\n", "# parameters \n", "i = 1000000\n", "i = min(len(data), i)\n", "n_worker = 1\n", "testdata = data[0:i,:]\n", "_chunksize = int(len(testdata)/n_worker)\n", "_remainder = len(testdata) % n_worker\n", "print(f\"num data: {i}\")\n", "print(f\"Chunk size: {_chunksize}\")\n", "print(f\"_remainder: {_remainder}\")\n", "\n", "# single process\n", "t = time.time()\n", "res = model.pdf(testdata)\n", "elapsed = time.time() - t\n", "print(f\"basic func took {elapsed} seconds\")\n", "# print(res)\n", "\n", "# --- multiprocess ---\n", "# split the data set into chunks\n", "chunklist = []\n", "zeroindex_offset = 0\n", "for i in range(n_worker):\n", "    prev_end = 0\n", "    start_index = _chunksize*i\n", "    rem = 0 if i < n_worker - zeroindex_offset else _remainder\n", "    end_index = start_index + _chunksize - zeroindex_offset + rem\n", "    print(f\"chunk range: {start_index} to {end_index}\")\n", "    chunk = testdata[start_index:end_index,:]\n", "    chunklist.append(chunk)    \n", "\n", "with mp.Pool(n_worker) as pool:\n", "    t = time.time()\n", "    res = pool.map(model.pdf, chunklist)\n", "    elapsed = time.time() - t \n", "    print(f\"pooled basic_func {elapsed} seconds\")\n", "    "]}, {"cell_type": "code", "execution_count": 144, "metadata": {}, "outputs": [], "source": ["import sys\n", "\n", "def get_size(obj, seen=None):\n", "    \"\"\"Recursively finds size of objects\"\"\"\n", "    size = sys.getsizeof(obj)\n", "    if seen is None:\n", "        seen = set()\n", "    obj_id = id(obj)\n", "    if obj_id in seen:\n", "        return 0\n", "    # Important mark as seen *before* entering recursion to gracefully handle\n", "    # self-referential objects\n", "    seen.add(obj_id)\n", "    if isinstance(obj, dict):\n", "        size += sum([get_size(v, seen) for v in obj.values()])\n", "        size += sum([get_size(k, seen) for k in obj.keys()])\n", "    elif hasattr(obj, '__dict__'):\n", "        size += get_size(obj.__dict__, seen)\n", "    elif hasattr(obj, '__iter__') and not isinstance(obj, (str, bytes, bytearray)):\n", "        size += sum([get_size(i, seen) for i in obj])\n", "    return size\n", "\n", "get_size(model)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["t = time()\n", "# transform predictions in R^2 domain back to [0,1]^2 domain \n", "z_valuesU = np.array(z_hat_meshR / (norm.pdf(data_predictR.Xhat) * norm.pdf(data_predictR.Yhat)))\n", "\n", "# make meshgrid in [0,1]^2 domain for plotting\n", "x_meshU, y_meshU = np.meshgrid(x_valuesU, y_valuesU)\n", "z_meshU = z_valuesU.reshape(x_hat_meshR.shape)\n", "elapsed2 = time() - t"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# contour plot of resulting kde esimate of copula pdf in [0,1]^2\n", "from matplotlib import cm\n", "plt.figure(figsize=(12,10))\n", "plt.contourf(x_meshU, y_meshU, z_meshU, 1000, cmap=cm.coolwarm)\n", "plt.colorbar()\n", "\n", "# 3d plot of resulting kde estimate of copula pdf in [0,1]^2\n", "from mpl_toolkits import mplot3d\n", "\n", "fig = plt.figure(figsize=(10,10))\n", "ax = fig.gca(projection='3d')\n", "surf = ax.plot_surface(x_meshU, y_meshU, z_meshU,\n", "                       linewidth=0, antialiased=False, cmap=cm.coolwarm)\n", "ax.view_init(45, 135)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We have so far obtained a grid of pdf values using the Gaussian Transformation method. We now use interpolation between these points to obtain a continuous approximation of the copula pdf. The interpolation method has a integral function implemented so this can be used to calculate cdf values and conditional cdf values. "]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import scipy\n", "\n", "t = time()\n", "\n", "# Use degree 1 for interpolation as larger degrees can overshoot to lead to negative pdf values\n", "# do not enforce bbox as it causes the interpolation to overshoot into negative values _BBOX = [0,1,0,1]\n", "interp_model = scipy.interpolate.RectBivariateSpline(x_valuesU, y_valuesU, z_meshU, kx=1, ky=1)\n", "elapsed = time() - t\n", "print(\"Interpolation model fit in {}s\".format(elapsed))"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# ---------------------------------------------------------------------------------\n", "# plot the interpolated model\n", "# ---------------------------------------------------------------------------------\n", "from matplotlib import cm\n", "\n", "z_mesh_interp = interp_model(x_valuesU, y_valuesU)\n", "z_mesh_interp = z_mesh_interp\n", "\n", "plt.figure(figsize=(12,10))\n", "plt.contourf(x_meshU, y_meshU, z_mesh_interp, 300, cmap=cm.coolwarm)\n", "plt.colorbar()\n", "\n", "# 3d plot of resulting kde estimate of copula pdf in [0,1]^2\n", "from mpl_toolkits import mplot3d\n", "\n", "fig = plt.figure(figsize=(20,20))\n", "ax = fig.gca(projection='3d')\n", "surf = ax.plot_surface(x_meshU, y_meshU, z_mesh_interp,\n", "                       linewidth=0, antialiased=False, cmap=cm.coolwarm)\n", "ax.view_init(45, 135)\n", "\n", "print(\"SSE: {}\".format(interp_model.get_residual()))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Above we investigated an alternative approach: the original method used the KDE + Interpolation to estimate the copula pdf. From this we would exploit the interpolation model's ease of integral calculation to find the cdf and the conditional copula values. \n", "\n", "The alternative to this was to use the KDE + Interpolation to estimate the cdf directly by using the KDE cdf member function. This allows us to avoid having to integrate the pdf to find the cdf and the conditional which requires extra compuational step for normalisation. \n", "\n", "However, it has been shown above that the latter approach is unable to capture the lower tail dependencies well enough. \n", "\n", "We therefore resort to using the former method. "]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# test that the cdf gives us the pdf\n", "if USE_CDF: \n", "    z_values = interp_model.ev(x_meshU, y_meshU, dx=1, dy=1)\n", "    z_values_norm = z_values / np.max(z_values)\n", "    z_values_norm = z_values_norm.reshape(x_meshU.shape)\n", "    plt.figure(figsize=(12,10))\n", "    plt.contourf(x_meshU, y_meshU, z_values_norm, 500, cmap=cm.coolwarm)\n", "    plt.colorbar()\n", "\n", "    # 3d plot of resulting kde estimate of copula pdf in [0,1]^2\n", "    from mpl_toolkits import mplot3d\n", "\n", "    fig = plt.figure(figsize=(20,20))\n", "    ax = fig.gca(projection='3d')\n", "    surf = ax.plot_surface(x_meshU, y_meshU, z_values_norm,\n", "                           linewidth=0, antialiased=False, cmap=cm.coolwarm)\n", "    ax.view_init(45, 135)\n", "\n", "    print(\"SSE: {}\".format(interp_model.get_residual()))\n", "    "]}, {"cell_type": "markdown", "metadata": {}, "source": ["We now consider how to compute the mispricing index efficiently. From the definition of the copula density, note that:  \n", "\n", "$$\n", "    c(u,v) = \n", "    \\frac{ \\partial^2 C(u,v) }{ \\partial u \\partial v } = \n", "    \\frac{\\partial}{\\partial u}\\bigg( \\frac{\\partial C(u,v)}{\\partial v} \\bigg) = \n", "    \\frac{\\partial}{\\partial u}\\bigg( C(u|v) \\bigg)\n", "$$\n", "\n", "Therefore,\n", "\n", "$$\n", "   \\therefore C(u|v) = \\int_0^{u} c(s,v)ds\n", "$$\n"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# Calculate MI for u |v first\n", "from scipy import integrate\n", "\n", "# let's say u = 0.9, v = 0.9\n", "u_val = 0.000795\n", "v_val = 0.000795\n", "\n", "d = 0.0001\n", "rng = np.arange(0,u_val,d)\n", "fullrng = np.arange(0,1,d)\n", "\n", "# C(u=0.9|v=0.9)\n", "c_val_for_norm = interp_model.ev(fullrng, v_val)\n", "c_values = interp_model.ev(rng, v_val)\n", "\n", "c = np.trapz(c_values, rng)\n", "norm_const = np.trapz(c_val_for_norm, fullrng)\n", "integral = c / norm_const\n", "\n", "print(\"Non-normalised sum: {}\".format(c))\n", "print(\"Normalisation const: {}\".format(norm_const))\n", "print(\"Normalised sum: {}\".format(integral))\n", "\n", "# for plotting\n", "pdf_norm_domain = c_values / norm_const\n", "pdf_full_norm_domain = c_val_for_norm / norm_const\n", "cdf_full_norm_domain = integrate.cumtrapz(pdf_full_norm_domain, fullrng) \n", "cdf_norm_domain = integrate.cumtrapz(pdf_norm_domain, rng)\n", "\n", "fig,ax = plt.subplots(2,1, figsize=(10,10))\n", "ax[0].plot(fullrng, pdf_full_norm_domain)\n", "ax[0].plot(rng, pdf_norm_domain)\n", "ax[1].plot(fullrng[1:], cdf_full_norm_domain)\n", "ax[1].plot(rng[1:], cdf_norm_domain)\n", "ax[0].grid()\n", "ax[1].grid()\n", "ax[0].set_title(\"pdf\")\n", "ax[1].set_title(\"cdf\")\n", "\n", "# Check for violation\n", "_TOL= 1e-5\n", "_max_cdf = np.max(cdf_full_norm_domain)\n", "_min_pdf = np.min(pdf_norm_domain)\n", "assert _max_cdf <= 1.0 + _TOL, \"Violation: CDF exceeds 1 with value {}\".format(_max_cdf)\n", "assert integral <= 1.0 + _TOL, \"Violation: Conditional exceeds 1 with value {}\".format(integral)\n", "assert _min_pdf >= 0.0, \"Violation: negative pdf with value {}\".format(_min_cdf)\n"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [], "source": ["# now do it to account for both \n", "# integrate for particular value of u or v\n", "from scipy import integrate\n", "\n", "def _get_range_linear(val=None, width=1e-3):\n", "    '''\n", "        Using this is not recommended as any u,v values with resolution\n", "        smaller than 'width' will not yield sensible results. \n", "        Args: \n", "            val: The range upper limit. If None, defaults to 1. \n", "        Returns: \n", "            Numpy array within range [0,val] with spacing = width\n", "    '''\n", "    _MINRNG = 0\n", "    _MAXRNG = 1 if val is None else val\n", "    return np.arange(_MINRNG, _MAXRNG, width)\n", "\n", "def _get_range_gaussian_transformation(val=None, width=1e-3): \n", "    '''\n", "        Using this is recommended as opposed to the _get_range_linear as\n", "        this has higher sample density in the extremes and therefore avoids\n", "        the issues with results where val is smaller than the 'width' value. \n", "        Args: \n", "            val: The range upper limit. If None, defaults to 1. \n", "        Returns: \n", "            Numpy array within range [0,val] with spacing = width\n", "    '''\n", "    _MINRNG = -4\n", "    _MAXRNG_default = 4\n", "    # take account for None and inf values\n", "    cond = ((val is None) or (norm.ppf(val) > _MAXRNG_default))\n", "    _MAXRNG = _MAXRNG_default if cond else norm.ppf(val)\n", "    linrng = np.arange(_MINRNG, _MAXRNG, width)\n", "    return norm.cdf(linrng)\n", "\n", "def _integrate(yval, xval): \n", "    '''Wrapper for numpy numerical integration via trapezoidal rule'''\n", "    return np.trapz(yval, xval)\n", "\n", "def MI(u, v, res=1e-3, model=None, debug=False): \n", "    '''\n", "        Calculate the mispricing index for a given u,v value\n", "        Args: \n", "            u: Scalar bounded in range [0,1]. \n", "            v: Scalar bounded in range [0,1]. \n", "            res: width for numerical integration. Larger value \n", "            leads to faster computation for less accuracy\n", "    '''\n", "    _FULLRNG = _get_range_gaussian_transformation(width=res)\n", "    _urng = _get_range_gaussian_transformation(u, width=res)\n", "    _vrng = _get_range_gaussian_transformation(v, width=res)\n", "    \n", "    # c(u,v=v') c(u=u',v)\n", "    c_uv, c_vu = model.ev(_urng, v), model.ev(u, _vrng)\n", "    \n", "    # copula values for normalisation const\n", "    z_uv, z_vu = model.ev(_FULLRNG, v), model.ev(u, _FULLRNG)\n", "    \n", "    # integrate \n", "    C_uv, C_vu = _integrate(c_uv, _urng), _integrate(c_vu, _vrng)\n", "    Z_uv, Z_vu = _integrate(z_uv, _FULLRNG), _integrate(z_vu, _FULLRNG)\n", "    \n", "    if debug: \n", "        print(\"non-normalised values: \\n C(u|v): {}, C(v|u): {}\".format(C_uv, C_vu))\n", "        print(\"normalisation const: \\n Z(u|v): {}, Z(v|u): {}\".format(Z_uv, Z_vu))\n", "        print(\"normalised values: \\n C(u|v): {}, C(v|u): {}\".format(C_uv/Z_uv, C_vu/Z_vu))\n", "    \n", "    # MI values \n", "    return C_uv/Z_uv, C_vu/Z_vu\n", "\n", "# test MI function\n", "u = 0.000795\n", "v = 0.9\n", "mi_uv, mi_vu = MI(u,v,res=1e-3,model=interp_model, debug=True)\n", "\n", "print(\"Mispricing Index C(u|v): {}\".format(mi_uv))\n", "print(\"Mispricing Index C(v|u): {}\".format(mi_vu))\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We now have a method for calculating the mispricing index. Let's verify that this generates valid signals. "]}, {"cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [], "source": ["def MI_series(x, model, pair):\n", "    return MI(x[pair[0]], x[pair[1]], model=model)\n", "\n", "def signaltype_cond(x, siglvl): \n", "    '''\n", "    Long signal is Buy U and Sell V (+(U-V))\n", "    Short signal is Sell U and Buy V (-(U-V))\n", "    '''\n", "    if (x['C_uv'] > UPPER and x['C_vu'] < LOWER): \n", "        return 1\n", "    elif  (x['C_uv'] < LOWER and x['C_vu'] > UPPER): \n", "        return -1 \n", "    else: \n", "        return 0\n", "\n", "test_marginals = marginal_values[list(pair)]\n", "# calculate MI values\n", "df_MI_values = \\\n", "    test_marginals.apply(MI_series, args=(interp_model, pair), axis=1).apply(pd.Series)\n", "df_MI_values.columns = ['C_uv', 'C_vu']\n", "\n", "# identify rows where the MI surpasses the threshold\n", "siglvl = 0.1\n", "UPPER, LOWER = 1 - siglvl, siglvl\n", "df_MI_values['signal'] = \\\n", "    df_MI_values.apply(\n", "        lambda x: ((x['C_uv'] > UPPER and x['C_vu'] < LOWER) or \\\n", "                   (x['C_uv'] < LOWER and x['C_vu'] > UPPER)), \n", "        axis=1)\n", "\n", "df_MI_values['cumulative_signal'] = \\\n", "    df_MI_values.apply(signaltype_cond, args=(siglvl,), axis=1).cumsum()\n", "\n", "# join to main df for comparison\n", "prices = close.prices / close.prices.iloc[0]\n", "df = prices[list(pair)].join(df_MI_values, how='left')\n", "df_signal = df[df.signal==True]\n", "\n", "# Calculate pnl\n", "_TRADE_SIZE = 100 \n", "df['pnl'] = (df.cumulative_signal * (df[pair[0]] - df[pair[1]]) * _TRADE_SIZE).cumsum()\n"]}, {"cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [], "source": ["_MARKER = 'o'\n", "_MARKER_SIZE = 2.5\n", "_XLIM = [datetime(2011,6,1), datetime(2012,1,1)]\n", "\n", "# df = df[(df.index > _XLIM[0]) & (df.index < _XLIM[1])]\n", "\n", "fig, ax = plt.subplots(3, figsize=(20,40))\n", "col2drop = [\"C_uv\", \"C_vu\", \"signal\", \"cumulative_signal\", \"pnl\"]\n", "df.drop(labels=col2drop, axis=1).plot(ax=ax[0], marker=_MARKER, ms=_MARKER_SIZE)\n", "df.drop(labels=col2drop, axis=1).diff().plot(ax=ax[1], marker=_MARKER, ms=_MARKER_SIZE)\n", "\n", "returns = df.drop(labels=col2drop, axis=1).diff()\n", "returns_signal = returns[df.signal == True].join(df[[\"C_uv\", \"C_vu\"]], how='left')\n", "\n", "diff = pd.DataFrame(returns[pair[0]] - returns[pair[1]])\n", "diff.columns = [\"diff\"]\n", "ax[2].vlines(diff.index.values, [0], diff.values)\n", "ax[1].vlines(diff.index.values, [0], diff.values)\n", "_YLIM = [-0.015,0.015]\n", "ax[2].set_ylim(_YLIM)\n", "\n", "for row in returns_signal.iterrows():\n", "    index, data = row\n", "    ax[1].annotate('x<-'+str(round(data[\"C_uv\"],2)), xy=(index, data[pair[0]]))\n", "    ax[1].annotate('x<-'+str(round(data[\"C_vu\"],2)), xy=(index, data[pair[1]]))\n", "\n", "diff_signal = diff[df.signal==True].join(df[[\"C_uv\", \"C_vu\"]], how='left')\n", "for row in diff_signal.iterrows(): \n", "    index, data = row\n", "    d = (round(data[\"C_uv\"],2), round(data[\"C_vu\"],2))\n", "    ax[2].annotate('x<-' + str(d), xy=(index, data[\"diff\"]))\n", "    ax[1].annotate('x<-' + str(d), xy=(index, data[\"diff\"]))\n", "\n", "_offset = 0.05\n", "for row in df_signal.iterrows(): \n", "    index, data = row\n", "    ax[0].vlines(index, data[pair[1]]-_offset, data[pair[1]] + _offset)\n", "    ax[0].vlines(index, data[pair[0]]-_offset, data[pair[0]] + _offset)\n", "\n", "# major_ticks = \n", "# for axis in ax: \n", "#     axis.grid()\n", "\n", "#     axis.set_xticks(major_ticks)\n", "#     axis.set_xticks(minor_ticks, minor=True)\n", "#     axis.set_yticks(major_ticks)\n", "#     axis.set_yticks(minor_ticks, minor=True)"]}, {"cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [], "source": ["# reset df_ object\n", "try: \n", "    del df_\n", "except NameError as e: \n", "    pass \n", "\n", "fig, ax = plt.subplots(figsize=(20,20))\n", "df_ = df.join(diff, how='left')\n", "ax.vlines(df_.index.values, [0], df_['diff'].values)\n", "for row in df_.iterrows(): \n", "    index, data = row\n", "    d = (round(data[\"C_uv\"],2), round(data[\"C_vu\"],2))\n", "    ax.annotate('x<-' + str(d), xy=(index, data[\"diff\"]))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Take a look at cumulative MI signal"]}, {"cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [], "source": ["_MARKER = 'o'\n", "_MARKER_SIZE = 2.5\n", "_XLIM = [datetime(2011,6,1), datetime(2012,1,1)]\n", "\n", "# limit the data to specified range in _XLIM\n", "# df = df[(df.index > _XLIM[0]) & (df.index < _XLIM[1])]\n", "\n", "\n", "fig, ax = plt.subplots(4, figsize=(20,40))\n", "col2drop = [\"C_uv\", \"C_vu\", \"signal\", \"cumulative_signal\", \"pnl\"]\n", "\n", "# First plot - price series\n", "df.drop(labels=col2drop, axis=1).plot(ax=ax[0], marker=_MARKER, ms=_MARKER_SIZE)\n", "\n", "# Second plot - Returns difference plot\n", "returns = df.drop(labels=col2drop, axis=1).diff()\n", "returns_signal = returns[df.signal == True].join(df[[\"C_uv\", \"C_vu\"]], how='left')\n", "diff = pd.DataFrame(returns[pair[0]] - returns[pair[1]])\n", "diff.columns = [\"diff\"]\n", "ax[1].vlines(diff.index.values, [0], diff.values)\n", "_YLIM = [-0.015,0.015]\n", "ax[1].set_ylim(_YLIM)\n", "\n", "# Third plot - cumulative signal\n", "df.cumulative_signal.plot(ax=ax[2], marker=_MARKER, ms=_MARKER_SIZE)\n", "\n", "# Fourth plot - estimated Pnl value (excluding trade commission)\n", "df.pnl.plot(ax=ax[3], marker=_MARKER, ms=_MARKER_SIZE)\n", "ax[3].legend(['Pairs Trading', 'Benchmark'])\n", "\n", "# annotations \n", "diff_signal = diff[df.signal==True].join(df[[\"C_uv\", \"C_vu\"]], how='left')\n", "for row in diff_signal.iterrows(): \n", "    index, data = row\n", "    d = (round(data[\"C_uv\"],2), round(data[\"C_vu\"],2))\n", "    ax[1].annotate('x<-' + str(d), xy=(index, data[\"diff\"]))\n", "\n", "_offset = 0.1\n", "for row in df_signal.iterrows(): \n", "    index, data = row\n", "    ax[0].vlines(index, data[pair[1]]-_offset, data[pair[1]] + _offset)\n", "    ax[0].vlines(index, data[pair[0]]-_offset, data[pair[0]] + _offset)\n"]}, {"cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [], "source": ["Roll_Max = df.pnl.cummax()\n", "Daily_Drawdown = df.pnl/Roll_Max - 1.0\n", "Max_Daily_Drawdown = Daily_Drawdown.cummin()\n", "print(f\"Max Daily Drawdown: {Max_Daily_Drawdown.iloc[-1]}%\")\n", "\n", "daily_return = df.pnl.pct_change(1)\n", "sharpe_ratio = 255**0.5*daily_return.iloc[2:].mean() / daily_return.iloc[2:].std()\n", "print(f\"Sharpe ratio: {sharpe_ratio}\")"]}, {"cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [], "source": ["daily_return.iloc[2:]"]}, {"cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [], "source": ["plt.plot(close.prices['SPY'] / close.prices['SPY'].iloc[0])"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}], "metadata": {"celltoolbar": "Slideshow", "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.6.8"}}, "nbformat": 4, "nbformat_minor": 2}