{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d5ac471c",
   "metadata": {},
   "source": [
    "# Kurtosis Portfolio Optimization\n",
    "\n",
    "In this notebook we show how to use the semidefinite cone and second order cone to model the optimization of portfolio kurtosis.\n",
    "\n",
    "## <font color='red'> Warning: Don't use this model for n >= 30 without Mosek and I don't recommend to calculate $L_2$, $D_2$, $S_2$ and $\\Sigma_4$ for n >= 30 with Python (Use C++ and Armadillo library) </font>\n",
    "\n",
    "## 1. Kurtosis Minimization\n",
    "\n",
    "The minimization of portfolio kurtosis, as shown in __[Cajas (2022)](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4202967)__, can be posed as a convex problem as follows:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\min_{x, \\, X, \\, z, \\, g} \\; &  g \\\\\n",
    "& \\left [ \\begin{array}{cc} X & x \\\\ x^{\\prime} & 1 \\end{array} \\right ] \\succeq 0 \\\\\n",
    "& z = L_{2} \\text{vec}(X) \\\\\n",
    "& \\left \\| \\left [ S_{2} \\Sigma_{4} S^{\\prime}_{2} \\right ]^{1/2} z \\right \\|_{2} \\leq g \\\\\n",
    "& \\sum^{n}_{i=1}x_{i} = 1 \\\\\n",
    "& x \\geq 0 \\\\\n",
    "& X \\succeq 0 \\\\\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where:\n",
    "<br>$g$: is an auxiliary variable that represents the square root of portfolio kurtosis\n",
    "<br>$\\Sigma_{4}$: is the square co-kurtosis matrix.\n",
    "<br>$x$: is the vector of portfolio weights of size $n \\times 1$\n",
    "<br>$X$: is an auxiliary variable of size $n \\times n$ that is positive semidefinite.\n",
    "<br>$D_{2}$: is the duplication efined in __[Magnus (1980)](https://pure.uvt.nl/ws/portalfiles/portal/649691/26951_6623.pdf)__.\n",
    "<br>$L_{2}$: is the duplication elimination matrix defined in __[Magnus (1980)](https://pure.uvt.nl/ws/portalfiles/portal/649691/26951_6623.pdf)__.\n",
    "<br>$S_{2}$ is the duplication summation matrix defined as $S_{2} = D^{\\prime}_{2}D_{2}L_{2}$\n",
    "<br>$\\text{vec}$: is the vectorization operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ec087b00",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[*********************100%***********************]  14 of 14 completed\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BAX</th>\n",
       "      <th>BMY</th>\n",
       "      <th>CMCSA</th>\n",
       "      <th>CNP</th>\n",
       "      <th>CPB</th>\n",
       "      <th>GE</th>\n",
       "      <th>GOOG</th>\n",
       "      <th>MO</th>\n",
       "      <th>MSFT</th>\n",
       "      <th>NI</th>\n",
       "      <th>SEE</th>\n",
       "      <th>T</th>\n",
       "      <th>TGT</th>\n",
       "      <th>VZ</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2016-01-05</th>\n",
       "      <td>0.4036%</td>\n",
       "      <td>1.9693%</td>\n",
       "      <td>0.0179%</td>\n",
       "      <td>0.9305%</td>\n",
       "      <td>0.3678%</td>\n",
       "      <td>0.0977%</td>\n",
       "      <td>0.0998%</td>\n",
       "      <td>2.0213%</td>\n",
       "      <td>0.4562%</td>\n",
       "      <td>1.5881%</td>\n",
       "      <td>0.9759%</td>\n",
       "      <td>0.6987%</td>\n",
       "      <td>1.7539%</td>\n",
       "      <td>1.3735%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016-01-06</th>\n",
       "      <td>0.2412%</td>\n",
       "      <td>-1.7557%</td>\n",
       "      <td>-0.7727%</td>\n",
       "      <td>-1.2473%</td>\n",
       "      <td>-0.1736%</td>\n",
       "      <td>-1.5940%</td>\n",
       "      <td>0.1400%</td>\n",
       "      <td>1.0589%</td>\n",
       "      <td>-1.8165%</td>\n",
       "      <td>0.5547%</td>\n",
       "      <td>-1.5647%</td>\n",
       "      <td>0.3108%</td>\n",
       "      <td>-1.0155%</td>\n",
       "      <td>-0.9034%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016-01-07</th>\n",
       "      <td>-1.6573%</td>\n",
       "      <td>-2.7699%</td>\n",
       "      <td>-1.1047%</td>\n",
       "      <td>-1.9769%</td>\n",
       "      <td>-1.2206%</td>\n",
       "      <td>-4.2314%</td>\n",
       "      <td>-2.3170%</td>\n",
       "      <td>-1.7407%</td>\n",
       "      <td>-3.4782%</td>\n",
       "      <td>-2.2066%</td>\n",
       "      <td>-3.1557%</td>\n",
       "      <td>-1.6148%</td>\n",
       "      <td>-0.2700%</td>\n",
       "      <td>-0.5493%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016-01-08</th>\n",
       "      <td>-1.6037%</td>\n",
       "      <td>-2.5425%</td>\n",
       "      <td>0.1098%</td>\n",
       "      <td>-0.2241%</td>\n",
       "      <td>0.5707%</td>\n",
       "      <td>-1.7950%</td>\n",
       "      <td>-1.6410%</td>\n",
       "      <td>0.1720%</td>\n",
       "      <td>0.3067%</td>\n",
       "      <td>-0.1538%</td>\n",
       "      <td>-0.1448%</td>\n",
       "      <td>0.0895%</td>\n",
       "      <td>-3.3839%</td>\n",
       "      <td>-0.9719%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016-01-11</th>\n",
       "      <td>-1.6851%</td>\n",
       "      <td>-1.0215%</td>\n",
       "      <td>0.0915%</td>\n",
       "      <td>-1.1791%</td>\n",
       "      <td>0.5674%</td>\n",
       "      <td>0.4569%</td>\n",
       "      <td>0.2183%</td>\n",
       "      <td>2.0948%</td>\n",
       "      <td>-0.0573%</td>\n",
       "      <td>1.6435%</td>\n",
       "      <td>-0.1450%</td>\n",
       "      <td>1.2224%</td>\n",
       "      <td>1.4570%</td>\n",
       "      <td>0.5800%</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                BAX      BMY    CMCSA      CNP      CPB       GE     GOOG  \\\n",
       "Date                                                                        \n",
       "2016-01-05  0.4036%  1.9693%  0.0179%  0.9305%  0.3678%  0.0977%  0.0998%   \n",
       "2016-01-06  0.2412% -1.7557% -0.7727% -1.2473% -0.1736% -1.5940%  0.1400%   \n",
       "2016-01-07 -1.6573% -2.7699% -1.1047% -1.9769% -1.2206% -4.2314% -2.3170%   \n",
       "2016-01-08 -1.6037% -2.5425%  0.1098% -0.2241%  0.5707% -1.7950% -1.6410%   \n",
       "2016-01-11 -1.6851% -1.0215%  0.0915% -1.1791%  0.5674%  0.4569%  0.2183%   \n",
       "\n",
       "                 MO     MSFT       NI      SEE        T      TGT       VZ  \n",
       "Date                                                                       \n",
       "2016-01-05  2.0213%  0.4562%  1.5881%  0.9759%  0.6987%  1.7539%  1.3735%  \n",
       "2016-01-06  1.0589% -1.8165%  0.5547% -1.5647%  0.3108% -1.0155% -0.9034%  \n",
       "2016-01-07 -1.7407% -3.4782% -2.2066% -3.1557% -1.6148% -0.2700% -0.5493%  \n",
       "2016-01-08  0.1720%  0.3067% -0.1538% -0.1448%  0.0895% -3.3839% -0.9719%  \n",
       "2016-01-11  2.0948% -0.0573%  1.6435% -0.1450%  1.2224%  1.4570%  0.5800%  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "####################################\n",
    "# Downloading Data\n",
    "####################################\n",
    "!pip install --quiet yfinance\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import yfinance as yf\n",
    "from itertools import permutations\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mtick\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "yf.pdr_override()\n",
    "pd.options.display.float_format = '{:.4%}'.format\n",
    "\n",
    "# Date range\n",
    "start = '2016-01-01'\n",
    "end = '2021-12-30'\n",
    "\n",
    "# Tickers of assets\n",
    "assets = ['TGT', 'CMCSA', 'CPB', 'MO', 'T', 'BAX', 'BMY',\n",
    "          'MSFT', 'SEE', 'VZ', 'CNP', 'NI', 'GE', 'GOOG']\n",
    "assets.sort()\n",
    "\n",
    "# Downloading data\n",
    "data = yf.download(assets, start = start, end = end)\n",
    "data = data.loc[:,('Adj Close', slice(None))]\n",
    "data.columns = assets\n",
    "\n",
    "# Calculating returns\n",
    "Y = data[assets].pct_change().dropna()\n",
    "\n",
    "display(Y.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c25573b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "####################################\n",
    "# Auxiliary functions\n",
    "####################################\n",
    "\n",
    "# Function that calculates D_2\n",
    "def duplication_matrix(n):\n",
    "    out = np.zeros((int(n*(n+1)/2), n**2))\n",
    "    for j in range(1, n+1):\n",
    "        for i in range(j, n+1):\n",
    "            u = np.zeros((int(n*(n+1)/2), 1))\n",
    "            u[round((j-1)*n + i - ((j-1)*j)/2)-1] = 1.0\n",
    "            E = np.zeros((n,n))\n",
    "            E[i-1,j-1] = 1.0\n",
    "            E[j-1,i-1] = 1.0\n",
    "            out += u @ E.reshape(-1,1).T;\n",
    "    return out.T\n",
    "\n",
    "# Function that calculates L_2\n",
    "def duplication_elimination_matrix(n):\n",
    "    out = np.zeros((int(n*(n+1)/2), n**2))\n",
    "    for j in range(n):\n",
    "        e_j = np.zeros((1, n))\n",
    "        e_j[0, j] = 1.0\n",
    "        for i in range(j, n):\n",
    "            u = np.zeros((int(n*(n+1)/2), 1))      \n",
    "            row = round(j*n+i-((j+1)*j)/2)\n",
    "            u[row] = 1.0\n",
    "            e_i = np.zeros((1, n))\n",
    "            e_i[0,i] = 1.0\n",
    "            out += np.kron(u, np.kron(e_j, e_i))\n",
    "    return out\n",
    "\n",
    "# Function that calculates S_4\n",
    "def kurt_matrix(Y):\n",
    "    P = Y.to_numpy()\n",
    "    T, n = P.shape\n",
    "    mu = np.mean(P, axis=0).reshape(1,-1)\n",
    "    mu = np.repeat(mu, T, axis=0)\n",
    "    x = P - mu\n",
    "    ones = np.ones((1,n))\n",
    "    z = np.kron(ones, x) * np.kron(x, ones);\n",
    "    S4 = 1/T * z.T @ z\n",
    "    return S4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0f4caac4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>weights</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>BAX</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BMY</th>\n",
       "      <td>21.0918%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CMCSA</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CNP</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CPB</th>\n",
       "      <td>22.6310%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GE</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GOOG</th>\n",
       "      <td>8.5225%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MO</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSFT</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NI</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEE</th>\n",
       "      <td>5.0681%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TGT</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>VZ</th>\n",
       "      <td>42.6866%</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       weights\n",
       "BAX    0.0000%\n",
       "BMY   21.0918%\n",
       "CMCSA  0.0000%\n",
       "CNP    0.0000%\n",
       "CPB   22.6310%\n",
       "GE     0.0000%\n",
       "GOOG   8.5225%\n",
       "MO     0.0000%\n",
       "MSFT   0.0000%\n",
       "NI     0.0000%\n",
       "SEE    5.0681%\n",
       "T      0.0000%\n",
       "TGT    0.0000%\n",
       "VZ    42.6866%"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "####################################\n",
    "# Finding the Min Kurtosis Portfolio\n",
    "####################################\n",
    "\n",
    "import cvxpy as cp\n",
    "#import mosek\n",
    "from scipy.linalg import sqrtm\n",
    "\n",
    "# Parameters \n",
    "T, n = Y.shape\n",
    "\n",
    "D_2 = duplication_matrix(n)\n",
    "L_2 = duplication_elimination_matrix(n)\n",
    "S_2 = D_2.T @ D_2 @ L_2\n",
    "\n",
    "Sigma_4 = kurt_matrix(Y)\n",
    "Sigma_4_sqrt = sqrtm(S_2 @ Sigma_4 @ S_2.T)\n",
    "\n",
    "# Variables\n",
    "rows = int(n*(n+1)/2)\n",
    "x = cp.Variable((n,1))\n",
    "z = cp.Variable((rows,1))\n",
    "X = cp.Variable((n,n), PSD=True)\n",
    "g = cp.Variable()\n",
    "\n",
    "risk = g\n",
    "constraints = [cp.sum(x) == 1,\n",
    "               x >= 0,\n",
    "               z == L_2 @ cp.reshape(cp.vec(X), (n*n,1)),\n",
    "              ]\n",
    "\n",
    "M3 = cp.bmat([[X, x], [x.T, np.ones((1, 1))]])\n",
    "\n",
    "constraints += [\n",
    "                cp.SOC(g, Sigma_4_sqrt @ z),\n",
    "                M3 >> 0,\n",
    "               ]\n",
    "\n",
    "obj = cp.Minimize(risk * 1000)\n",
    "prob = cp.Problem(obj, constraints)\n",
    "\n",
    "#prob.solve(solver='MOSEK', mosek_params= {'MSK_IPAR_NUM_THREADS': 2})\n",
    "prob.solve()\n",
    "w = pd.DataFrame(x.value, index=assets, columns=['weights'])\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8abf1eb1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>weights</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>BAX</th>\n",
       "      <td>8.4725%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BMY</th>\n",
       "      <td>12.6225%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CMCSA</th>\n",
       "      <td>1.5316%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CNP</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CPB</th>\n",
       "      <td>15.1033%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GE</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GOOG</th>\n",
       "      <td>9.8518%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MO</th>\n",
       "      <td>7.1117%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSFT</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NI</th>\n",
       "      <td>2.9515%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEE</th>\n",
       "      <td>4.3759%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T</th>\n",
       "      <td>0.0001%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TGT</th>\n",
       "      <td>6.4104%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>VZ</th>\n",
       "      <td>31.5685%</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       weights\n",
       "BAX    8.4725%\n",
       "BMY   12.6225%\n",
       "CMCSA  1.5316%\n",
       "CNP    0.0000%\n",
       "CPB   15.1033%\n",
       "GE     0.0000%\n",
       "GOOG   9.8518%\n",
       "MO     7.1117%\n",
       "MSFT   0.0000%\n",
       "NI     2.9515%\n",
       "SEE    4.3759%\n",
       "T      0.0001%\n",
       "TGT    6.4104%\n",
       "VZ    31.5685%"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "####################################\n",
    "# Finding the Min Variance Portfolio\n",
    "####################################\n",
    "\n",
    "!pip install --quiet riskfolio-lib\n",
    "import riskfolio as rp\n",
    "\n",
    "# Building the portfolio object\n",
    "port = rp.Portfolio(returns=Y)\n",
    "\n",
    "# Calculating optimum portfolio\n",
    "\n",
    "# Select method and estimate input parameters:\n",
    "\n",
    "method_mu='hist' # Method to estimate expected returns based on historical data.\n",
    "method_cov='hist' # Method to estimate covariance matrix based on historical data.\n",
    "\n",
    "port.assets_stats(method_mu=method_mu, method_cov=method_cov, d=0.94)\n",
    "\n",
    "# Estimate optimal portfolio:\n",
    "\n",
    "port.solvers = [cp.MOSEK]\n",
    "model='Classic' # Could be Classic (historical), BL (Black Litterman) or FM (Factor Model)\n",
    "rm = 'MV' # Risk measure used, this time will be EVaR\n",
    "obj = 'MinRisk' # Objective function, could be MinRisk, MaxRet, Utility or Sharpe\n",
    "hist = True # Use historical scenarios for risk measures that depend on scenarios\n",
    "rf = 0 # Risk free rate\n",
    "l = 0 # Risk aversion factor, only useful when obj is 'Utility'\n",
    "\n",
    "w2 = port.optimization(model=model, rm=rm, obj=obj, rf=rf, l=l, hist=hist)\n",
    "\n",
    "display(w2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "01cea5bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "####################################\n",
    "# Montecarlo Simulation\n",
    "####################################\n",
    "\n",
    "# Montecarlo simulation of portfolio weights\n",
    "rs = np.random.RandomState(seed=123)\n",
    "s1 = rs.dirichlet([0.1] * len(assets), 10000)\n",
    "s2 = rs.dirichlet([0.25] * len(assets), 10000)\n",
    "s3 = rs.dirichlet([0.5] * len(assets), 10000)\n",
    "s4 = rs.dirichlet([0.75] * len(assets), 10000)\n",
    "s5 = rs.dirichlet([1.0] * len(assets), 10000)\n",
    "s6 = rs.dirichlet([1.5] * len(assets), 10000)\n",
    "s7 = rs.dirichlet([2.0] * len(assets), 10000)\n",
    "s8 = rs.dirichlet([3.0] * len(assets), 10000)\n",
    "sample = np.concatenate([np.identity(len(assets)), s1, s2, s3, s4, s5, s6, s7, s8], axis=0)\n",
    "\n",
    "# Calculating mean, standard deviation and square root kurtosis of each portfolio\n",
    "m = sample.shape[0]\n",
    "M_1 = np.mean(Y.to_numpy(), axis=0).reshape(1,-1)\n",
    "M_2 = Y.cov().to_numpy()\n",
    "\n",
    "c_mean = M_1 @ sample.T\n",
    "c_var = np.zeros(m)\n",
    "c_kurt = np.zeros(m)\n",
    "\n",
    "for i in range(0, m):\n",
    "    c_var[i] = (sample[i] @ M_2 @ sample[i].T) ** (0.5)\n",
    "    c_kurt[i] = (np.kron(sample[i], sample[i]) @ Sigma_4 @ np.kron(sample[i],sample[i]).T) ** (1/4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a108ffdb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 301780 characters]\n",
      "text/plain": [
       "<Figure size 864x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "####################################\n",
    "# Plotting Portfolios\n",
    "####################################\n",
    "\n",
    "fig, ax = plt.subplots(2, 1, figsize=(12,10))\n",
    "ax = np.ravel(ax)\n",
    "\n",
    "# Plotting Portfolios in mean-standard deviation plane\n",
    "cax0 = ax[0].scatter(c_var, c_mean, c=c_mean/c_var, cmap='Spectral')\n",
    "ax[0].scatter((w.to_numpy().T @ M_2 @ w.to_numpy()) ** (0.5),\n",
    "              M_1 @ w.to_numpy(),\n",
    "              marker='*',\n",
    "              s=2**8,\n",
    "              color='fuchsia',\n",
    "              label='Minimum Kurtosis')\n",
    "ax[0].scatter((w2.to_numpy().T @ M_2 @ w2.to_numpy()) ** (0.5),\n",
    "              M_1 @ w2.to_numpy(),\n",
    "              marker='*',\n",
    "              s=2**8,\n",
    "              color='tab:red',\n",
    "              label='Minimum Variance')\n",
    "\n",
    "# Plotting Portfolios in mean-square root kurtosis plane\n",
    "cax1 = ax[1].scatter(c_kurt, c_mean, c=c_mean/c_kurt**2, cmap='Spectral')\n",
    "ax[1].scatter((np.kron(w.to_numpy(), w.to_numpy()).T @ Sigma_4 @ np.kron(w.to_numpy(), w.to_numpy())) ** (0.25),\n",
    "              M_1 @ w.to_numpy(),\n",
    "              marker='*',\n",
    "              s=2**8,\n",
    "              color='fuchsia',\n",
    "              label='Minimum Kurtosis')\n",
    "ax[1].scatter((np.kron(w2.to_numpy(), w2.to_numpy()).T @ Sigma_4 @ np.kron(w2.to_numpy(), w2.to_numpy())) ** (0.25),\n",
    "              M_1 @ w2.to_numpy(),\n",
    "              marker='*',\n",
    "              s=2**8,\n",
    "              color='tab:red',\n",
    "              label='Minimum Variance')\n",
    "\n",
    "fig.colorbar(cax0, ax=ax[0], label='Expected Return / Standard Deviation')\n",
    "fig.colorbar(cax1, ax=ax[1], label='Expected Return / Square Root Kurtosis')\n",
    "\n",
    "ax[0].xaxis.set_major_formatter(mtick.PercentFormatter(xmax=1))\n",
    "ax[0].yaxis.set_major_formatter(mtick.PercentFormatter(xmax=1))\n",
    "ax[0].set_xlabel('Daily Standard Deviation')\n",
    "ax[0].set_ylabel('Daily Espected Return')\n",
    "ax[0].legend()\n",
    "\n",
    "ax[1].xaxis.set_major_formatter(mtick.PercentFormatter(xmax=1))\n",
    "ax[1].yaxis.set_major_formatter(mtick.PercentFormatter(xmax=1))\n",
    "ax[1].set_xlabel('Daily Fourth Root Kurtosis')\n",
    "ax[1].set_ylabel('Daily Espected Return')\n",
    "ax[1].legend()\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f48b5e9e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 105208 characters]\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "####################################\n",
    "# Plotting Portfolios Composition\n",
    "####################################\n",
    "\n",
    "fig, ax = plt.subplots(2, 1, figsize=(10,10))\n",
    "ax = np.ravel(ax)\n",
    "\n",
    "rp.plot_pie(w=w2,\n",
    "            title='Minimum Variance Portfolio',\n",
    "            others=0.05,\n",
    "            nrow=25,\n",
    "            ax=ax[0])\n",
    "\n",
    "rp.plot_pie(w=w,\n",
    "            title='Minimum Kurtosis Portfolio',\n",
    "            others=0.05,\n",
    "            nrow=25,\n",
    "            ax=ax[1])\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "31e4892f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 133424 characters]\n",
      "text/plain": [
       "<Figure size 720x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "####################################\n",
    "# Plotting Portfolios Histogram\n",
    "####################################\n",
    "\n",
    "fig, ax = plt.subplots(2, 1, figsize=(10,8))\n",
    "ax = np.ravel(ax)\n",
    "\n",
    "rp.plot_hist(returns=Y,\n",
    "             w=w2,\n",
    "             alpha=0.05,\n",
    "             bins=50,\n",
    "             ax=ax[0])\n",
    "\n",
    "rp.plot_hist(returns=Y,\n",
    "             w=w,\n",
    "             alpha=0.05,\n",
    "             bins=50,\n",
    "             ax=ax[1])\n",
    "\n",
    "ax[0].set_title('Minimum Variance Portfolio Histogram')\n",
    "ax[1].set_title('Minimum Kurtosis Portfolio Histogram')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "156c9909",
   "metadata": {},
   "source": [
    "We can see that the minimum kurtosis portfolio reduce the tails of portfolio returns distribution compared to minimun variance portfolio because more conservative downside risk measures like Entropic Value at Risk and Worst realization are lower for the minimum kurtosis portfolio.\n",
    "\n",
    "This portfolio optimization model will be added soon to the CVXPY based library __[Riskfolio-Lib](https://github.com/dcajasn/Riskfolio-Lib)__."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
