{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Independent and identically distributed returns\n\nA very common assumption in modelling of financial returns is independence in time.\nIn other words, there is no use in considering today (or any other past day) returns to predict\ntomorrow returns. Another very popular assumption is that returns have the same distribution.\nCombining both of these assumptions we basically claim that there is some underlying multivariate\ndistribution and the actual returns we observe are just independent samples from it.\n\n:code:`deepdow` can be applied in a market with the above described IID dynamics. This example will\nanalyze how good :code:`deepdow` is at learning optimal portfolios given these dynamics.\n\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>:code:`deepdow` contains many tools that help extract hidden features from the original ones. However,\n    under the IID setup we are in a position where the lagged returns have no predictive value. We\n    have zero features which is very unusual from the point of view of standard machine learning\n    applications. However, nothing prevents us from learning a set of constants that uniquely determine\n    the final allocation. As in the more standard case with features, these parameters\n    can be learned via minimization of an empirical loss like mean returns, sharpe ratio, etc.</p></div>\n\nThis example is divided into the following sections\n\n1. Creation of synthetic return series\n2. Convex optimization given the ground truth\n3. Sample estimators of covmat and mean\n4. Implementation and setup of a :code:`deepdow` network\n5. Training and evaluation\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Preliminaries\nLet us start with importing all important dependencies.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import cvxpy as cp\nfrom deepdow.benchmarks import Benchmark\nfrom deepdow.callbacks import EarlyStoppingCallback\nfrom deepdow.data import InRAMDataset, RigidDataLoader\nfrom deepdow.layers import NumericalMarkowitz\nfrom deepdow.losses import MeanReturns, StandardDeviation\nfrom deepdow.experiments import Run\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport torch"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let us also set seeds for both :code:`numpy` and :code:`torch`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "torch.manual_seed(4)\nnp.random.seed(21)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Creation of synthetic return series\nTo ensure that we generate realistic samples we are going to use S&P500 stocks. Namely, we precomputed\ndaily mean returns and covariance matrix for a subset of all the stocks and one can just simply\nload them.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mean_all = pd.read_csv('sp500_mean.csv', index_col=0, header=None).iloc[:, 0]\ncovmat_all = pd.read_csv('sp500_covmat.csv', index_col=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let us now randomly select some assets\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_assets = 15\n\nasset_ixs = np.random.choice(len(mean_all), replace=False, size=n_assets)\nmean = mean_all.iloc[asset_ixs]\ncovmat = covmat_all.iloc[asset_ixs, asset_ixs]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "These are the expected returns\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sns.heatmap(mean.to_frame())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And this is the covariance matrix\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sns.heatmap(covmat)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Additionally let us inspect them in the risk and return plane\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "df_risk_ret = pd.DataFrame({'risk': np.diag(covmat) ** (1 / 2),\n                            'return': mean.values})\n\nx_lim = (0, df_risk_ret['risk'].max() * 1.1)\ny_lim = (df_risk_ret['return'].min() * 1.1, df_risk_ret['return'].max() * 1.1)\nax = df_risk_ret.plot.scatter(x='risk', y='return')\nax.set_xlim(x_lim)\nax.set_ylim(y_lim)\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From now on we see the mean and covmat as the ground truth parameters that fully determine the\nmarket dynamics. We now generate some samples from a multivariate normal distribution with\nmean and covmat as the defining parameters.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_timesteps = 1500\nreturns = pd.DataFrame(np.random.multivariate_normal(mean.values,\n                                                     covmat.values,\n                                                     size=n_timesteps,\n                                                     check_valid='raise'), columns=mean.index)\n\nsns.violinplot(data=returns.iloc[:, ::3])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Convex optimization given the ground truth\nSince we posses the ground truth mean and covmat we can find optimal weight allocation given\nsome **objective** and **constraints**. We are going to consider the following objectives\n\n- Minimum variance\n- Maximum return\n- Maximum utility\n\nSince the combination of our objective and constraints will result in a convex optimization problem\nwe can be sure that the optimal solution exists and is unique. Additionally, we will use\n:code:`cvxpy` to find the solution numerically. For simplicity, we will refer to all of these\noptimization problems as *Markowitz optimization problems*. See below the implementation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class Markowitz:\n    \"\"\"Solutions to markowitz optimization problems.\n\n    Parameters\n    ----------\n    mean : np.ndarray\n        1D array representing the mean of returns. Has shape `(n_assets,)`.\n\n    covmat : np.ndarray\n        2D array representing the covariance matrix of returns. Has shape `(n_assets, n_assets)`.\n\n    \"\"\"\n\n    def __init__(self, mean, covmat):\n        if mean.ndim != 1 and covmat.ndim != 2:\n            raise ValueError('mean needs to be 1D and covmat 2D.')\n\n        if not (mean.shape[0] == covmat.shape[0] == covmat.shape[1]):\n            raise ValueError('Mean and covmat need to have the same number of assets.')\n\n        self.mean = mean\n        self.covmat = covmat\n\n        self.n_assets = self.mean.shape[0]\n\n    def minvar(self, max_weight=1.):\n        \"\"\"Compute minimum variance portfolio.\"\"\"\n        w = cp.Variable(self.n_assets)\n        risk = cp.quad_form(w, self.covmat)\n        prob = cp.Problem(cp.Minimize(risk),\n                          [cp.sum(w) == 1,\n                           w <= max_weight,\n                           w >= 0])\n\n        prob.solve()\n\n        return w.value\n\n    def maxret(self, max_weight=1.):\n        \"\"\"Compute maximum return portfolio.\"\"\"\n        w = cp.Variable(self.n_assets)\n        ret = self.mean @ w\n        prob = cp.Problem(cp.Maximize(ret),\n                          [cp.sum(w) == 1,\n                           w <= max_weight,\n                           w >= 0])\n\n        prob.solve()\n\n        return w.value\n\n    def maxutil(self, gamma, max_weight=1.):\n        \"\"\"Maximize utility.\"\"\"\n        w = cp.Variable(self.n_assets)\n        ret = self.mean @ w\n        risk = cp.quad_form(w, self.covmat)\n        prob = cp.Problem(cp.Maximize(ret - gamma * risk),\n                          [cp.sum(w) == 1,\n                           w <= max_weight,\n                           w >= 0])\n        prob.solve()\n\n        return w.value\n\n    def compute_portfolio_moments(self, weights):\n        \"\"\"Compute mean and standard deviation of some allocation.\n\n        Parameters\n        ----------\n        weights : np.array\n            1D array representing weights of a portfolio. Has shape `(n_assets,)`.\n\n        Returns\n        -------\n        mean : float\n            Mean (return)\n        std : float\n            Standard deviation (risk)\n        \"\"\"\n        pmean = np.inner(weights, self.mean)\n        pmean = pmean.item()\n        pvar = weights @ self.covmat @ weights\n        pvar = pvar.item()\n\n        return pmean, pvar ** (1 / 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now compute all the relevant optimal portfolios. We consider two cases with respect to\nthe :code:`max_weight` - constrained (:code:`max_weight<1`) and unconstrained\n(:code:`max_weight=1`).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "markowitz = Markowitz(mean.values, covmat.values)\nmax_weight = 2 / n_assets\ngamma = 10\n\noptimal_portfolios_u = {'minvar': markowitz.minvar(max_weight=1.),\n                        'maxret': markowitz.maxret(max_weight=1.),\n                        'maxutil': markowitz.maxutil(gamma=gamma, max_weight=1.),\n                        }\noptimal_portfolios_c = {'minvar': markowitz.minvar(max_weight=max_weight),\n                        'maxret': markowitz.maxret(max_weight=max_weight),\n                        'maxutil': markowitz.maxutil(gamma=gamma, max_weight=max_weight)}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's write some visualization machinery!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "color_mapping = {'minvar': 'r',\n                 'maxret': 'g',\n                 'maxutil': 'yellow'}\n\nmpl_config = {'s': 100,\n              'alpha': 0.65,\n              'linewidth': 0.5,\n              'edgecolor': 'black'}\n\nmarker_mapping = {'true_c': '*',\n                  'true_u': 'o',\n                  'emp_c': 'v',\n                  'emp_u': 'p',\n                  'deep_c': 'P',\n                  'deep_u': 'X'}\n\n\ndef plot_scatter(title='', **risk_ret_portfolios):\n    ax = df_risk_ret.plot.scatter(x='risk', y='return')\n    ax.set_xlim(x_lim)\n    ax.set_ylim(y_lim)\n    plt.title(title)\n    plt.tight_layout()\n\n    all_points = []\n    all_names = []\n\n    for name, portfolios in risk_ret_portfolios.items():\n        for objective, w in portfolios.items():\n            y, x = markowitz.compute_portfolio_moments(w)\n            all_points.append(ax.scatter(x,\n                                         y,\n                                         c=color_mapping[objective],\n                                         marker=marker_mapping[name],\n                                         **mpl_config))\n            all_names.append(\"{}_{}\".format(objective, name))\n\n    plt.legend(all_points,\n               all_names,\n               scatterpoints=1,\n               loc='lower left',\n               ncol=len(risk_ret_portfolios),\n               fontsize=8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we are ready to analyze the results. In the below figure we compare\nconstrained and unconstrained optimization.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_scatter(title='Ground truth: Unconstrained vs constrained',\n             true_u=optimal_portfolios_u,\n             true_c=optimal_portfolios_c)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Sample estimators of covmat and mean\nThe main question now is whether we are able to find the optimal allocation just from the\n:code:`returns`. One obvious way to do it is to explicitly estimate mean and covariance matrix\nand use the same optimizers as before to find the solution. If number of samples is big enough\nthis should work well!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "markowitz_emp = Markowitz(returns.mean().values, returns.cov().values)\n\nemp_portfolios_u = {'minvar': markowitz_emp.minvar(max_weight=1.),\n                    'maxret': markowitz_emp.maxret(max_weight=1.),\n                    'maxutil': markowitz_emp.maxutil(gamma=gamma, max_weight=1.)}\n\nemp_portfolios_c = {'minvar': markowitz_emp.minvar(max_weight=max_weight),\n                    'maxret': markowitz_emp.maxret(max_weight=max_weight),\n                    'maxutil': markowitz_emp.maxutil(gamma=gamma, max_weight=max_weight)}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let us see how the unconstrained portfolios compare with the ground truth.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_scatter(title='Ground truth vs empirical estimates: Unconstrained',\n             true_u=optimal_portfolios_u,\n             emp_u=emp_portfolios_u)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Similarly, let's look at the constrained version.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_scatter(title='Ground truth vs empirical estimates: Constrained',\n             true_c=optimal_portfolios_c,\n             emp_c=emp_portfolios_c)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can see that under this setup, the empirical estimates are yielding close to optimal\nportfolios.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Implementation and setup of a :code:`deepdow` network\nWe would like to construct a network that accepts no input features, however, has\nmultiple learnable parameters that determine the result of a forward pass. Namely,\nwe are going to use the `layers_numericalmarkowitz` layer. Our goal is to learn\na mean return vector, covariance matrix and gamma totally from scratch. Note that\nthis is fundamentally different from using sample estimators as done in the previous section.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class Net(torch.nn.Module, Benchmark):\n    \"\"\"Learn covariance matrix, mean vector and gamma.\n\n    One can enforce max weight per asset.\n    \"\"\"\n\n    def __init__(self, n_assets, max_weight=1.):\n        super().__init__()\n\n        self.force_symmetric = True\n        self.matrix = torch.nn.Parameter(torch.eye(n_assets), requires_grad=True)\n        self.exp_returns = torch.nn.Parameter(torch.zeros(n_assets), requires_grad=True)\n        self.gamma_sqrt = torch.nn.Parameter(torch.ones(1), requires_grad=True)\n\n        self.portfolio_opt_layer = NumericalMarkowitz(n_assets, max_weight=max_weight)\n\n    def forward(self, x):\n        \"\"\"Perform forward pass.\n\n        Parameters\n        ----------\n        x : torch.Tensor\n            Tensor of shape `(n_samples, n_channels, lookback, n_assets)`.\n\n        Returns\n        -------\n        weights_filled : torch.Tensor\n            Of shape (n_samples, n_assets) representing the optimal weights as determined by the\n            convex optimizer.\n\n        \"\"\"\n        n = len(x)\n\n        cov_sqrt = torch.mm(self.matrix, torch.t(self.matrix)) if self.force_symmetric else self.matrix\n\n        weights = self.portfolio_opt_layer(self.exp_returns[None, ...],\n                                           cov_sqrt[None, ...],\n                                           self.gamma_sqrt,\n                                           torch.zeros(1).to(device=x.device, dtype=x.dtype))\n        weights_filled = torch.repeat_interleave(weights, n, dim=0)\n\n        return weights_filled"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let us set some parameters. Note that :code:`lookback` is irrelevant since our network expects\nno features. We need to define it anyway to be able to use :code:`deepdow` dataloaders.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lookback, gap, horizon = 2, 0, 5  # lookback does not matter in this case, just used for consistency\nbatch_size = 256\nn_epochs = 100  # We employ early stopping"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We are ready to create a :code:`DeepDow` dataset.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "X_list, y_list = [], []\n\nfor i in range(lookback, n_timesteps - horizon - gap + 1):\n    X_list.append(returns.values[i - lookback: i, :])\n    y_list.append(returns.values[i + gap: i + gap + horizon, :])\n\nX = np.stack(X_list, axis=0)[:, None, ...]\ny = np.stack(y_list, axis=0)[:, None, ...]\n\ndataset = InRAMDataset(X, y, asset_names=returns.columns)\ndataloader = RigidDataLoader(dataset, batch_size=batch_size)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The main feature of :code:`deepdow` is that it only cares about the final allocation that minimizes\nsome function of **empirical** portfolio returns. Unlike with the sample estimators in the previous\nsections, one does not need to explicitly model the dynamics of the market and find allocation via\nthe two step procedure. Below we define empirical counterparts of the convex optimization objectives\n(losses). Note that all losses in :code:`deepdow` have the *the lower the better* logic.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "all_losses = {'minvar': StandardDeviation() ** 2,\n              'maxret': MeanReturns(),\n              'maxutil': MeanReturns() + gamma * StandardDeviation() ** 2}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Training and evaluation\nNow it is time to train!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "deep_portfolios_c = {}\ndeep_portfolios_u = {}\n\nfor mode in ['u', 'c']:\n    for loss_name, loss in all_losses.items():\n        network = Net(n_assets, max_weight=max_weight if mode == 'c' else 1.)\n        run = Run(network,\n                  loss,\n                  dataloader,\n                  val_dataloaders={'train': dataloader},\n                  callbacks=[EarlyStoppingCallback('train', 'loss', patience=3)])\n\n        run.launch(n_epochs=n_epochs)\n\n        # Results\n        w_pred = network(torch.ones(1, n_assets)).detach().numpy().squeeze()  # the input does not matter\n\n        if mode == 'c':\n            deep_portfolios_c[loss_name] = w_pred\n        else:\n            deep_portfolios_u[loss_name] = w_pred"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Unconstrained case:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_scatter(title='Ground truth vs empirical estimates vs deep learning: Unconstrained',\n             true_u=optimal_portfolios_u,\n             emp_u=emp_portfolios_u,\n             deep_u=deep_portfolios_u)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Constrained case:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_scatter(title='Ground truth vs empirical estimates vs deep learning: Constrained',\n             true_c=optimal_portfolios_c,\n             emp_c=emp_portfolios_c,\n             deep_c=deep_portfolios_c)"
      ]
    }
  ],
  "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.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}