09-Markowitz-Portfolio-Selection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pickle\n",
    "import cvxpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to Portfolio Theory\n",
    "\n",
    "In classical portfolio theory periodic returns are normally distributed. We will assume daily returns of two stocks $R_A$ and $R_B$ denote prices of two stocks $A$ and $B$, respectively, such that they are jointly normally distributed with the following properties. Marginally, \n",
    "$$R_A \\sim N(\\mu_A,\\sigma_A^2) \\text{, and } R_B \\sim N(\\mu_B,\\sigma_B^2),$$\n",
    "and, jointly, their covariance and correlation are\n",
    "$$\\sigma_\\text{AB} = \\text{Cov}(R_A,R_B) \\text{, and } \\rho_\\text{AB} = \\frac{\\sigma_\\text{AB}}{\\sigma_A\\sigma_B}$$\n",
    "\n",
    "If we consider managing a portfolio consisting of these two stocks, then the return for the portfolio $R_p$ is \n",
    "$$R_p = R_A x_A + R_B x_B,$$\n",
    "where $x_A$ and $x_B$ are fixed proportions of portfolio such that $x_A + x_B = 1$.\n",
    "\n",
    "What are the properties of $R_p$? We know $R_p$ is also normally distributed with following mean and variance:\n",
    "$$ \n",
    "\\mu_p = \\text{E}(R_p) = \\text{E}(R_B)\\, x_A + \\text{E}(R_B)\\, x_B = \\mu_A x_A + \\mu_B x_B\\\\\n",
    "\\sigma_p^2 = \\text{Var}(R_p) = \\text{E}\\left((R_p - \\text{E}(R_p))^2\\right) = \\sigma_A^2\\,x_A^2 + \\sigma_B^2\\,x_B^2 + 2\\sigma_{AB}x_A x_B.\n",
    "$$\n",
    "Note that $\\sigma_{AB} = \\rho_{AB}\\sigma_A\\sigma_B$.\n",
    "\n",
    "Finally, the distribution of the portfolio return is\n",
    "$$ R_p \\sim N(\\mu_p, \\sigma_p^2) $$\n",
    "\n",
    "Let's simulate behavior of returns $R_A$, $R_B$, and $R_p$ (portfolio weight distributed evenly): i.e. $x_A = x_B = 0.5$. Also, choose some parameters for our two stocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = pd.Series({'A': 0.5, 'B': 0.5})\n",
    "mu = pd.Series({'A':0.15, 'B':0.1})\n",
    "sig = pd.Series({'A':0.1, 'B':0.05})\n",
    "\n",
    "## correlation between stocks\n",
    "rhoAB = -0.15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate theoretical portfolio returns\n",
    "\n",
    "Calculate the volatility ($\\sigma^2$) and expected returns ($\\mu$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu['p'] = mu.A * x.A + mu.B * x.B\n",
    "sig['p'] = np.sqrt((sig.A**2 * x.A**2) + (sig.B**2 * x.B**2) + (2 * rhoAB * sig.A *sig.B * x.A * x.B))\n",
    "\n",
    "[one for one in zip(sig**2,mu)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate random returns\n",
    "\n",
    "Generate random observations from bivariate normal distribution specified by "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "param_mu = [mu.A, mu.B]\n",
    "param_sigma = [[sig.A**2, rhoAB*sig.A*sig.B],\n",
    "               [rhoAB*sig.A*sig.B, sig.B**2]]\n",
    "\n",
    "R_AB = np.random.multivariate_normal(param_mu, param_sigma, 1000)\n",
    "RA = R_AB[:,0]\n",
    "RB = R_AB[:,1]\n",
    "\n",
    "print('RA: mean is', RA.mean(), ', std dev is', RA.std())\n",
    "print('RB: mean is', RB.mean(), ', std dev is', RB.std())\n",
    "print('Correlation between RA and RB is', np.corrcoef(RA, RB)[1,0])\n",
    "\n",
    "pd.tools.plotting.scatter_matrix(pd.DataFrame(R_AB, columns=['RA', 'RB']));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Volatility of returns are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $\\mu_p$ is a simple average of the returns; however, $\\sigma_p^2 < 0.5\\sigma_A^2 + 0.5\\sigma_B^2$. This means that volatility is less than if $A$ and $B$ were independent. This owes to the characteristic that $\\rho<0$, so that the cross-term subtracts from $\\sigma_p^2$. In fact, this is theoretical basis for the advantage of diversification."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Efficient frontier of portfolios\n",
    "\n",
    "What happens if we change the proportion of the portfolios $x_A$ and $x_B$? For a range of possibilities, we compute the resulting return $\\mu_p$ and volatility $\\sigma_p^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xa = np.linspace(-1, 2, num=101)\n",
    "pf = pd.DataFrame({'xa':xa, 'xb':1-xa})\n",
    "\n",
    "pf['mu'] = mu.A * pf.xa + mu.B * pf.xb\n",
    "pf['sig2'] = ((sig.A**2 * pf.xa**2) + (sig.B**2 * pf.xb**2) + (2 * rhoAB * sig.A *sig.B * pf.xa * pf.xb))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Long and short positions\n",
    "\n",
    "Note that `xa` is allowed to be negative. A positive $x_A$ value indicates what is called a _long position_ for stock A and indicates a purchase. On the otherhand, a negative $x_A$ indicates a _short position_ stock A and indicates a sale. The sale of A allows for higher long positions on other stocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = plt.subplot(1,1,1)\n",
    "plt.plot(pf.sig2, pf.mu, color='b')\n",
    "plt.plot(pf.sig2[pf.xa<0], pf.mu[pf.xa<0], color='r')\n",
    "plt.plot(pf.sig2[pf.xb<0], pf.mu[pf.xb<0], color='g')\n",
    "plt.scatter(sig**2, mu)\n",
    "p.set_xlabel('Portfolio volatility')\n",
    "p.set_ylabel('Portfolio return')\n",
    "plt.show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above plot is called the efficient frontier. The line represents the lower boundary of achievable volatility for each level of portfolio return using component stocks. Two stocks $A$ and $B$ are shown as dots in the interior of the parabola, and the dot on the parabola is the return and volatility of a portfolio in which $x_A = x_b = 0.5$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Minimum variance portfolio (MVP)\n",
    "\n",
    "What is the point at the tip of the parabola? This point indicates what is called the minimum variance portfolio that, regardless of the return, has the minimum volatility. Minimum variance portfolio (MVP) has the minimum volatility of any portfolio achievable.\n",
    "\n",
    "It can be shown that minimum variance portfolio would have the following allocation $x_A$:\n",
    "\n",
    "$$x_A = \\frac{\\sigma_B^2-\\sigma_\\text{AB}^2}{\\sigma_A^2 + \\sigma_B^2 – \\sigma_\\text{AB}^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mvp_A = (sig.B**2 - rhoAB * sig.A *sig.B )/(sig.A**2 + sig.B**2 - rhoAB * sig.A *sig.B ) \n",
    "mvp_B = 1 - mvp_A\n",
    "mvp_mu = mu.A * mvp_A + mu.B * mvp_B \n",
    "mvp_sig2 = ((sig.A**2 * mvp_A**2) + (sig.B**2 * mvp_B**2) + (2 * rhoAB * sig.A *sig.B * mvp_A * mvp_B))\n",
    "mvp_mu, mvp_sig2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the minimum variance portfolio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = plt.subplot(1,1,1)\n",
    "plt.plot(pf.sig2, pf.mu, color='b')\n",
    "plt.plot(pf.sig2[pf.xa<0], pf.mu[pf.xa<0], color='r')\n",
    "plt.plot(pf.sig2[pf.xb<0], pf.mu[pf.xb<0], color='g')\n",
    "plt.scatter(mvp_sig2, mvp_mu, color='r', marker='+')\n",
    "p.set_xlabel('Portfolio volatility')\n",
    "p.set_ylabel('Portfolio return')\n",
    "plt.show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mathematical problem to solve for the minimum variance portfolio can be stated as follows:\n",
    "$$ \\min_{x_A,x_B}\\ \\ \\sigma_p^2 = \\sigma_A^2 x_A^2 + \\sigma_B^2 x_B^2 + \\sigma_{AB} x_A x_B\\\\\n",
    "\\text{such that }x_A + x_B = 1$$\n",
    "Here, $\\sigma_p^2 = \\sigma_A^2 x_A^2 + \\sigma_B^2 x_B^2 + \\sigma_{AB} x_A x_B$ is called the objective function, and $x_A + x_B = 1$ is called the constraint.\n",
    "\n",
    "The constraint has to be satisfied exactly, and objective wants to be minimized with the freedom of changing the values of $x_A$ and $x_B$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `cvxpy` for numerical optimization\n",
    "\n",
    "Calculating analytical solutions is not always possible. There are ways to compute solutions numerically. Although out of scope of our class, a large field in applied math called optimization. The tools from this area allow us to compute solutions for problems such as the one we have here. (For more information, visit the "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Higher dimensional portfolios\n",
    "\n",
    "In higher dimensions, the portfolio equantion is written in terms of vectors and matrices. For example, suppose we construct a portfolio consisting of $s$-assets. The portfolio allocation can be written as,\n",
    "$$\\mathbf{1}^\\intercal x = 1,$$\n",
    "where $x$ is an $s$-vector whose sum adds up to 1: i.e., $x_1 + x_2 + \\cdots + x_s = 1$.\n",
    "\n",
    "The portfolio volatility in matrix-vector form looks like:\n",
    "$$\n",
    "\\begin{pmatrix}x_1 & x_2 & \\cdots & x_s\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\sigma_1^2 & \\sigma_{12} & \\cdots & \\sigma_{1s}\\\\\n",
    "\\sigma_{21} & \\sigma_2^2 & \\cdots & \\sigma_{2s}\\\\\n",
    "\\vdots & \\vdots & & \\vdots \\\\\n",
    "\\sigma_{s1} & \\sigma_{s2} & \\cdots & \\sigma_s^2\\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_s\\end{pmatrix}.\n",
    "$$\n",
    "To verify the two stock example result, we check $\\sigma_p^2$:\n",
    "$$\n",
    "\\sigma_p^2 = \\begin{pmatrix}x_1 & x_2\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\sigma_1^2 & \\sigma_{12}\\\\\n",
    "\\sigma_{21} & \\sigma_2^2\\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}x_1 \\\\ x_2\\end{pmatrix} = \n",
    "\\begin{pmatrix}x_1 & x_2\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\sigma_1^2 x_1 + \\sigma_{12} x_2\\\\\n",
    "\\sigma_{21} x_1 + \\sigma_2^2 x_2\\\\\n",
    "\\end{pmatrix} = \n",
    "\\sigma_1^2 x_1^2 + \\sigma_2^2 x_2^2 + \\sigma_{12} x_1 x_2,\n",
    "$$\n",
    "which matches the result from before.\n",
    "\n",
    "With this notation, we can re-write the problem.\n",
    "\n",
    "$$ \\min_{x\\in\\mathbb{R}^s}\\ \\ x^\\intercal \\Sigma x\\\\\n",
    "\\text{subject to }\\mathbf{1}^\\intercal x = 1,$$\n",
    "\n",
    "It can be shown that the optimal solution to this problem is,\n",
    "$$ x^* = (\\mathbf{1}^\\intercal\\Sigma^{-1}\\mathbf{1})^{-1}\\Sigma^{-1}\\mathbf{1} $$.\n",
    "\n",
    "There are many variations of this problem. In fact the minimum variance portfolio is special case of a more general problem\n",
    "\n",
    "$$ \\min_{x\\in\\mathbb{R}^s}\\ \\ x^\\intercal \\Sigma x\\\\\n",
    "\\text{subject to }\\mu^\\intercal x\\geq \\mu^* \\text{, and } \\mathbf{1}^\\intercal x = 1,$$\n",
    "\n",
    "Note that in all of the above problems, the measure of $\\mu^*$, expected returns, and covariance matrix have to be provided. In practice, these quantities are estimated from data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download data from Quandl\n",
    "\n",
    "Use the Quandl API to download some stock data. In this example, we will use the constituent stocks from the [Dow Jones Industrial Average](https://en.wikipedia.org/wiki/Dow_Jones_Industrial_Average) to form our portfolio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import mykeys as m\n",
    "# import quandl\n",
    "# \n",
    "# quandl.ApiConfig.api_key = None\n",
    "# quandl.ApiConfig.api_key = m.apikeys['quandl']\n",
    "# \n",
    "# # get the table for daily stock prices and,\n",
    "# # filter the table for selected tickers, columns within a time range\n",
    "# # set paginate to True because Quandl limits tables API to 10,000 rows per call\n",
    "# \n",
    "# symbols = ['AAPL','AXP','BA','BAC','CAT',\n",
    "#            'CSCO','CVX','DD','DWDP','DIS','GE',\n",
    "#            'HD','HPQ','IBM','INTC','JNJ',\n",
    "#            'JPM','KFT','KO','MCD','MMM',\n",
    "#            'MRK','MSFT','PFE','PG','T',\n",
    "#            'TRV','UTX','VZ','WMT','XOM']\n",
    "# \n",
    "# data = quandl.get_table('WIKI/PRICES', ticker = symbols, \n",
    "#                         qopts = { 'columns': ['ticker', 'date', 'adj_close'] }, \n",
    "#                         date = { 'gte': '2000-01-01', 'lte': '2018-05-01' }, \n",
    "#                         paginate=True)\n",
    "# pickle.dump(data, open( \"dowjones_data.pkl\", \"wb\" ))\n",
    "# \n",
    "# data.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since it is faster to load locally saved data, save the downloaded data into the file.\n",
    "\n",
    "### Load saved data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "wget https://www.quandl.com/api/v3/databases/WIKI/codes?api_key=DjK3WD8BVCXtBo6QfW-5\n",
    "unzip -o codes\\?api_key\\=DjK3WD8BVCXtBo6QfW-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "symbols = pd.read_csv('WIKI-datasets-codes.csv', header=None)\n",
    "symbols.columns = ['code', 'description']\n",
    "symbols['ticker'] = symbols.code.str.split('/', expand=True)[1]\n",
    "symbols['description'] = symbols.description.str.split('\\) Prices', expand=True)[1]\n",
    "symbols.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pickle.load(open(\"dowjones_data.pkl\", \"rb\")).set_index('date')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stock data can have irregularities such as missing data due stocks being added and removed from the index. Some examples are\n",
    "* Alcoa Corp. (AA) was removed in 2013\n",
    "* Apple (AAPL) was added in 2015\n",
    "* E.I. du Pont de Nemours & Company (DD) was removed and replaced with Dow du Pont (DWDP) as a continuation in 2017\n",
    "\n",
    "For simplicity the stocks we will use are based on the most recent DJIA constituent companies until 2017-08-31 so all 30 stocks can be included. Note that this is not the same as the Dow Jones Index throughout the time period since some symbols have entered the DJIA later than year 2000: e.g., Verizon was added in 2004.\n",
    "\n",
    "We have set the index to be 'date' column. For the next line, unset datetime index in order to operate on the column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_range = data.reset_index().groupby('ticker')['date'].agg([min, max])\n",
    "data_range"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Subsetting date ranges with datetime index\n",
    "Setting `date` column as the index is useful for time series data. The row index can be used to subset the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data['2000-01-03':'2017-08-31']\n",
    "data = data.drop(['DWDP'])\n",
    "\n",
    "## new data_range\n",
    "data_range = data.reset_index().groupby('ticker')['date'].agg([min, max])\n",
    "data_range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Long to wide-format data\n",
    "\n",
    "Since we will compute the covariance matrix, turn the data into long format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datawide = data.reset_index().pivot(index='date',columns='ticker',values='adj_close')\n",
    "datawide"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Log returns from stock prices\n",
    "\n",
    "Stock market data are measured in terms of price per share. We need to compute the returns from the prices. \n",
    "Given the prices $P_t$ and $P_{t-1}$ where $t$ indicates time, the return at time $t$ is defined as\n",
    "$$ R_t = \\frac{P_t - P_{t-1}}{P_{t-1}} = \\frac{P_t}{P_{t-1}} - 1 $$\n",
    "\n",
    "It is well known that linear approximation of $\\log(1+x)\\approx x$ when $x$ is small. Since daily returns of stocks are small, we approximate that\n",
    "$$ r_t = \\log(1 + R_t) = \\log\\left(\\frac{P_t}{P_{t-1}}\\right) = \\log(P_t) - \\log(P_{t-1})$$\n",
    "So, in order to compute the log-returns, compute the difference of log prices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logret = np.log(datawide).diff()\n",
    "logret.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the first time period is NaN since there is no data corresponding to $-1$.\n",
    "\n",
    "Note that $100\\cdot r_t$% represent daily percentage returns.\n",
    "\n",
    "### Estimate expected returns\n",
    "\n",
    "Estimate the daily expected returns by computing the means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = logret[1:].mean()\n",
    "mu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate covariance matrix (volatility structure)\n",
    "\n",
    "Estimate the covarince matrix of returns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = logret.cov()\n",
    "sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize volatility structure\n",
    "\n",
    "We can visualize the relationship between the stocks by a heatmap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = plt.figure(figsize=(18, 13))\n",
    "\n",
    "ax = p.add_subplot(1,2,1)\n",
    "plt.imshow(logret.cov())\n",
    "plt.colorbar(orientation='horizontal')\n",
    "ax.set_xticks(np.arange(0, 29, 1))\n",
    "ax.set_yticks(np.arange(0, 29, 1))\n",
    "ax.set_xticklabels(sigma.columns.values, rotation=90)\n",
    "ax.set_yticklabels(sigma.columns.values)\n",
    "\n",
    "ax = p.add_subplot(1,2,2)\n",
    "plt.imshow(logret.corr())\n",
    "plt.colorbar(orientation='horizontal')\n",
    "ax.set_xticks(np.arange(0, 29, 1))\n",
    "ax.set_yticks(np.arange(0, 29, 1))\n",
    "ax.set_xticklabels(sigma.columns.values, rotation=90)\n",
    "ax.set_yticklabels(sigma.columns.values)\n",
    "plt.show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cvx\n",
    "\n",
    "s, _ = sigma.shape\n",
    "\n",
    "w = cvx.Variable(s)\n",
    "risk = cvx.quad_form(w, sigma.as_matrix())\n",
    "prob = cvx.Problem(cvx.Minimize(risk), \n",
    "               [cvx.sum_entries(w) == 1])\n",
    "prob.solve()\n",
    "w.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examples for portfolio optimization: http://nbviewer.jupyter.org/github/cvxgrp/cvx_short_course/blob/master/applications/portfolio_optimization.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though a much more difficult problem to solve, we can use `cvxpy` package to solve optimization problems easily. Let's try to solve the more difficult problem with a target expected return:\n",
    "\n",
    "$$ \\min_{x\\in\\mathbb{R}^s}\\ \\ x^\\intercal \\Sigma x\\\\\n",
    "\\text{subject to }\\mu^\\intercal x\\geq \\mu^* \\text{, and } \\mathbf{1}^\\intercal x = 1,$$\n",
    "\n",
    "Note that in all of the above problems, the measure of $\\mu^*$, expected returns, and covariance matrix have to be provided. In practice, these quantities are estimated from data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cvx\n",
    "\n",
    "s,_ = sigma.shape\n",
    "\n",
    "w = cvx.Variable(s)\n",
    "risk = cvx.quad_form(w, sigma.as_matrix())\n",
    "prob = cvx.Problem(cvx.Minimize(risk), \n",
    "               [\n",
    "                   cvx.sum_entries(w) == 1,\n",
    "                   mu.as_matrix()*w >= 0.001\n",
    "               ])\n",
    "prob.solve()\n",
    "w.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify that the solution satisfies the constraints: i.e. $\\mu^\\intercal w=0.001$ and $\\mathbf{1}^\\intercal w=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.value.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.dot(mu.as_matrix(),w.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earned = np.dot(logret.fillna(method='bfill').as_matrix(), w.value)\n",
    "# f = plt.figure(figsize=(15,15))\n",
    "# plt.plot((1+earned[1:]).cumprod())\n",
    "# plt.show(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earnedcp = (1+earned[1:]).cumprod()\n",
    "ecp = pd.DataFrame(earnedcp.T, columns=['Portfolio Value'])\n",
    "ecp.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reality check\n",
    "\n",
    "Obviously, this cannot be realistic. What aspects were unrealistic?"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}