portfolio_selection_optimization.ipynb.txt

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "# Portfolio Selection Optimization\n",
    "This model is an example of the classic [Markowitz portfolio selection optimization model](https://en.wikipedia.org/wiki/Markowitz_model). We want to find the fraction of the portfolio to invest among a set of stocks that balances risk and return. It is a Quadratic Programming (QP) model with vector and matrix data for returns and risk, respectively. This is best suited to a matrix formulation, so we use the Gurobi Python *matrix* interface. The basic model is fairly simple, so we also solve it parametrically to find the efficient frontier.\n",
    "\n",
    "**Download the Repository** <br /> \n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). \n",
    "\n",
    "**Gurobi License** <br /> \n",
    "In order to run this Jupyter Notebook properly, you must have a Gurobi license. If you do not have one, you can request an [evaluation license](https://www.gurobi.com/downloads/request-an-evaluation-license/?utm_source=3PW&utm_medium=OT&utm_campaign=WW-MU-MUI-OR-O_LEA-PR_NO-Q3_FY20_WW_JPME_Lost_Luggage_Distribution_COM_EVAL_GitHub&utm_term=Lost%20Luggage%20Distribution&utm_content=C_JPM) as a *commercial user*, or download a [free license](https://www.gurobi.com/academia/academic-program-and-licenses/?utm_source=3PW&utm_medium=OT&utm_campaign=WW-MU-EDU-OR-O_LEA-PR_NO-Q3_FY20_WW_JPME_Lost_Luggage_Distribution_COM_EVAL_GitHub&utm_term=Lost%20Luggage%20Distribution&utm_content=C_JPM) as an *academic user*.\n",
    "\n",
    "\n",
    "## Model Formulation\n",
    "### Parameters\n",
    "\n",
    "We use the [Greek values](https://en.wikipedia.org/wiki/Greeks_\\(finance\\)) that are traditional in finance:\n",
    "\n",
    "- $\\delta$: n-element vector measuring the change in price for each stock\n",
    "- $\\sigma$: n x n matrix measuring the covariance among stocks\n",
    "\n",
    "There is one additional parameter when solving the model parametrically:\n",
    "\n",
    "- r: target return\n",
    "\n",
    "\n",
    "### Decision Variables\n",
    "- $x \\ge 0$: n-element vector where each element represents the fraction of the porfolio to invest in each stock\n",
    "\n",
    "### Objective Function\n",
    "Minimize the total risk, a convex quadratic function:\n",
    "\n",
    "\\begin{equation}\n",
    "\\min x^t \\cdot \\sigma \\cdot x\n",
    "\\end{equation}\n",
    "\n",
    "### Constraints\n",
    "\n",
    "Allocate the entire portfolio: the total investments should be 1.0 (100%), where $e$ is a unit vector (all 1's):\n",
    "\n",
    "\\begin{equation}\n",
    "e \\cdot x = 1\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Return: When we solve the model parametrically for different return values $r$, we add a constraint on the target return:\n",
    "\n",
    "\\begin{equation}\n",
    "\\delta \\cdot x = r\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Python Implementation\n",
    "### Stock data\n",
    "Use [yfinance](https://pypi.org/project/yfinance/) library to get the latest 2 years of _actual stock data_ from the 20 most profitable US companies, [according to Wikipedia in April 2021](https://en.wikipedia.org/wiki/List_of_largest_companies_in_the_United_States_by_revenue#List_of_companies_by_profit)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import yfinance as yf\n",
    "\n",
    "stocks = ['BRK-A', 'AAPL', 'MSFT', 'JPM', 'GOOG', 'BAC', 'INTC', 'WFC',\n",
    "          'C', 'VZ', 'FB', 'PFE', 'JNJ', 'WMT', 'XOM',\n",
    "          'FNMA', 'T', 'UNH', 'CMCSA', 'V' ]\n",
    "\n",
    "data = yf.download(stocks, period='2y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute Greeks\n",
    "Using the downloaded stock data, find the delta (return), sigma (covariance) and standard deviation values for stock prices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "closes = np.transpose(np.array(data.Close)) # matrix of daily closing prices\n",
    "absdiff = np.diff(closes)                   # change in closing price each day\n",
    "reldiff = np.divide(absdiff, closes[:,:-1]) # relative change in daily closing price\n",
    "delta = np.mean(reldiff, axis=1)            # mean price change\n",
    "sigma = np.cov(reldiff)                     # covariance (standard deviations)\n",
    "std = np.std(reldiff, axis=1)               # standard deviation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Minimize risk by solving QP model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "from math import sqrt\n",
    "\n",
    "# Create an empty model\n",
    "m = gp.Model('portfolio')\n",
    "\n",
    "# Add matrix variable for the stocks\n",
    "x = m.addMVar(len(stocks))\n",
    "\n",
    "# Objective is to minimize risk (squared).  This is modeled using the\n",
    "# covariance matrix, which measures the historical correlation between stocks\n",
    "portfolio_risk = x @ sigma @ x\n",
    "m.setObjective(portfolio_risk, GRB.MINIMIZE)\n",
    "\n",
    "# Fix budget with a constraint\n",
    "m.addConstr(x.sum() == 1, 'budget')\n",
    "\n",
    "# Verify model formulation\n",
    "m.write('portfolio_selection_optimization.lp')\n",
    "\n",
    "# Optimize model to find the minimum risk portfolio\n",
    "m.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display minimum risk portfolio using Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "minrisk_volatility = sqrt(m.ObjVal)\n",
    "minrisk_return = delta @ x.X\n",
    "pd.DataFrame(data=np.append(x.X, [minrisk_volatility, minrisk_return]),\n",
    "             index=stocks + ['Volatility', 'Expected Return'],\n",
    "             columns=['Minimum Risk Portfolio'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the efficient frontier\n",
    "Solve the QP parametrically to find the lowest risk portfolio for different expected returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an expression representing the expected return for the portfolio\n",
    "portfolio_return = delta @ x\n",
    "target = m.addConstr(portfolio_return == minrisk_return, 'target')\n",
    "\n",
    "# Solve for efficient frontier by varying target return\n",
    "frontier = np.empty((2,0))\n",
    "for r in np.linspace(delta.min(), delta.max(), 25):\n",
    "    target[0].rhs = r\n",
    "    m.optimize()\n",
    "    frontier = np.append(frontier, [[sqrt(m.ObjVal)],[r]], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot results\n",
    "Use the matplot library to plot the optimized solutions, along with the individual stocks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "#plt.figure(figsize=(10,10))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,8))\n",
    "\n",
    "# Plot volatility versus expected return for individual stocks\n",
    "ax.scatter(x=std, y=delta,\n",
    "           color='Blue', label='Individual Stocks')\n",
    "for i, stock in enumerate(stocks):\n",
    "    ax.annotate(stock, (std[i], delta[i]))\n",
    "\n",
    "# Plot volatility versus expected return for minimum risk portfolio\n",
    "ax.scatter(x=minrisk_volatility, y=minrisk_return, color='DarkGreen')\n",
    "ax.annotate('Minimum\\nRisk\\nPortfolio', (minrisk_volatility, minrisk_return),\n",
    "            horizontalalignment='right')\n",
    "\n",
    "# Plot efficient frontier\n",
    "ax.plot(frontier[0], frontier[1], label='Efficient Frontier', color='DarkGreen')\n",
    "\n",
    "# Format and display the final plot\n",
    "ax.axis([frontier[0].min()*0.7, frontier[0].max()*1.3, delta.min()*1.2, delta.max()*1.2])\n",
    "ax.set_xlabel('Volatility (standard deviation)')\n",
    "ax.set_ylabel('Expected Return')\n",
    "ax.legend()\n",
    "ax.grid()\n",
    "plt.show()"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}