ipynb/nikolai_chow-qe_notebook_-_optimization_with_stochastic_dominance_constraints.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "id": "e417a7c4",
            "metadata": {
                "slideshow": {
                    "slide_type": "slide"
                }
            },
            "source": [
                "# Using linear programming to solve a portfolio optimization problem with stochastic dominance constraints"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "0f64c362",
            "metadata": {
                "slideshow": {
                    "slide_type": "slide"
                }
            },
            "source": [
                "## Overview\n",
                "\n",
                "In this notebook, we will briefly discuss the following:\n",
                "\n",
                "1. The concept of stochastic dominance\n",
                "\n",
                "2. A portfolio optimization problem with stochastic dominance constraints that can be solved by linear programing\n",
                "\n",
                "3. An example in \"Dentcheva, D., & Ruszczyński, A. (2003). Optimization with stochastic dominance constraints. SIAM Journal on Optimization.\" with Python"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "1fc14da2",
            "metadata": {},
            "source": [
                "## 1. The concept of stochastic dominance\n",
                "\n",
                "Let’s start with some notations and definitions. \n",
                "\n",
                "Let random variable $X\\in \\mathcal{L}^1(\\Omega,\\mathcal{F},P)$ such that $X: \\Omega \\rightarrow \\mathbb{R}$, where $(\\Omega,\\mathcal{F},P)$ is the probability space.\n",
                "\n",
                "The cumulative distribution function (CDF) of $X$ is defined as:\n",
                "\n",
                "$$F(\\eta) = P[X \\leq \\eta]$$\n",
                "\n",
                "The second CDF of $X$ is defined as:\n",
                "\n",
                "$$F_2(\\eta) = \\int^\\eta_{-\\infty}F(t)dt$$ for $\\eta \\in R$.\n",
                "\n",
                "\n",
                "Similarly, for random variable $Y\\in \\mathcal{L}^1(\\Omega,\\mathcal{F},P)$ such that $Y: \\Omega \\rightarrow \\mathbb{R}$, we have $G(\\eta)$ and $G_2(\\eta)$.\n",
                "\n",
                "With the above setting, stochastic dominance relationships are defined as:\n",
                "\n",
                "1. $G$ is said to be first order stochastic dominated by $F$  ($F\\succeq_1 G$) if\n",
                "\n",
                "$$F(\\eta)\\leq G(\\eta) \\;\\text{for all}\\; \\eta \\in  \\mathbb{R}.$$\n",
                "\n",
                "\n",
                "2. $G$ is said to be second order stochastic dominated by $F$  ($F\\succeq_2 G$) if\n",
                "\n",
                "$$F_2(\\eta)\\leq G_2(\\eta) \\;\\text{for all}\\; \\eta \\in  \\mathbb{R}.$$"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "00c00efc",
            "metadata": {},
            "source": [
                "Stochastic dominance is a tool for comparing two random variables or two probability distributions. It allows us to determine which distribution is \"better\" in a certain sense.\n",
                "\n",
                "First order stochastic dominance (FSD), or $F\\succeq_1 G$, indicates that for any given reference outcome $\\eta$, the probability of getting outcomes lower than $\\eta$ is always weakly smaller for $F$ compared to $G$ for all $\\eta \\in  \\mathbb{R}$. This means $F$ is less likely to produce outcomes that are worse than $\\eta$ compared to $G$. Second order stochastic dominance (SSD), or $F\\succeq_2 G$, is similar to first order stochastic dominance, but different in the sense that it concerns about the integral of cumulative probabilities. \n",
                "\n",
                "In the graphs shown below, we compare two CDFs. The CDFs are denoted by $F$ and $G$. In the left graph, $F$ is normal CDF with a mean ($\\mu_F$) of 0.2 and a standard deviation ($\\sigma_F$) of 0.1 while $G$ is normal CDF with a mean ($\\mu_G$) of 0.1 and a standard deviation ($\\sigma_G$) of 0.1. $F$ is represented by the red line, while $G$ is represented by the blue line. In the right graph, $F$ is still a normal CDF with a mean of 0.2 and a standard deviation of 0.1, but $G$ has a mean of 0.1 and a standard deviation of 0.2. As before, $F$ is represented by the red line and $G$ is represented by the blue line."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "105d18ff",
            "metadata": {},
            "source": [
                "![Picture5](https://lh6.googleusercontent.com/bxxuLvf3yl4bf5YjTFDz8R7wNmDChoA6mzhfTVjI9PvROm2qKX_HO604FV2lVgo3QMw=w2400)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "f07ffe56",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "source": [
                "### 1.1. Why should we care?\n",
                "\n",
                "There are several reasons why stochastic dominance is a valuable tool in economics and finance. Here are two key theorems from Hadar and Russell (1969), Hanoch and Levy (1975), and Tesfatsion (1976) that highlight the importance of stochastic dominance:\n",
                "\n",
                "1. The relation $F\\succeq_1 Y$ is equivalent to $E[u(X)]\\geq E[u(Y)]$\n",
                "for all utility function $u$ with $u'\\geq 0$\n",
                "\n",
                "\n",
                "2. The relation $F\\succeq_2 G$ is equivalent to $E[u(X)]\\geq E[u(Y)]$\n",
                "for all utility function $u$ with $u'\\geq 0$ and $u''\\leq 0$\n",
                "\n",
                "In short, stochastic dominance is a useful tool for analyzing decision-making under risk because it has two key advantages:\n",
                "\n",
                "1. It allows us to compare the expected utility of different choices by comparing the CDFs of the corresponding random variables. This allows us to make conclusions about the relative utility of different options without specifying a particular utility function.\n",
                "\n",
                "2. The conclusions we draw from comparing CDFs are valid for a class of utility functions, rather than just a single specific function. \n",
                "\n",
                "If you want to learn more about the topic of stochastic dominance, you can refer to the following sources:\n",
                "\n",
                "- Levy (2016) provides a thorough discussion of the concept and its applications.\n",
                "- Rothschild and Stiglitz (1978) is a classic paper that is useful if you want to discuss the topic with decision theorists.\n",
                "- Atkinson (1970), Cho et al. (2007), Maasoumi et al. (2005), Guerre et al. (2009) and Chui et al. (2020) are all examples of articles that explore different applications of stochastic dominance in various contexts."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "1e1fbbef",
            "metadata": {},
            "source": [
                "## 2. A portfolio optimization problem with stochastic dominance constraints\n",
                "\n",
                "One direction of stochastic dominance research is to find out the way to form SSD-efficient portfolios. In general, this is not an easy task. The reason is that $\\succeq_2$ is a partial order, meaning that it is often not possible to use a real value function to rank the choices under consideration. As a result, it is not possible to use traditional optimization techniques to find portfolios with the highest \"SSD rank\". In this section, we will examine one specific portfolio optimization problem that incorporates the concept of SSD."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "2d626367",
            "metadata": {
                "slideshow": {
                    "slide_type": "slide"
                }
            },
            "source": [
                "### 2.1. A portfolio optimization problem\n",
                "\n",
                "Let $R_1,..., R_N $ be random returns for assets $1,...,N$,  with joint cumulative distribution function $\\mathcal{F}(\\textbf{r})$ where $\\textbf{r}\\in R_1 \\times ... \\times R_N$. Assume that $E[|R_n|]< \\infty$ for all $n=1,...,N$, and further let $w_1,...,w_N$ as portfolio weights on assets $1,...,N$. Without short selling, the set of possible asset allocations is defined as:\n",
                "\n",
                "$$W=\\Big\\{w\\in \\mathbb{R}^N: \\sum^N_{i=1} w_i=1, w_n\\geq 0 \\;\\; \\text{for}\\;\\; n=1,...,N\\Big\\}.$$\n",
                "\n",
                "The CDF $\\mathcal{F}_w(\\eta)$ for asset allocation $w\\in W$ is defined as:\n",
                "\n",
                "$$\\mathcal{F}_w(\\eta)=\\int_{\\{\\textbf{r}^Tw\\leq \\eta\\}}d\\mathcal{F}(\\textbf{r})$$\n",
                "\n",
                "where the corresponding portfolio return is denoted as $R_w=R_1w_1+...+R_Nw_N$. For empirical purposes, we assume that by selecting appropriate values for $a$ and $b$, we can ensure that the $\\mathcal{F}_w(a) = 0$ and $\\mathcal{F}_w(b) = 1$ for all $w \\in W$."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "b69315ef",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "source": [
                "Letting $Y$ as target random return with finite expected value and corresponding CDF $G$, such that $G(a)=0$ and $G(b)=1$. Under these settings, one could consider the following problem: \n",
                "\n",
                "\\begin{align}\n",
                "    \\max E[R_w] \\nonumber\\\\\n",
                "    \\text{subject to} \\;\\; \\mathcal{F}_w\\succeq_{2} G\\nonumber\\\\\n",
                "    w \\in W \n",
                "\\end{align}\n",
                "where $E[R_w]$ is the expected portoflio return.\n",
                "\n",
                "The above problem is first considered in Dentcheva and Ruszczyński (2003). In plain english, this problem is about finding a portfolio with highest expected return given that all risk-averse investors (or investors with concave utility function) will prefer the portfolio over the target random return (e.g. the return of S&P 500). Solving this problem could tell us whether the given porfolio $Y$ is SSD-efficent and how to get a \"better\" portfolio if it exists. To understand issues related to optimality, it may be helpful to consult the work of Dentcheva and Ruszczyński (2003). Intuitively, solution must exist if $Y = R_w$ for some $w \\in W$ since the target random return is no worse than itself in SSD sense. \n",
                "\n",
                "Solving the above problem can be challenging since it involves computing intergals of CDF $\\mathcal{F}_w(\\eta)$ and $\\mathcal{F}_w(\\eta)$ will vary significantly when the weights on assets are different. To address this, Dentcheva and Ruszczyński (2003) use a three steps approach:\n",
                "\n",
                "1. Transform the problem into an equivalent problem that doesn't require computing intergals of marginal CDF.\n",
                "2. Show that if $Y$ has a discrete distribution with known realizations, we don't need to compare intergals of CDF for all $\\eta \\in [a,b]$.\n",
                "3. Further transform the problem into a linear programming problem."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "dafbb17c",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "source": [
                "### 2.2. An equivalent problem\n",
                "\n",
                "Bawa et al.(1985) states the following:\n",
                "\n",
                "$$F_2(\\eta)\\leq G_2(\\eta) \\;\\text{for all}\\; \\eta \\in [a,b]$$\n",
                "\n",
                "$$\\Longleftrightarrow$$\n",
                "\n",
                "\\begin{aligned}\n",
                "E[(\\eta - X)_+]\\leq E[(\\eta -Y)_+] \\ & \\ \\ \\text{for all} \\ \\eta\\in [a,b]\n",
                "\\end{aligned}\n",
                "\n",
                "where $(\\eta - X)_+=\\max\\{\\eta - X,0\\}$. $E[(\\eta - X)_+]$ is the lower partial moment of X and $E[(\\eta -Y)_+]$ is the lower partial moment of Y. Using this idea, one could transform the original problem into an equivalent problem as: "
            ]
        },
        {
            "cell_type": "markdown",
            "id": "1be4ade6",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "source": [
                "\\begin{align}\n",
                "    \\max &\\; \\;\\;  E[R_w] \\nonumber\\\\\n",
                "    \\text{subject to} \n",
                "    &\\;\\;\\; E[(\\eta - R_w)_+]\\leq E[(\\eta -Y)_+] \\ & \\ \\ \\text{for all} \\ \\eta\\in [a,b]\\\\ \n",
                "   &\\;\\;\\;  w \\in W. \\nonumber\n",
                "\\end{align}"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "8b1ce13e",
            "metadata": {},
            "source": [
                "### 2.3. Show that we don't need to compare intergals of CDFs for all $\\eta \\in [a,b]$\n",
                "\n",
                "Further assume Y has a discrete distribution with realizations $y_i$, $i =1,...,m$ then Dentcheva and Ruszczyński (2003) shows:\n",
                "\n",
                "\\begin{aligned}\n",
                "E[(\\eta - R_w)_+]\\leq E[(\\eta -Y)_+] \\ & \\ \\ \\text{for all} \\ \\eta\\in [a,b]\n",
                "\\end{aligned}\n",
                "\n",
                "$$\\Longleftrightarrow$$\n",
                "\n",
                "\\begin{aligned}\n",
                "E[(y_i - R_w)_+]\\leq E[(y_i -Y)_+], \\ & \\ \\ \\       i=1,...,m \n",
                "\\end{aligned}"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "745b12ba",
            "metadata": {},
            "source": [
                "This equivalent relationship allow us to compare these partial moments only for a finite number of return levels."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "a642cee6",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "source": [
                "### 2.4. Further transform the problem into a linear programming problem\n",
                "\n",
                "Recall that $R_w: \\Omega \\rightarrow \\mathbb{R}$. Assume there are finitly many elementary events such that $\\Omega=\\{\\nu_1,...,\\nu_m\\}$ and introduce a decision vector $S:[a,b] \\times \\Omega \\rightarrow  \\mathbb{R}$ , then the lower partial moment $E[(\\eta - R_w)_+]$ for an asset allocation  $𝑤\\in W$ will be the solution of the following linear programming problem:\n",
                "\n",
                "\\begin{align}\n",
                "\\min_{S(\\eta,\\nu)} &\\;\\;\\;E[ S(\\eta,\\nu)] \\nonumber \\\\\n",
                "\\text{subject to} &\\;\\;\\; S(\\eta,\\nu_i) \\geq \\eta - R_w(\\nu_i) &\\text{for} \\;i=1,...,m \\nonumber \\\\\n",
                "&\\;\\;\\; S(\\eta,\\nu_i) \\geq 0 &\\text{for} \\;i=1,...,m\n",
                "\\end{align}\n",
                "\n",
                "To understand why, it's important to note that the expectation $E[ S(\\eta,\\nu)]$ can be written as the sum of all $S(\\eta,\\nu_i)$ multiplied by their respective probabilities, that is  $E[ S(\\eta,\\nu)] = P[\\{\\nu_1\\}] S(\\eta,\\nu_1) + ... + P[\\{\\nu_m\\}] S(\\eta,\\nu_m)$. Therefore, minimizing $E[ S(\\eta,\\nu)]$ is equivalent to minimizing $S(\\eta,\\nu_i)$ for all $\\nu_i \\in \\Omega$ under the current setting. As $S(\\eta,\\nu_i)$ is greater than or equal to both $\\eta - R_w(\\nu_i)$ and zero,  minimizing $S(\\eta,\\nu_i)$ will enforce the solution $S(\\eta,\\nu_i)^*$ equal to $max\\{ \\eta - R_w(\\nu_i),0\\}$. Therefore $E[ S(\\eta,\\nu)^*]=E[(\\eta - R_w)_+]$. Using this idea, the original problem could be further transformed into the following:\n",
                "\n",
                "\\begin{align}\n",
                "\\max  &\\;\\;\\; E[R_w]\\\\\n",
                "\\text{subject to}  &\\;\\;\\; S(y_i,\\nu_k)\\geq y_i -R_w(\\nu_i) &i=1,...,m, k=1,..., m\\nonumber\\\\\n",
                "&\\;\\;\\; E[S(y_i,\\nu)]\\leq E[(y_i-Y)_+] &i=1,...,m\\nonumber\\\\\n",
                "&\\;\\;\\;S(y_i,\\nu_k)\\geq 0 &i=1,...,m, k=1,..., m\\nonumber\\\\\n",
                "&\\;\\;\\;w \\in W \\nonumber\\\\\n",
                "\\end{align}\n",
                "\n",
                "Kindly note that in this problem, the solution $S(η,\\nu_i)^{**}$ will not necessarily enforce $E[S(y_i,\\nu)^{**}]$ to be equal to $E[(y_i-R_w)_+]$, since we are not explicitly minimizing $E[S(\\eta,\\nu)]$. However, we know that $E[(y_i-R_w)_+] \\leq E[S(y_i,\\nu)^{**}]$, therefore solving the above problem will ensure that $E[(y_i-R_w)_+] \\leq E[(y_i-Y)_+]$.\n"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "5fd477e6",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "source": [
                "Before we show the final product, we need some more notations.  \n",
                "\n",
                "- Let $p_k=P[\\{\\nu_k\\}]$, $r_{nk}=R_n(\\nu_k), s_{ik}=S(y_i,\\nu_k)$, and $v_i=E[(y_i-Y)_+]$.\n",
                "\n",
                "Then the original problem could be formulated as the following linear programming problem:\n",
                "\n",
                "\\begin{aligned}\n",
                "\\max_{w,s} &\\; \\sum_{k=1}^{m}\\sum_{n=1}^{N} p_kr_{nk}w_n\\nonumber\\\\\n",
                "\\text{subject to}&\\; \\sum_{n=1}^{N} r_{nk}w_n+s_{ik}\\geq y_i, &i=1,...,m, k=1,..., m\\nonumber\\\\\n",
                "&\\sum_{k=1}^{m} p_k s_{ik}\\leq v_i, &i=1,...,m\\nonumber\\\\\n",
                "&s_{ik}\\geq0&i=1,...,m, k=1,..., m\\nonumber\\\\\n",
                "&\\sum^N_{n=1} w_n =1\\nonumber\\\\\n",
                "&w_n \\geq 0 &n=1,...,N \\nonumber\\\\\n",
                "\\end{aligned}"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "9fdb6fe7",
            "metadata": {},
            "source": [
                "## 3. An example in Dentcheva and Ruszczyński (2003)  with Python\n",
                "\n",
                "Now we turn to the an example at p.562 in Dentcheva and Ruszczyński (2003). We will use the following imports:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "id": "403e64e8",
            "metadata": {
                "slideshow": {
                    "slide_type": "slide"
                }
            },
            "outputs": [],
            "source": [
                "import numpy as np\n",
                "import pandas as pd\n",
                "import itertools\n",
                "import matplotlib.pyplot as plt\n",
                "import time\n",
                "\n",
                "from scipy.optimize import linprog\n",
                "from matplotlib.patches import Polygon\n",
                "from itertools import repeat\n",
                "%matplotlib inline"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "0ef69616",
            "metadata": {},
            "source": [
                "In this notebook, we use the linear programming solver in SciPy. Sometimes we also use solver in CVXOPT and CVXPY which have different advantages. \n",
                "\n",
                "Now let's take a look on the data:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "id": "8843c5a0",
            "metadata": {
                "slideshow": {
                    "slide_type": "slide"
                }
            },
            "outputs": [],
            "source": [
                "#Number of w is N=8 and s = m*m = 22*22 = 484 is 492\n",
                "# intialise data of lists.\n",
                "\n",
                "data = {'Asset 1':[  7.5,   8.4,  6.1,  5.2,  5.5,\n",
                "                     7.7,  10.9, 12.7, 15.6, 11.7,\n",
                "                     9.2,  10.3,    8,  6.3,  6.1,\n",
                "                     7.1,   8.7,    8,  5.7,  3.6,  3.1,  4.5],\n",
                "        'Asset 2':[ -5.8,     2,  5.6, 17.5,  0.2, \n",
                "                    -1.8,  -2.2, -5.3,  0.3, 46.5,\n",
                "                    -1.5,  15.9, 36.6, 30.9, -7.5,\n",
                "                     8.6,  21.2,  5.4, 19.3,  7.9, 21.7,-11.1],\n",
                "        'Asset 3':[-14.8, -26.5, 37.1, 23.6, -7.4,\n",
                "                     6.4,  18.4, 32.3, -5.1, 21.5,\n",
                "                    22.4,   6.1, 31.6, 18.6,  5.2,\n",
                "                    16.5,  31.6, -3.2, 30.4,  7.6,   10,  1.2],\n",
                "        'Asset 4':[-18.5, -28.4, 38.5, 26.6, -2.6,\n",
                "                     9.3,  25.6, 33.7, -3.7, 18.7,\n",
                "                    23.5,     3, 32.6, 16.1,  2.3,\n",
                "                    17.9,  29.2, -6.2, 34.2,    9, 11.3, -0.1],\n",
                "        'Asset 5':[-30.2, -33.8, 31.8,   28,  9.3, \n",
                "                    14.6,  30.7, 36.7,   -1, 21.3,\n",
                "                    21.7,  -9.7, 33.3,  8.6, -4.1,\n",
                "                    16.5,  20.4,  -17, 59.4, 17.4, 16.2, -3.2],\n",
                "        'Asset 6':[  2.3,   0.2, 12.3, 15.6,    3,  \n",
                "                     1.2,   2.3,  3.1,  7.3, 31.1,\n",
                "                       8,    15, 21.3, 15.6,  2.3,\n",
                "                     7.6,  14.2,  8.3, 16.1,  7.6,   11, -3.5],\n",
                "        'Asset 7':[-14.9, -23.2, 35.4,  2.5, 18.1, \n",
                "                    32.6,   4.8, 22.6, -2.3, -1.9,\n",
                "                    23.7,   7.4, 56.2, 69.4, 24.6,\n",
                "                    28.3,  10.5,-23.4, 12.1,-12.2, 32.6,  7.8],\n",
                "        'Asset 8':[ 67.7,  72.2,  -24,   -4,   20, \n",
                "                    29.5,  21.2, 29.6,-31.2,  8.4,\n",
                "                   -12.8, -17.5,  0.6, 21.6, 24.4,\n",
                "                   -13.9,  -2.3, -7.8, -4.2, -7.4, 14.6,   -1]}"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 3,
            "id": "ccd41972",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "outputs": [],
            "source": [
                "# Create DataFrame\n",
                "df = pd.DataFrame(data)\n",
                "\n",
                "# Number of assets\n",
                "NumA = len(df.columns)\n",
                "# Number of states\n",
                "NumS = len(df)\n",
                "\n",
                "# List for col index\n",
                "cols = [f'Asset {i}' for i in range(1, NumA+1)]"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "f989c615",
            "metadata": {},
            "source": [
                "We start by setting up the target portfolio. It is just a equally weighted portfolio."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 4,
            "id": "792aa23a",
            "metadata": {},
            "outputs": [],
            "source": [
                "# Create an array of weights, each set to 1/NumA\n",
                "weight = np.full(NumA, 1/NumA)\n",
                "\n",
                "# Calculate the weighted average of the columns in the dataframe\n",
                "y_k = np.average(np.asarray(df[cols]), weights=weight, axis=1)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "9b3f78d4",
            "metadata": {},
            "source": [
                "### 3.1. The goal\n",
                "\n",
                "The goal is use the following line at the end:\n",
                "\n",
                "res_ex1 = linprog(-c_ex1, A_ub=A_ex1, b_ub=b_ex1, A_eq=A_ex2, b_eq=b_ex2,\n",
                "                  bounds=bounds_ex1, method='revised simplex')\n",
                "\n",
                "where \"-c_ex1\" is the coefficients of the objective function, \"A_ex1\" is the inequality constraint matrix, \"b_ex1\" is related to the upper bound of the inequality constraint matrix, \"A_ex2\" is the equality constraint matrix, \"b_ex2\" is the equality constraint vector, and \"bounds_ex1\" defines the minimum and maximum of decision variables.\n",
                "\n",
                "For more details, read: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html\n"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "ff457b2a",
            "metadata": {},
            "source": [
                "### 3.2. Pain in the ass (or Setting up objective and constraints)\n",
                "\n",
                "Here we setup the parameters for the objective function:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 5,
            "id": "f36b9d96",
            "metadata": {},
            "outputs": [],
            "source": [
                "# Calculate the probability of each state \n",
                "p_k = 1/NumS\n",
                "\n",
                "# Calculate the probability of each state times the total return of assets\n",
                "#e.g pk(r11+r21+...), pk(r21+r22+...)\n",
                "pk_S_rnk = p_k * df[cols].sum()\n",
                "\n",
                "# Convert the resulting series to a list\n",
                "pk_S_rnk = pk_S_rnk.values.tolist()\n",
                "\n",
                "#The s_ik variables\n",
                "s_ik = [0]*NumS*NumS\n",
                "\n",
                "# Setting the objective for linear programing \n",
                "# Adding the list together\n",
                "c_ex1 = np.array(pk_S_rnk + s_ik)  "
            ]
        },
        {
            "cell_type": "markdown",
            "id": "274cca69",
            "metadata": {},
            "source": [
                "Now we turn to inequality constraints:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 6,
            "id": "f3e3730e",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "outputs": [],
            "source": [
                "r_list = df.values.tolist()\n",
                "\n",
                "#Copy the return for m states\n",
                "r_list = r_list * NumS \n",
                "\n",
                "# Stack two below matrices together\n",
                "b1 = np.array(r_list)\n",
                "b2 = np.identity(NumS*NumS)\n",
                "\n",
                "# The first line of the constraints\n",
                "A_ex1_1 = np.column_stack((b1, b2))"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 7,
            "id": "1e9cc2b8",
            "metadata": {},
            "outputs": [],
            "source": [
                "#b3 is for w_i\n",
                "b3 = np.array(df.values.tolist())*0\n",
                "b4 = np.zeros((NumS,NumS*NumS))\n",
                "\n",
                "for x in range(NumS):\n",
                "   b4[x,NumS*x:NumS*(x+1)] = 1  \n",
                "    \n",
                "A_ex1_2 = p_k*np.column_stack((b3, b4))\n",
                "\n",
                "# The inequality constraints for linear programing\n",
                "A_ex1 = np.row_stack((-A_ex1_1,A_ex1_2))"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 8,
            "id": "a8bf4c23",
            "metadata": {
                "slideshow": {
                    "slide_type": "slide"
                }
            },
            "outputs": [],
            "source": [
                "#Working on the bounds of inequality constraints\n",
                "#v_i is patial moments\n",
                "v_i = np.array([])\n",
                "\n",
                "\n",
                "for x in range(NumS):\n",
                "    yp = y_k[x]-y_k\n",
                "    yp[yp < 0] = 0\n",
                "    v_i = np.append(v_i, np.mean(yp))\n",
                "\n",
                "\n",
                "# Bounds of inequality constraints used for lP\n",
                "b_ex1 = np.array([])\n",
                "for i in range(NumS):\n",
                "    for x in range(NumS):\n",
                "        b_ex1 = np.append(b_ex1, y_k[i])\n",
                "\n",
                "\n",
                "b_ex1 = np.append(-b_ex1,v_i)"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "1abac7fd",
            "metadata": {},
            "source": [
                "Last but not least, the equality constraint:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 9,
            "id": "3c4316de",
            "metadata": {
                "slideshow": {
                    "slide_type": "subslide"
                }
            },
            "outputs": [],
            "source": [
                "# Working on equality constraint\n",
                "A_ex2 = np.array([np.append(np.ones(NumA), np.zeros(NumS*NumS))])\n",
                "\n",
                "#print(A_ex2)\n",
                "b_ex2 = np.array([1])\n",
                "\n",
                "# Bounds on decision variables\n",
                "bounds_ex1 =[]\n",
                "for x in range(NumA+NumS*NumS):\n",
                "    bounds_ex1 = bounds_ex1 + [(0, None)]"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "297819bb",
            "metadata": {},
            "source": [
                "### 3.3. Solving the linear programming problem\n",
                "Now, we can use the solver to figure out the best solution of our linear programming problem:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 10,
            "id": "553151e5",
            "metadata": {
                "slideshow": {
                    "slide_type": "slide"
                }
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "The optimal portfolio has the form\n",
                        "[0.         0.         0.06803595 0.18800344 0.         0.39137554\n",
                        " 0.23092417 0.12166089]\n"
                    ]
                }
            ],
            "source": [
                "# Solve the problem\n",
                "\n",
                "res_ex1 = linprog(-c_ex1, A_ub=A_ex1, b_ub=b_ex1, A_eq=A_ex2, b_eq=b_ex2,\n",
                "                  bounds=bounds_ex1, method='revised simplex')\n",
                "\n",
                "print(\"The optimal portfolio has the form\")\n",
                "print(res_ex1.x[0:8])"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "bb728022",
            "metadata": {},
            "source": [
                "## Reference\n",
                "\n",
                "Atkinson, A. B. (1970). On the measurement of inequality. Journal of economic theory, 2(3), 244-263.\n",
                "\n",
                "Bawa, V. S., Bodurtha Jr, J. N., Rao, M. R., & Suri, H. L. (1985). On determination of stochastic dominance optimal sets. The Journal of Finance, 40(2), 417-431.\n",
                "\n",
                "Cho, Y. H., Linton, O., & Whang, Y. J. (2007). Are there Monday effects in stock returns: A stochastic dominance approach. Journal of Empirical Finance, 14(5), 736-755.\n",
                "\n",
                "Chui, D., Cheng, W. W., Chow, S. C., & Ya, L. I. (2020). Eastern Halloween effect: A stochastic dominance approach. Journal of International Financial Markets, Institutions and Money, 68, 101241.\n",
                "\n",
                "Dentcheva, D., & Ruszczyński, A. (2003). Optimization with stochastic dominance constraints. SIAM Journal on Optimization, 14(2), 548-566.\n",
                "\n",
                "Fang, Y., & Post, T. (2022). Optimal portfolio choice for higher-order risk averters. Journal of Banking & Finance, 137, 106429.\n",
                "\n",
                "Guerre, E., Perrigne, I., & Vuong, Q. (2009). Nonparametric identification of risk aversion in first‐price auctions under exclusion restrictions. Econometrica, 77(4), 1193-1227.\n",
                "\n",
                "Hadar, J., & Russell, W. R. (1969). Rules for ordering uncertain prospects. The American economic review, 59(1), 25-34.\n",
                "\n",
                "Hanoch, G., & Levy, H. (1975). The efficiency analysis of choices involving risk. In Stochastic Optimization Models in Finance (pp. 89-100). Academic Press.\n",
                "\n",
                "Levy, H. (2016). Stochastic dominance: Investment decision making under uncertainty. New York: Springer.\n",
                "\n",
                "Maasoumi, E., Millimet, D. L., & Rangaprasad, V. (2005). Class size and educational policy: Who benefits from smaller classes?. Econometric Reviews, 24(4), 333-368.\n",
                "\n",
                "Rothschild, M., & Stiglitz, J. E. (1978). Increasing risk: I. A definition. In Uncertainty in Economics (pp. 99-121). Academic Press.\n",
                "\n",
                "Tesfatsion, L. (1976). Stochastic dominance and the maximization of expected utility. The Review of Economic Studies, 43(2), 301-315."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "22fa1f0e",
            "metadata": {},
            "source": [
                "### $\\bigstar$ I want to thank John Stachurski for inviting me to write this note. I would also like to thank John, Simon Grant and Timothy Kam for their valuable comments. Any errors that appear in this note are my sole responsibility."
            ]
        }
    ],
    "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
}