{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d5ac471c",
   "metadata": {},
   "source": [
    "# Approximate 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",
    "## 1. Approximate Kurtosis Minimization\n",
    "\n",
    "The minimization of portfolio kurtosis can be approximated, as shown in __[Cajas (2023)](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4472793)__, as a convex optimization problem as follows:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "& \\min_{x} & &  \\left \\| g  \\right \\|_{2}  & \\\\\n",
    "& \\text{s.t.} & & \\left [ \\begin{array}{cc} X & x \\\\ x^{\\prime} & 1 \\end{array} \\right ] \\succeq 0 &\\\\\n",
    "& & & g_{i} = \\text{Tr}\\left ( B_{i} X \\right ) & \\forall i=1, \\ldots , K \\\\\n",
    "& & & \\sum^{n}_{i=1} x_{i} = 1 & \\\\\n",
    "& & & x \\geq 0 & \\\\\n",
    "& & & X = X^{\\prime} & \\\\\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where:\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 symmetric.\n",
    "<br>$K$: is the number of matrixes $B_{i}$ that are considered in the approximation of kurtosis, I recommend $K=2n$.\n",
    "<br>$B_{i} =  \\left ( \\text{vec} \\left(I_{n} \\right )^{\\prime} \\otimes I_m \\right)(I_n \\otimes \\sqrt{\\lambda_{i}} q_{i}) = \\text{vec}^{-1} \\left ( \\sqrt{\\lambda_{i}} q_{i} \\right )$.\n",
    "<br>$\\text{vec}$: is the vectorization operator.\n",
    "<br>$I_{n}$: is the identity matrix of size $n \\times n$.\n",
    "<br>$\\lambda_{i}$: is the eigenvalue $i$ of matrix $\\text{vec}^{n,n} \\left ( \\Sigma_{4} \\right )$.\n",
    "<br>$\\Sigma_{4}$: is the cokurtosis square matrix.\n",
    "<br>$\\text{vec}^{n,n}(.)$ is the rearrangement operator or block vectorization operator defined in __[Loan (1992)](https://ecommons.cornell.edu/handle/1813/5484)__ and __[Ojeda (2015)](https://arxiv.org/abs/1310.8149)__."
   ]
  },
  {
   "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.9692%</td>\n",
       "      <td>0.0180%</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.9758%</td>\n",
       "      <td>0.6987%</td>\n",
       "      <td>1.7539%</td>\n",
       "      <td>1.3734%</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.4783%</td>\n",
       "      <td>-2.2066%</td>\n",
       "      <td>-3.1557%</td>\n",
       "      <td>-1.6148%</td>\n",
       "      <td>-0.2700%</td>\n",
       "      <td>-0.5492%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016-01-08</th>\n",
       "      <td>-1.6037%</td>\n",
       "      <td>-2.5425%</td>\n",
       "      <td>0.1099%</td>\n",
       "      <td>-0.2241%</td>\n",
       "      <td>0.5707%</td>\n",
       "      <td>-1.7949%</td>\n",
       "      <td>-1.6410%</td>\n",
       "      <td>0.1720%</td>\n",
       "      <td>0.3067%</td>\n",
       "      <td>-0.1539%</td>\n",
       "      <td>-0.1448%</td>\n",
       "      <td>0.0895%</td>\n",
       "      <td>-3.3838%</td>\n",
       "      <td>-0.9720%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016-01-11</th>\n",
       "      <td>-1.6851%</td>\n",
       "      <td>-1.0215%</td>\n",
       "      <td>0.0914%</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.6436%</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.9692%  0.0180%  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.1099% -0.2241%  0.5707% -1.7949% -1.6410%   \n",
       "2016-01-11 -1.6851% -1.0215%  0.0914% -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.9758%  0.6987%  1.7539%  1.3734%  \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.4783% -2.2066% -3.1557% -1.6148% -0.2700% -0.5492%  \n",
       "2016-01-08  0.1720%  0.3067% -0.1539% -0.1448%  0.0895% -3.3838% -0.9720%  \n",
       "2016-01-11  2.0948% -0.0573%  1.6436% -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 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\n",
    "\n",
    "# Function that calculates block vectorization operator\n",
    "def block_vec_pq(A, p, q):\n",
    "    mp, nq = A.shape\n",
    "    if mp % p == 0 and nq % q == 0:\n",
    "        m = int(mp / p)\n",
    "        n = int(nq / q)\n",
    "        bvec_A = np.empty((0, p * q))\n",
    "        for j in range(n):\n",
    "            Aj = np.empty((0, p * q))\n",
    "            for i in range(m):\n",
    "                Aij = (\n",
    "                    A[i * p : (i + 1) * p, j * q : (j + 1) * q]\n",
    "                    .reshape(-1, 1, order=\"F\")\n",
    "                    .T\n",
    "                )\n",
    "                Aj = np.vstack([Aj, Aij])\n",
    "            bvec_A = np.vstack([bvec_A, Aj])\n",
    "    return bvec_A"
   ]
  },
  {
   "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>20.5587%</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>19.7034%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GE</th>\n",
       "      <td>0.0000%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GOOG</th>\n",
       "      <td>4.2250%</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>4.5958%</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>50.9171%</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       weights\n",
       "BAX    0.0000%\n",
       "BMY   20.5587%\n",
       "CMCSA  0.0000%\n",
       "CNP    0.0000%\n",
       "CPB   19.7034%\n",
       "GE     0.0000%\n",
       "GOOG   4.2250%\n",
       "MO     0.0000%\n",
       "MSFT   0.0000%\n",
       "NI     0.0000%\n",
       "SEE    4.5958%\n",
       "T      0.0000%\n",
       "TGT    0.0000%\n",
       "VZ    50.9171%"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "##################################################\n",
    "# Finding the Min Approximate 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",
    "K = 2*n\n",
    "\n",
    "# Calculate square cokurtosis matrix\n",
    "Sigma_4 = kurt_matrix(Y)\n",
    "\n",
    "# Calculate block vectorization of square cokurtosis matrix\n",
    "A = block_vec_pq(Sigma_4, n, n)\n",
    "# Finding eigenvalues and eigenvectors\n",
    "s_A, V_A = np.linalg.eig(A)\n",
    "s_A = np.clip(s_A, 0, np.inf)\n",
    "# Calculating matrixes Bi\n",
    "Bi = []\n",
    "for i in range(K):\n",
    "    B = s_A[i] ** 0.5 * V_A[:, i]\n",
    "    B = B.reshape((n, n), order=\"F\").real\n",
    "    Bi.append(B)\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((K, 1))\n",
    "\n",
    "risk = cp.pnorm(g,p=2)\n",
    "M = cp.bmat([[X, x], [x.T, np.ones((1, 1))]])\n",
    "constraints = [cp.sum(x) == 1,\n",
    "               x >= 0,\n",
    "               M >> 0]\n",
    "\n",
    "for i in range(K):\n",
    "    constraints += [g[i, 0] == cp.trace(Bi[i] @ X)]\n",
    "\n",
    "obj = cp.Minimize(risk * 1000)\n",
    "prob = cp.Problem(obj, constraints)\n",
    "\n",
    "#prob.solve(solver='MOSEK', mosek_params= {mosek.iparam.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.4726%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BMY</th>\n",
       "      <td>12.6221%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CMCSA</th>\n",
       "      <td>1.5312%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CNP</th>\n",
       "      <td>0.0001%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CPB</th>\n",
       "      <td>15.1031%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GE</th>\n",
       "      <td>0.0001%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GOOG</th>\n",
       "      <td>9.8516%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MO</th>\n",
       "      <td>7.1123%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSFT</th>\n",
       "      <td>0.0001%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NI</th>\n",
       "      <td>2.9518%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEE</th>\n",
       "      <td>4.3766%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T</th>\n",
       "      <td>0.0003%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TGT</th>\n",
       "      <td>6.4106%</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>VZ</th>\n",
       "      <td>31.5675%</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       weights\n",
       "BAX    8.4726%\n",
       "BMY   12.6221%\n",
       "CMCSA  1.5312%\n",
       "CNP    0.0001%\n",
       "CPB   15.1031%\n",
       "GE     0.0001%\n",
       "GOOG   9.8516%\n",
       "MO     7.1123%\n",
       "MSFT   0.0001%\n",
       "NI     2.9518%\n",
       "SEE    4.3766%\n",
       "T      0.0003%\n",
       "TGT    6.4106%\n",
       "VZ    31.5675%"
      ]
     },
     "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, 500180 characters]",
      "text/plain": [
       "<Figure size 1200x1000 with 4 Axes>"
      ]
     },
     "metadata": {},
     "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 Approximate 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 Approximate 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, 159176 characters]",
      "text/plain": [
       "<Figure size 1000x1000 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 Approximate 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, 257752 characters]",
      "text/plain": [
       "<Figure size 1000x800 with 2 Axes>"
      ]
     },
     "metadata": {},
     "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 Approximate Kurtosis Portfolio Histogram')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "156c9909",
   "metadata": {},
   "source": [
    "We can see that the minimum approximate 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 is available in the CVXPY based library __[Riskfolio-Lib](https://github.com/dcajasn/Riskfolio-Lib)__.\n",
    "\n",
    "Finally, I hope you liked this example."
   ]
  }
 ],
 "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
