Quiz/m8/optimization_with_tcosts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization with Transaction costs\n",
    "\n",
    "In this lesson, we’ll show you how to incorporate transaction costs into portfolio optimization.  This will give your backtest a more realistic measure of your alpha’s performance.  In addition, we’ll show you some additional ways to design your optimization with efficiency in mind.  This is really helpful when  backtesting, because having reasonably shorter runtimes allows you to test and iterate on your alphas more quickly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "!{sys.executable} -m pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "import gzip\n",
    "import bz2\n",
    "from statsmodels.formula.api import ols\n",
    "from scipy.stats import gaussian_kde\n",
    "import scipy\n",
    "import scipy.sparse\n",
    "import patsy\n",
    "from statistics import median\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Barra data\n",
    "\n",
    "We’ll be using factor data that is generated by Barra.  This will be good practice because Barra data is used throughout the industry.  \n",
    "\n",
    "Note that we've pre-processed the raw barra data files and stored the data into pickle files. The alternative would be to load the original data, and perform the parsing each time.  Since parsing and pre-processing takes time, we recommend doing the pre-processing once and saving the pre-processed data for later use in your backtest.\n",
    "\n",
    "Choose the number of years to use for the backtest.  The data is available for years 2003 to 2008 inclusive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "barra_dir = '../../data/project_8_barra/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls {barra_dir}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = {}\n",
    "for year in [2003]:\n",
    "    fil = barra_dir + \"pandas-frames.\" + str(year) + \".pickle\"\n",
    "    data.update(pickle.load( open( fil, \"rb\" ) ))\n",
    "    \n",
    "covariance = {}\n",
    "for year in [2003]:\n",
    "    fil = barra_dir + \"covariance.\" + str(year) + \".pickle\"\n",
    "    covariance.update(pickle.load( open(fil, \"rb\" ) ))\n",
    "    \n",
    "daily_return = {}\n",
    "for year in [2003, 2004]:\n",
    "    fil = barra_dir + \"price.\" + str(year) + \".pickle\"\n",
    "    daily_return.update(pickle.load( open(fil, \"rb\" ) ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the frames variale is a dictionary, where the keys are strings representing each business day."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View the Barra data\n",
    "\n",
    "We'll take a look at the value stored for a single day (it's a data frame).\n",
    "\n",
    "As a general reminder of best practices, remember to check what unit of measure your data is in. In some cases, the unit of measure isn’t available in the documentation, so you’ll want to inspect the data to see what makes sense.\n",
    "\n",
    "For instance, there are volatility fields that are large enough that we can assume they are in percentage units, as opposed to decimal values.  In other cases, when we look at daily volume, we may not have documentation about whether the units are in number of shares or in dollars. One way to find this out is to spot check a single stock on a single day, and cross-reference with another source, such as Bloomberg or Yahoo Finance.\n",
    "Remember to inspect the data before you use it, as it will help you derive more meaningful results in your portfolio optimization, and in your backtest.\n",
    "\n",
    "Remember to inspect the data before you use it, as it will help you derive more meaningful results in your portfolio optimization, and in your backtest.\n",
    "\n",
    "In the exercise, we'll re-scale the data before using it, and there will be comments to point out when we re-scale the data.  So don't worry about adjusting anything here, just take a look to get familiar with the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['20030102'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['20030102'].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Factors\n",
    "\n",
    "Note that the data fields that start with the prefix U-S-F-A-S-T are factor exposures, one column for each factor.  We will use some of these as alpha factors, and the rest as risk factors.  The reason this makes sense is that, for the time periods in which we’re back-testing, some of these factors were able to produce better than average returns.  Barra works with its clients (funds) and gathers information about alphas that worked in the past.  These were calculated on historical data to produce the factor exposure data found in the Barra data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![alt text](./images_optimization/barra_field_factor_exposure.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Factors\n",
    "\n",
    "Here's a partial list of the barra factors in our dataset and their definitions.  These are collected from documentation by Barra.  There are style factors and industry factors.  The industry factors will be used as risk factors.  You can consider using the style factors as alpha factors.  Any factors not used as alpha factors can be included in the risk factors category.\n",
    "\n",
    "#### Style factors\n",
    "* beta: Describes market risk that cannot be explained by the Country factor. The Beta factor is typically the most important style factor. We calculate Beta by time-series regression of stock excess returns against the market return.\n",
    "* 1 day reversal\n",
    "* dividend yield: Describes differences in stock returns attributable to stock's historical and predicted dividend-to-price ratios.\n",
    "* downside risk (maximum drawdown)\n",
    "* earnings quality:  Describes stock return differences due to the accrual components of earnings.\n",
    "* earnings yield: Describes return differences based on a company’s earnings relative to its price. Earnings Yield is considered by many investors to be a strong value signal. The most important descriptor in this factor is the analyst-predicted 12-month earnings-to-price ratio.\n",
    "* growth: Differentiates stocks based on their prospects for sales or earnings growth. The most important descriptor in this factor is the analyst predicted long-term earnings growth. Other descriptors include sales and earnings growth over the previous five years.\n",
    "* leverage: Describes return differences between high and low-leverage stocks. The descriptors within this style factor include market leverage, book leverage, and debt-to-assets ratio.\n",
    "* liquidity: Describes return differences due to relative trading activity. The descriptors for this factor are based on the fraction of total shares outstanding that trade over a recent window.\n",
    "* long-term reversal: Describes common variation in returns related to a long-term (five years ex. recent thirteen months) stock price behavior.\n",
    "* management quality\n",
    "* Mid capitalization: Describes non-linearity in the payoff to the Size factor across the market-cap spectrum. This factor is based on a single raw descriptor: the cube of the Size exposure. However, because this raw descriptor is highly collinear with the Size factor, it is orthogonalized with respect to Size. This procedure does not affect the fit of the model, but does mitigate the confounding effects of collinearity, while preserving an intuitive meaning for the Size factor. As described by Menchero (2010), the Mid Capitalization factor roughly captures the risk of a “barbell portfolio” that is long mid-cap stocks and short small-cap and large-cap stocks.\n",
    "* Momentum – Differentiates stocks based on their performance over the trailing 12 months. When computing Momentum exposures, we exclude the most recent returns in order to avoid the effects of short-term reversal. The Momentum factor is often the second strongest factor in the model, although sometimes it may surpass Beta in importance.\n",
    "* Profitability – Combines profitability measures that characterize efficiency of a firm's operations and total activities.\n",
    "* Residual Volatility – Measures the idiosyncratic volatility anomaly. It has three descriptors: (a) the volatility of daily excess returns, (b) the volatility of daily residual returns, and (c) the cumulative range of the stock over the last 12 months. Since these descriptors tend to be highly collinear with the Beta factor, the Residual Volatility factor is orthogonalized with respect to the Beta and Size factors.\n",
    "* seasonality\n",
    "* sentiment\n",
    "* Size – Represents a strong source of equity return covariance, and captures return differences between large-cap and small-cap stocks. We measure Size by the log of market capitalization.\n",
    "* Short term reversal\n",
    "* Value\n",
    "* Prospect -- is a function of skewness and maximum drawdown.\n",
    "* Management Quality -- is a function of the following: \n",
    "    * Asset Growth: Annual reported company assets are regressed against time over the past five fiscal years. The slope coefficient is then divided by the average annual assets to obtain the asset growth.\n",
    "    * Issuance Growth Annual reported company number of shares outstanding regressed against time over the past five fiscal years. The slope coefficient is then divided by the average annual number of shares outstanding.\n",
    "    * Capital Expenditure Growth: Annual reported company capital expenditures are regressed against time over the past five fiscal years. The slope coefficient is then divided by the average annual capital expenditures to obtain the capital expenditures growth.\n",
    "    * Capital Expenditure: The most recent capital expenditures are scaled by the average of capital expenditures over the last five fiscal years.\n",
    "\n",
    "\n",
    "\n",
    "#### Industry Factors\n",
    "* aerospace and defense\n",
    "* airlines\n",
    "* aluminum and steel\n",
    "* apparel\n",
    "* Automotive\n",
    "* banks\n",
    "* beta (market)\n",
    "* beverage and tobacco\n",
    "* biotech & life science\n",
    "* building products\n",
    "* chemicals\n",
    "* construction & engineering\n",
    "* construction & machinery\n",
    "* construction materials\n",
    "* commercial equipment\n",
    "* computer & electronics\n",
    "* commercial services\n",
    "* industrial conglomerates\n",
    "* containers (forest, paper, & packaging)\n",
    "* distributors\n",
    "* diversified financials\n",
    "* electrical equipment\n",
    "* electrical utility\n",
    "* food & household products & personal\n",
    "* food & staples retailing\n",
    "* gas & multi-utilities\n",
    "* healthcare equipment and services\n",
    "* health services\n",
    "* home building\n",
    "* household durables\n",
    "* industry machinery\n",
    "* non-life insurance\n",
    "* leisure products\n",
    "* leisure services\n",
    "* life insurance\n",
    "* managed healthcare\n",
    "* multi-utilities\n",
    "* oil & gas conversion\n",
    "* oil & gas drilling\n",
    "* oil & gas equipment\n",
    "* oil and gas export\n",
    "* paper\n",
    "* pharmaceuticals\n",
    "* precious metals\n",
    "* personal products\n",
    "* real estate\n",
    "* restaurants\n",
    "* road & rail\n",
    "* semiconductors\n",
    "* semiconductors equipment\n",
    "* software\n",
    "* telecommunications\n",
    "* transportation\n",
    "* wireless\n",
    "* SPTY\\* and SPLTY\\* are various industries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['20030102'].columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## covariance of factors\n",
    "\n",
    "Let's look at the covariance of the factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covariance.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "View the data for a single day.  Notice that the factors are listed in two columns, followed by the covariance between them.  We'll use this data later to create a factor covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covariance['20030102'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Daily returns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "daily_return.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "daily_return['20030102'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add date for returns\n",
    "\n",
    "We'll be dealing with two different dates; to help us keep track, let's add an additional column in the daily_return dataframes that stores the date of the returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_date = '20030102'\n",
    "tmp = daily_return[tmp_date]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_n_rows = tmp.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.Series([tmp_date]*tmp_n_rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp['DlyReturnDate'] = pd.Series([tmp_date]*tmp_n_rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: add daily return date to each dataframe in daily_return dictionary\n",
    "\n",
    "Name the column `DlyReturnDate`.\n",
    "**Hint**: create a list containing copies of the date, then create a pandas series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for DlyReturnDate, df in daily_return.items():\n",
    "    # TODO\n",
    "    n_rows = # ...\n",
    "    df['DlyReturnDate'] = # ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check results\n",
    "\n",
    "daily_return['20030102'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adjust dates to account for trade execution\n",
    "\n",
    "The data stored in `data` and `covariance` are used to choose the optimal portfolio, whereas the data in `daily_return` represents the the returns that the optimized portfolio would realize, but only after we've received the data, then chosen the optimal holdings, and allowed a day to trade into the optimal holdings.  In other words, if we use the data from `data` and `covariance` that is collected at the end of Monday, we'll use portfolio optimization to choose the optimal holdings based on this data, perhaps after hours on Monday.  Then on Tuesday, we'll have a day to execute trades to adjust the portfolio into the optimized positions.  Then on Wednesday, we'll realize the returns using those optimal holdings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of what we want\n",
    "data_date_l = sorted(data.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "return_date_l = sorted(daily_return.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(data_date_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(return_date_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "return_date_l_shifted = return_date_l[2:len(data) + 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(return_date_l_shifted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data date\n",
    "data_date_l[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# returns date\n",
    "return_date_l_shifted[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = data['20030102'].merge(daily_return['20030102'], on=\"Barrid\")\n",
    "tmp.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge data and daily returns into single dataframe\n",
    "\n",
    "Use a loop to merge the `data` and `daily_return` tables on the `barrid` column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frames ={}\n",
    "# TODO\n",
    "dlyreturn_n_days_delay = # ...\n",
    "\n",
    "# TODO\n",
    "date_shifts = zip(\n",
    "        \"\"\"<your code here>\"\"\",\n",
    "        \"\"\"<your code here>\"\"\")\n",
    "\n",
    "# TODO\n",
    "for data_date, price_date in date_shifts:\n",
    "    frames[price_date] = # ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's work with a single day's data. Later, we'll put this into a loop\n",
    "\n",
    "\n",
    "Notice how the keys are now dates of the returns.  So the earliest date in \"frames\" dictionary is two business days after the earliest date in \"data\" dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frames.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = frames['20030106']\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Filter the stocks so that the estimation universe has stocks with at least 1 billion in market cap.  As an aside, it doesn't make much of a difference whether we choose a \">\" or \">=\", since the threshold we choose is just meant to get a set of relatively liquid assets.\n",
    "\n",
    "**Hint**: use `.copy(deep=True)` to make an independent copy of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "estu = # ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estu.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For all the columns in the dataframe, the ones with the prefix \"USFAST\" are factors.  We'll use a helper function to get the list of factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def factors_from_names(n):\n",
    "    return(list(filter(lambda x: \"USFASTD_\" in x, n)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors = factors_from_names(list(df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## factors exposures and factor returns\n",
    "\n",
    "Recall that a factor's factor return times its factor exposure gives the part of a stock's return that is explained by that factor.\n",
    "\n",
    "The Barra data contains the factor exposure of each factor.  We'll use regression to estimate the factor returns of each factor, on each day.  The observations will be the cross section of stock factor exposures, as well as the stock returns that are realized two trading days later.  Recall from an earlier lesson that this is a cross-sectional regression, because it's a cross section of stocks, for a single time period.\n",
    "\n",
    "$r_{i,t} = \\sum_{j=1}^{k} (\\beta_{i,j,t-2} \\times f_{j,t})$  \n",
    "where $i=1...N$ (N assets),   \n",
    "and $j=1...k$ (k factors).\n",
    "\n",
    "In the regression, the factor exposure, $\\beta_{i,j,t-2}$ is the independent variable, $r_{i,t}$ is the dependent variable, and the factor return $f_{j,t}$ is the coefficient that we'll estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating factor returns\n",
    "\n",
    "We'll estimate the factor returns $f_{j,t}$ of our chosen alpha factors, using the daily returns of the stocks $r_{i,t}$, where $i=1...N$ and the factor exposure $\\beta_{i,j,t-2}$ of each stock to each factor.  \n",
    "\n",
    "Note that we'll use a universe of stocks where the companies have a market capitalization of at least 1 billion.  The factor returns estimated would be slightly different depending on which stock universe is chosen, but choosing a market cap of 1 billion or more provides a reasonable estimate of what you'd expect to be tradable.  The estimated factor returns would be fairly close to what you'd find if you used the Russell 3000 index as the stock universe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## formula\n",
    "\n",
    "We'll use a helper function that creates a string that defines which are the independent and dependent variables for a model to use.  This string is called a \"formula.\"  We'll use this in the regression, and later again when we work with matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_formula(factors, Y):\n",
    "    L = [\"0\"]\n",
    "    L.extend(factors)\n",
    "    return Y + \" ~ \" + \" + \".join(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "form = get_formula(all_factors, \"DlyReturn\")    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the formula is saying `DlyReturn` is the dependent variable, whereas the `USFAST...` columns are the independent variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Run an ordinary least squares regression\n",
    "\n",
    "[ols documentation](https://www.statsmodels.org/dev/example_formulas.html)\n",
    "\n",
    "Here's an example of the syntax.\n",
    "```\n",
    "ols(formula='y ~ x1 + x2 + x3', data=dataframe)\n",
    "```\n",
    "\n",
    "Note that you're free to choose other regression models, such as ridge, lasso, or elastic net.  These may give you slightly different estimations of factor returns, but shouldn't be too different from each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: create the ols model\n",
    "model = # ...\n",
    "\n",
    "# TODO: fit the model\n",
    "results = # ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the factor data that we're using as the independent variables are the factor exposures, the coefficients estimated by the regression are the estimated factor returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: winsorize daily returns before calculating factor returns\n",
    "\n",
    "We're going to use regression to estimate the factor returns of all the factors.  To avoid using extreme values in the regression, we'll winsorize, or \"clip\" the returns.  We can check the data distribution using a density plot.\n",
    "\n",
    "Note that [numpy.where](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.where.html) has the form \n",
    "```\n",
    "numpy.where(<condition>, <value if true>, <value if false>)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wins(x,wins_lower,wins_upper):\n",
    "    #TODO\n",
    "    clipped_upper = # ...\n",
    "    clipped_both = # ...\n",
    "    return clipped_both"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A density plot will help us visually check the effect of winsorizing returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def density_plot(data): \n",
    "    density = gaussian_kde(data)\n",
    "    xs = np.linspace(np.min(data),np.max(data),200)\n",
    "    density.covariance_factor = lambda : .25\n",
    "    density._compute_covariance()\n",
    "    plt.plot(xs,density(xs))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# distribution without winsorizing\n",
    "test = frames['20040102']\n",
    "density_plot(test['DlyReturn'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# distribution after winsorizing\n",
    "test['DlyReturn_wins'] = wins(test['DlyReturn'],lower_bound=-0.1,upper_bound=0.1)\n",
    "density_plot(test['DlyReturn_wins'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Put the factor returns estimation into a function, so that this can be re-used for each day's data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_factor_returns(df, wins_lower=-.25, wins_upper=0.25): \n",
    "    ## TODO: build estimation universe based on filters \n",
    "    estu = # ...\n",
    "  \n",
    "    ## TODO: winsorize returns for fitting \n",
    "    estu['DlyReturn'] = # ...\n",
    "  \n",
    "    ## get a list of all the factors\n",
    "    all_factors = # ...\n",
    "    \n",
    "    ## define a 'formula' for the regression\n",
    "    form = # ...\n",
    "    \n",
    "    ## create the OLS model, passing in the formula and the estimation universe dataframe\n",
    "    model = # ...\n",
    "    \n",
    "    ## return the estimated coefficients\n",
    "    results = # ...\n",
    "    return(results.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choose alpha factors\n",
    "\n",
    "We'll choose the 1 day reversal, earnings yield, value, and sentiment factors as alpha factors.  We'll calculate the factor returns of these alpha factors to see how they performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_factors = [\"USFASTD_1DREVRSL\", \"USFASTD_EARNYILD\", \"USFASTD_VALUE\", \"USFASTD_SENTMT\"]\n",
    "print(alpha_factors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: estimate factor returns of alpha factors\n",
    "\n",
    "Loop through each day, and estimate the factors returns of each factor, that date, in the `frames` dictionary.  This may take a minute or more to run per year of data used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "facret = {}\n",
    "for date in frames:\n",
    "    # TODO: store factor returns as key-value pairs in a dictionary\n",
    "    facret[date] = # ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(facret['20040102'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "facret['20040102'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## put the factor returns into a dataframe\n",
    "\n",
    "The pandas series are stored inside a dictionary.  We'll put the factor returns into a dataframe where the rows are the dates and the columns are the factor returns (one column for each factor).\n",
    "\n",
    "First, let's get a list of dates, as Timestamp objects.  We'll use [pandas.to_datetime](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.to_datetime.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example of how to convert the keys of the dataframe into Timestamp objects\n",
    "pd.to_datetime('20040102', format='%Y%m%d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Store the timestamp objects in a list (can use a list comprehension, or for loop)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "dates_unsorted = # ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sort the dates in ascending order\n",
    "my_dates = sorted(dates_unsorted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll make an empty dataframe with the dates set as the row index.\n",
    "facret_df = pd.DataFrame(index = my_dates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "facret_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rows are the dates.  The columns will be the factor returns. \n",
    "\n",
    "To convert from Timestamp objects back into a string, we can use [Timestamp.strftime('%Y%m%d')](https://www.programiz.com/python-programming/datetime/strftime). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## example usage of Timestamp.strftime('%Y%m%d')\n",
    "my_dates[0].strftime('%Y%m%d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "For each date, and for each factor return, get the value from the dictionary and put it into the dataframe.\n",
    "\n",
    "We can use [pandas.DataFrame.at¶](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.at.html), \n",
    "\n",
    "```\n",
    "DataFrame.at[<index_value>,<column_name>] = <some_value>\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: for each date (rows), and for each factor (columns),\n",
    "# store factor return in the dataframe\n",
    "for dt in my_dates: \n",
    "    for alp in alpha_factors:\n",
    "        facret_df.at[dt, alp] = # ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Portfolio optimization for a single period\n",
    "\n",
    "When we get to the project, we'll want to define the portfolio optimization within a function.  For now, let's walk through the steps we'll take in separate cells, so that we can see what's going on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimization will want to know about the prior trading day's portfolio holdings, also called holdings.  The previous day's holdings will be used to estimate the size of the trades due to position changes, which in turn helps us estimate transaction costs.  We'll start with an initial holding of zero for a single stock.  The reason we'll use a single stock is that the estimation universe chosen on each day will include all stocks that have holdings on the previous day.  So we want to keep this list small when we first start out, else we'll keep many stocks that may fall below the 1 billion market cap threshold, just because they were chosen in the initialization of the backtest.\n",
    "\n",
    "We'll want to choose a stock that is likely to satisfy the 1 billion market cap threshold on any day.  So let's choose the stock with the largest market cap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we're going to set a single barra id to have a zero portfolio holding, \n",
    "# so let's pick any barra id of the stock with the largest issuer market cap\n",
    "estu.sort_values('IssuerMarketCap',ascending=False)[['Barrid','IssuerMarketCap']].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: Intialize previous holdings dataframe\n",
    "\n",
    "Create a new dataframe and initialize it with a dictionary, where the key is \"Barrid\" followed by a value that is a pandas.Series containing the barra id of the largest market cap in the stock universe.  \n",
    "\n",
    "Also set another key value pair to \"x.opt.previous\" and the value is set to a pandas.Series with the value 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "# create a dataframe of previous holdings, \n",
    "# initializing a single stock (barra id) to zero portfolio holding\n",
    "previous_holdings = pd.DataFrame(data = {\"Barrid\" : \"\"\"<your code here>\"\"\", \n",
    "                                         \"x.opt.previous\" : \"\"\"<your code here>\"\"\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "previous_holdings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get a single day's data to be used for the portfolio optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = my_dates[0]\n",
    "date = dt.strftime('%Y%m%d')\n",
    "print(date)\n",
    "df = frames[date]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add the previous holdings column to the dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## merge previous portfolio holdings \n",
    "df = df.merge(previous_holdings, how = 'left', on = 'Barrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clean missing and zero values.\n",
    "\n",
    "First replace missing values with zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def na2z(x):\n",
    "    return(np.nan_to_num(x))\n",
    "\n",
    "def names_numeric_columns(df): \n",
    "    return(df.select_dtypes(include=[np.number]).columns.tolist())\n",
    "\n",
    "def clean_nas(df):\n",
    "    for x in names_numeric_columns(df): \n",
    "        df[x] = na2z(df[x])\n",
    "    return(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = clean_nas(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: Clean specific risk\n",
    "Barra calculates specific risk for each asset.  If the value in the data is zero, this may be due to missing data rather than the specific risk actually being zero.  So we'll set zero values to the median, to make sure our model is more realistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: if SpecRisk is zero, set it to median\n",
    "df.loc[df['SpecRisk'] == \"\"\"<enter a number here>\"\"\"]['SpecRisk'] = median(\"\"\"<your code here>\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## universe\n",
    "\n",
    "We'll look at stocks that are 1 billion in market cap or greater.  An important point here is that we'll need to account for stocks that are already in our portfolio, even if the market cap of the stock is no longer 1 billion on the current day."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Quiz: think about what would happen if we had an existing position in a stock, then the market cap fell below the threshold and the stock was excluded from the stock universe.  What would happen to the position on that stock?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer\n",
    "The stock would not be included in the optimization, which means it would be given a zero position.  So this effectively says to sell all holdings in the asset once it falls below the market cap threshold.  That's not what we want to do."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modify the code to account for the prior day's positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: modify the given code to include the prior day's assets\n",
    "universe = df.loc[(df['IssuerMarketCap'] >= 1e9)].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: Nothing here should be allowed to look at returns when forming the portfolio.\n",
    "Make this impossible by removing the Daily returns data from the dataframe.  Drop the DlyReturn field from the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: drop DlyReturn column\n",
    "universe = # ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## this will extract all of the factors, including the alphas \n",
    "# list(universe) gets a list of the column names of the dataframe\n",
    "all_factors = factors_from_names(list(universe))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alpha factors\n",
    "\n",
    "Just a reminder that we chose four of these factors that represent previously effective alpha factors.  Since these factors became well known over time, they were added to the Barra data set.  For the time frame that we're running the back-test, these were effective alpha factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_factors #4 alpha factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: risk factors\n",
    "\n",
    "The risk factors we'll use are all the factors that are not alpha factors.  Complete the setdiff function so that it takes a superset, a subset, and returns the difference as a set.\n",
    "\n",
    "diff= SuperSet \\ Subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setdiff(superset, subset): \n",
    "    # TODO: create a set from subset\n",
    "    s = # ...\n",
    "\n",
    "    # TODO: get the difference\n",
    "    diffset = # ...\n",
    "    return(diffset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "risk_factors = setdiff(all_factors, alpha_factors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 77 risk factors\n",
    "len(risk_factors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save initial holdings in a variable for easier access.  We'll later use it in matrix multiplications, so let's convert this to a numpy array.  We'll also use another variable to represent the current holdings, which are to be run through the optimizer.  We'll set this to be a copy of the previous holdings.  Later the optimizer will continually update this to optimize the objective function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## initial holdings (before optimization)\n",
    "# optimal holding from prior day\n",
    "h0 = np.asarray( universe['x.opt.previous'] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = h0.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix of Risk Factor Exposures $\\textbf{B}$\n",
    "\n",
    "The dataframe contains several columns that we'll use as risk factors exposures.  Extract these and put them into a matrix.\n",
    "\n",
    "The data, such as industry category, are already one-hot encoded, but if this were not the case, then using `patsy.dmatrices` would help, as this function extracts categories and performs the one-hot encoding.  We'll practice using this package, as you may find it useful with future data sets.  You could also store the factors in a dataframe if you prefer to avoid using patsy.dmatrices.\n",
    "\n",
    "#### How to use patsy.dmatrices\n",
    "\n",
    "patsy.dmatrices takes in a formula and the dataframe.  The formula tells the function which columns to take.  The formula will look something like this:  \n",
    "`SpecRisk ~ 0 + USFASTD_AERODEF + USFASTD_AIRLINES + ...`  \n",
    "where the variable to the left of the ~ is the \"dependent variable\" and the others to the right are the independent variables (as if we were preparing data to be fit to a model).\n",
    "\n",
    "This just means that the pasty.dmatrices function will return two matrix variables, one that contains the single column for the dependent variable `outcome`, and the independent variable columns are stored in a matrix `predictors`.\n",
    "\n",
    "The `predictors` matrix will contain the matrix of risk factors, which is what we want.  We don't actually need the `outcome` matrix; it's just created because that's the way patsy.dmatrices works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that we chose \"SpecRisk\" simply because it's not one of the USFAST factors.\n",
    "# it will be discarded in the next step.\n",
    "formula = get_formula(risk_factors, \"SpecRisk\")\n",
    "formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the factors will be in the second returned variable (predictors)\n",
    "# the outcome variable contains the SpecRisk data, which we don't actually need here\n",
    "outcome, predictors = patsy.dmatrices(formula,universe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`predictors` contains the factor exposures of each asset to each factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictors.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Factor exposure matrix $\\textbf{B}$\n",
    "\n",
    "\n",
    "Remember, the factor exposure matrix has the exposure of each asset to each factor.  Thee number of rows is number of assets, and number of columns is the number of factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NROW(x):\n",
    "    return(np.shape(x)[0])\n",
    "\n",
    "def NCOL(x):\n",
    "    return(np.shape(x)[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Set the factor exposure matrix and its transpose, using one of the outputs from calling patsy.dmatrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: risk exposure matrix: \n",
    "B = # ...\n",
    "BT = B.transpose()\n",
    "\n",
    "k = NCOL(B) #number of factors (77)\n",
    "n = NROW(B) #number of assets (2000+)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Factor covariance matrix $\\textbf{F}$\n",
    "\n",
    "We can improve on the factor covariance matrix by reducing noise and also increasing computational efficiency.\n",
    "\n",
    "If we have, 70 risk factors in our risk model, then the covariance matrix of factors is a 70 by 70 square matrix.  The diagonal contains the variances of each factor, while the off-diagonals contain the pairwise covariances of two different risk factors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, it’s good to have a healthy suspicion of correlations and covariances, and to ask if correlation data adds information or just more noise.  One way to be conservative about the information in a covariance matrix is to shrink the covariances, or even reduce them to zero.  In other words, we could keep just the variances along the diagonal, and set the covariances in the off-diagonals to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case where we’re using the covariance matrix in a risk factor model, there’s also some additional intuition for why we can try using just the variances, and discard the covariances.  The goal of the optimizer is to reduce the portfolio’s exposure to these risk factors.  So if the optimizer reduces the portfolio’s exposure to risk factor “one”, and also reduces its exposure to risk factor “two”, then it’s less important to know exactly how factor one varies with factor two.\n",
    "\n",
    "You may wonder what are the benefits of throwing away the information about the covariances.  In addition to making your model more conservative, and limiting possible noise in your data, a diagonal matrix also makes matrix operations more efficient.  This theme of computational efficiency is one that you’ll come across in many use cases, including backtesting.  Backtesting is a computationally and time-intensive process, so the more efficient you can make it, the more quickly you can test your alphas, and iterate to make improvements.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Factor covariance matrix $\\textbf{F}$\n",
    "\n",
    "You can try getting all covariances into the matrix.  Notice that we'll run into some issues where the covariance data doesn't exist.\n",
    "\n",
    "One important point to remember is that we need to order the factors in the covariance matrix F so that they match up with the order of the factors in the factor exposures matrix B.\n",
    "\n",
    "Note that covariance data is in percentage units squared, so to use decimals, so we'll rescale it to convert it to decimal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## With all covariances\n",
    "def colnames(X):\n",
    "    if(type(X) == patsy.design_info.DesignMatrix): \n",
    "        return(X.design_info.column_names)\n",
    "    if(type(X) == pandas.core.frame.DataFrame): \n",
    "        return(X.columns.tolist())\n",
    "    return(None)\n",
    "\n",
    "## extract a diagonal element from the factor covariance matrix \n",
    "def get_cov_version1(cv, factor1, factor2): \n",
    "    try:\n",
    "        return(cv.loc[(cv.Factor1==factor1) & (cv.Factor2==factor2),\"VarCovar\"].iloc[0])\n",
    "    except:\n",
    "        print(f\"didn't find covariance for: factor 1: {factor1} factor2: {factor2}\")\n",
    "        return 0\n",
    "\n",
    "def diagonal_factor_cov_version1(date, B):\n",
    "    \"\"\"\n",
    "    Notice that we'll use the order of column names of the factor exposure matrix\n",
    "    to set the order of factors in the factor covariance matrix\n",
    "    \"\"\"\n",
    "    cv = covariance[date]\n",
    "    k = NCOL(B)\n",
    "    Fm = np.zeros([k,k])\n",
    "    for i in range(0,k):\n",
    "        for j in range(0,k): \n",
    "            fac1 = colnames(B)[i]\n",
    "            fac2 = colnames(B)[j]\n",
    "            # Convert from percentage units squared to decimal\n",
    "            Fm[i,j] = (0.01**2) * get_cov_version1(cv, fac1, fac2)\n",
    "    return(Fm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an example where the two factors don't have covariance data for the date selected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv = covariance['20031211']\n",
    "cv.loc[(cv.Factor1=='USFASTD_AERODEF') & (cv.Factor2=='USFASTD_ALUMSTEL')]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see where all the factor covariances aren't found in the data.\n",
    "\n",
    "## Which date?\n",
    "\n",
    "Recall that there's a DataDate column and DlyReturnDate column in the dataframe.  We're going to use a date to access the covariance data.  Which date should we use?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Answer here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "Choose the correct date, then use the `diagonal_factor_cov_version1` to get the factor covariance matrix of that date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "date = str(int(universe[\"\"\"<your code here>\"\"\"][1]))\n",
    "print(date, end =\" \")\n",
    "F_version1 = diagonal_factor_cov_version1(date, B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: Create matrix of factor variances\n",
    "\n",
    "Just use the factor variances and set the off diagonal covariances to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def colnames(X):\n",
    "    if(type(X) == patsy.design_info.DesignMatrix): \n",
    "        return(X.design_info.column_names)\n",
    "    if(type(X) == pandas.core.frame.DataFrame): \n",
    "        return(X.columns.tolist())\n",
    "    return(None)\n",
    "\n",
    "## extract a diagonal element from the factor covariance matrix \n",
    "def get_var(cv, factor): \n",
    "    # TODO\n",
    "    return(cv.loc[(cv.Factor1==\"\"\"<your code here>\"\"\") & (\"\"\"<your code here>\"\"\"==\"\"\"<your code here>\"\"\"),\"VarCovar\"].iloc[0])\n",
    "\n",
    "\n",
    "def diagonal_factor_cov(date, B):\n",
    "    \"\"\"\n",
    "    Notice that we'll use the order of column names of the factor exposure matrix\n",
    "    to set the order of factors in the factor covariance matrix\n",
    "    \"\"\"\n",
    "    # TODO: set the variances only\n",
    "    cv = covariance[date]\n",
    "    k = NCOL(B)\n",
    "    Fm = np.zeros([k,k])\n",
    "    for j in range(0,k): \n",
    "        fac = colnames(B)[j]\n",
    "        Fm[j,j] = # ...\n",
    "    return(Fm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## factor variances \n",
    "# gets factor vars into diagonal matrix\n",
    "# takes B to know column names of B; F will be multipled by B later\n",
    "# F is square; so row and col names must match column names of B.\n",
    "F = diagonal_factor_cov(date, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the off diagonals are all set to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## alpha combination\n",
    "\n",
    "As a simple alpha combination, combine the alphas with equal weight.  In the project, you're welcome to try other ways to combine the alphas.  For example, you could calculate some metric for each factor, which indicates which factor should be given more or less weight.\n",
    "\n",
    "## Scale factor exposures \n",
    "\n",
    "Note that the terms that we're calculating for the objective function will be in dollar units.  So the expected return $-\\alpha^T h$ will be in dollar units.  The $h$ vector of portfolio holdings will be in dollar units.  The vector of alpha factor exposures $\\alpha$ will represent the percent change expected for each stock.  Based on the ranges of values in the factor exposure data, which are mostly between -5 and +5 and centered at zero, **we'll make an assumption that a factor exposure of 1 maps to 1 basis point of daily return on that stock.**\n",
    "\n",
    "So we'll convert the factor values into decimals: 1 factor exposure value $\\rightarrow \\frac{1}{10,000}$ in daily returns.  In other words, we'll rescale the alpha factors by dividing by 10,000.\n",
    "\n",
    "This is to make the term representing the portfolio's expected return $\\alpha^T h$ be scaled so that it represents dollar units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_matrix(formula, data): \n",
    "    outcome, predictors = patsy.dmatrices(formula, data)\n",
    "    return(predictors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## matrix of alpha factors\n",
    "B_alpha = model_matrix(get_formula(alpha_factors, \"SpecRisk\"), data = universe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B_alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Sum across the rows, then re-scale so that the expression $\\mathbf{\\alpha}^T \\mathbf{h}$ is in dollar units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rowSums(m): \n",
    "    # TODO\n",
    "    return(np.sum(m, axis=\"\"\"<your code here>\"\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "scale = # ...\n",
    "alpha_vec = scale * rowSums(B_alpha) #sum across rows (collapse 4 columns into one)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_vec.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Original method of calculating common risk term\n",
    "\n",
    "Recall that the common risk term looks like this:\n",
    "$\\textbf{h}^T\\textbf{BFB}^T\\textbf{h}$\n",
    "\n",
    "Where h is the vector of portfolio holdings, B is the factor exposure matrix, and F is the factor covariance matrix.\n",
    "\n",
    "We'll walk through this calculation to show how it forms an N by N matrix, which is computationally expensive, and may lead to memory overflow for large values of N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.dot( np.dot( h.T, np.matmul( np.matmul(B,F),BT) ), h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = np.matmul(B,F)\n",
    "tmp.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this makes an N by matrix (large)\n",
    "tmp = np.matmul(tmp,BT)\n",
    "tmp.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = np.matmul(h.T,tmp)\n",
    "tmp.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = np.dot(tmp,h)\n",
    "tmp.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Efficiently calculate common risk term (avoid N by N matrix)\n",
    "\n",
    "Calculate the portfolio risk that is attributable to the risk factors:\n",
    "$\\mathbf{h}^T\\mathbf{BFB}^T\\mathbf{h}$\n",
    "\n",
    "Note that this can become computationally infeasible and/or slow.  Use matrix factorization and carefully choose the order of matrix multiplications to avoid creating an N by N matrix.\n",
    "\n",
    "#### square root of a matrix.\n",
    "\n",
    "We can find a matrix $\\mathbf{B}$ that's the matrix square root of another matrix $\\mathbf{A}$, which means that if we matrix multiply $\\mathbf{BB}$, we'd get back to the original matrix $\\mathbf{A}$.\n",
    "\n",
    "Find $\\mathbf{Q}$ such that $\\mathbf{Q}^T\\mathbf{Q}$ is the same as $\\mathbf{BFB}^T$.  Let's let $\\mathbf{G}$ denote the square root of matrix $\\mathbf{F}$, so that $\\mathbf{GG} = \\mathbf{F}$.\n",
    "\n",
    "Then the expression for the covariance matrix of assets, $\\mathbf{BFB}^T$, can be written as $\\mathbf{BGGB}^T$.  \n",
    "\n",
    "Let's let $\\mathbf{Q}=\\mathbf{GB}^T$ and let $\\mathbf{Q}^T=\\mathbf{BG}$, which means we can rewrite $\\mathbf{BGGB}^T = \\mathbf{Q}^T\\mathbf{Q}$, and the common risk term is $\\mathbf{h}^T\\mathbf{Q}^T\\mathbf{Qh}$\n",
    "\n",
    "Also, note that we don't have to calculate $\\mathbf{BFB}^T$ explicitly, because the actual value we wish to calculate in the objective function will apply the holdings $\\mathbf{h}$ to the covariance matrix of assets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: matrix square root of F\n",
    "\n",
    "We'll call this square root matrix $\\mathbf{G}$\n",
    "\n",
    "Use [scipy.linalg.sqrtm](https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.linalg.sqrtm.html)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "G = # ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Double check that multiplying the square root matrix to itself returns us back to the original matrix of factor variances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.matmul(G,G) - F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: calculate $\\textbf{Q}$ and $\\textbf{Q}^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "# Q = GB'\n",
    "# Q should be a short and wide matrix\n",
    "Q = # ...\n",
    "Q.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "# Q' = BG\n",
    "# Q should be a tall and narrow matrix\n",
    "QT = # ...\n",
    "QT.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# notice we could also use the transpose of Q to get Q'\n",
    "QT - Q.transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: Include portfolio holdings\n",
    "\n",
    "So the original formula of \n",
    "$h^TBFB^Th$ became  \n",
    "$h^TBGGB^Th$, where $GG = F$.  \n",
    "\n",
    "And then, if we let $Q^T=BG$ and $Q = GB^T$:  \n",
    "$h^TQ^TQh$\n",
    "\n",
    "\n",
    "Let $R = Q h$ and $R^T = h^T Q^T$:  \n",
    "\n",
    "The risk term becomes:  \n",
    "$R^TR$, where $R^T=h^TQ$ and $R=Q^Th$\n",
    "\n",
    "So an important point here is that we don't want to multiply $Q^TQ$ itself, because this creates the large N by N matrix. We want to multiply $h^TQ^T$ and $Qh$ separately, creating vectors of length k (k is number of risk factors)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "# R = Qh\n",
    "R = # ...\n",
    "R.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "# R' = Q'h'\n",
    "RT = # ...\n",
    "RT.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notice how we avoided creating a full N by N matrix\n",
    "\n",
    "Also, notice that if we have Q, we can take its transpose to get $Q^T$ instead of doing the matrix multiplication.  \n",
    "\n",
    "Similarly, if we have R, which is a vector, we notice that $R^TR$ is the same as taking the dot product.  In other words, it's squaring each element in the vector R, and adding up all the squared values.\n",
    "\n",
    "$R^TR = \\sum_{i}^{k}(r_i^2)$\n",
    "\n",
    "## Quiz: Put it all together: calculate common risk term efficiently"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: common risk term in term\n",
    "\n",
    "# TODO: calculate square root of F\n",
    "G = # ...\n",
    "\n",
    "# TODO: calculate Q\n",
    "Q = # ...\n",
    "\n",
    "# TODO: calculate R\n",
    "R = # ...\n",
    "\n",
    "# TODO: calculate common risk term\n",
    "common_risk = # ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specific Risk term\n",
    "\n",
    "The portfolio's variance that is specific to each asset is found by combining the holdings with the specific variance matrix: \n",
    "\n",
    "$h^TSh$, where $h^T$ is a 1 by N vector, S is an N by N matrix, and h is an N by 1 vector.\n",
    "\n",
    "Recall that S is a diagonal matrix, so all the off-diagonals are zero.  So instead of doing the matrix multiplication, we could save computation by working with the vector containing the diagonal values.\n",
    "\n",
    "$h^TSh = \\sum_i^{N}(h_i^2 \\times S_i)$ because $S$ is a diagonal matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## check the unit of measure of SpecRisk\n",
    "# Notice that these are in percent; multiply by .01 to get them back to decimals.aa\n",
    "universe['SpecRisk'][0:2] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: Specific Risk term\n",
    "\n",
    "Given specific risk (volatility), calculate specific variance.  First re-scale the specific risk data so that it's in decimal instead of percent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: specific variance : rescale it and then square to get specific variance\n",
    "specVar = # ...\n",
    "\n",
    "# TODO: specific risk term (include holdings)\n",
    "spec_risk_term = np.dot(\"\"\"<your code here>\"\"\", \"\"\"<your code here>\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximize portfolio returns\n",
    "\n",
    "Since the alpha vector $\\mathbf{\\alpha}$ is supposed to be indicative of future asset returns, when we look at a portfolio of assets, the weighted sum of these alphas $\\mathbf{\\alpha}^T \\mathbf{h}$ is predictive of the portfolio's future returns.  We want to maximize the portfolio's expected future returns, so we want to minimize the negative of portfolio's expected returns $-\\mathbf{\\alpha}^T \\mathbf{h}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO\n",
    "expected_return = # ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear price impact of trading\n",
    "\n",
    "Assume transaction cost is linearly related to the trade size as a fraction of the average daily volume.  Since we won't know the actual daily volume until the day that we're executing, we want to use past data as an estimate for future daily volume.  This would be kind of noisy if we simply use the prior day's daily volume, so we'd prefer a more stable estimate like a 30 day rolling average.\n",
    "\n",
    "A commonly used **estimate for linear market impact is that if a trade size is 1% of the ADV, this moves the price by 10 basis points (1/10,000).**\n",
    "\n",
    "$Trade size_{i,t}$ is the fraction of your trade relative to the average dollar volume estimated for that stock, for that day. \n",
    "\n",
    "$Trade_{i,t}$ = dollar amount to trade = $h_{t} - h_{t-1}$, which is the new holding of the asset minus the previous holding.\n",
    "\n",
    "$ADV_{i,t}$: (average dollar volume) is total dollar amount expected to be traded, based on a moving average of historical daily volume.\n",
    "\n",
    "$TradeSize_{i,t} = \\frac{Trade_{i,t}}{ADV_{i,t}}$: The size of the trade relative to the estimated daily volume.\n",
    "\n",
    "$\\% \\Delta Price_{i,t}$ = price change due to trading, as a fraction of the original price (it's a percent change)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll write out the ratio: change in price divided by the trade size.\n",
    "\n",
    "$ \\frac{\\% \\Delta price_{i,t}}{TradeSize_{i,t}} = \\frac{10 bps}{1\\%}$  \n",
    "\n",
    "$ \\frac{\\% \\Delta price_{i,t}}{TradeSize_{i,t}} = \\frac{10/10^4}{1/100}$\n",
    "\n",
    "$ \\frac{\\% \\Delta price_{i,t}}{TradeSize_{i,t}} = \\frac{10^{-3}}{10^{-2}}$\n",
    "\n",
    "$ \\frac{\\% \\Delta price_{i,t}}{TradeSize_{i,t}} = 10^{-1}$\n",
    "\n",
    "Now we'll move things around to solve for the change in price.\n",
    "\n",
    "$\\%  \\Delta price_{i,t} = 10^{-1} \\times TradeSize_{i,t}$\n",
    "\n",
    "We defined TradeSize to be the Trade divided by ADV.\n",
    "\n",
    "$\\%  \\Delta price_{i,t} = 10^{-1} \\times \\frac{Trade_{i,t}}{ADV_{i,t}}$\n",
    "\n",
    "Note that Trade is the current position minus the prior day's position\n",
    "\n",
    "$\\%  \\Delta price_{i,t} = 10^{-1} \\times \\frac{h_{i,t} - h_{i,t-1}}{ADV_{i,t}}$\n",
    "\n",
    "For convenience, we'll combine the constant $10^{-1}$ and $\\frac{1}{ADV_{i}}$ and call it lambda $\\lambda_{i}$\n",
    "\n",
    "$\\%  \\Delta price_{i,t} = \\lambda_{i,t} \\times (h_{i,t} - h_{i,t-1})$ where $\\lambda_{i,t} = 10^{-1}\\times \\frac{1}{ADV_{i,t}} = \\frac{1}{10 \\times ADV_{i,t}}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that since we're dividing by $ADV_{i,t}$, we'll want to handle cases when $ADV_{i,t}$ is missing or zero.  In those instances, we can set $ADV_{i,t}$ to a small positive number, such as 10,000, which, in practice assumes that the stock is illiquid.\n",
    "\n",
    "Represent the market impact as $\\Delta price_{i} = \\lambda_{i} (h_{i,t} - h_{i,t-1})$.  $\\lambda_{i}$ incorporates the $ADV_{i,t}$.  Review the lessons to see how to do this.\n",
    "\n",
    "Note that since we're dividing by $ADV_{i,t}$, we'll want to handle cases when $ADV_{i,t}$ is missing or zero.  In those instances, we can set $ADV_{i,t}$ to a small positive number, such as 10,000, which, in practice assumes that the stock is illiquid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "If the ADV field is missing or zero, set it to 10,000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: if missing, set to 10000\n",
    "universe.loc[np.isnan(universe[\"\"\"<your code here>\"\"\"]), 'ADTCA_30'] = \"\"\"<your code here>\"\"\" ## assume illiquid if no volume information\n",
    "\n",
    "# TODO: if zero, set to 10000\n",
    "universe.loc[universe['ADTCA_30'] == \"\"\"<your code here>\"\"\", 'ADTCA_30'] = \"\"\"<your code here>\"\"\" ## assume illiquid if no volume information\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: calculate Lambda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "adv = # ...\n",
    "Lambda = # ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: transaction cost term\n",
    "\n",
    "Transaction cost is change in price times dollar amount traded.  For a single asset \"i\":\n",
    "\n",
    "$tcost_{i,t} = (\\% \\Delta price_{i,t}) \\times (DollarsTraded_{i,t})$\n",
    "\n",
    "$tcost_{i,t} = (\\lambda_{i,t} \\times (h_{i,t} - h_{i,t-1}) ) \\times (h_{i,t} - h_{i,t-1})$\n",
    "\n",
    "Notice that we can simplify the notation so it looks like this:\n",
    "\n",
    "$tcost_{i,t} = \\lambda_{i,t} \\times (h_{i,t} - h_{i,t-1})^2$\n",
    "\n",
    "The transaction cost term to be minimized (for all assets) is:\n",
    "\n",
    "$tcost_{t} = \\sum_i^{N} \\lambda_{i,t} (h_{i,t} - h_{i,t-1})^2$  \n",
    "where $\\lambda_{i,t} = \\frac{1}{10\\times ADV_{i,t}}$\n",
    "\n",
    "For matrix notation, we'll use a capital Lambda, $\\Lambda_{t}$, instead of the lowercase lambda $\\lambda_{i,t}$.\n",
    "\n",
    "$tcost_{t} = (\\mathbf{h}_{t} - \\mathbf{h}_{t-1})^T \\mathbf{\\Lambda}_t (\\mathbf{h}_{t} - \\mathbf{h}_{t-1})$\n",
    "\n",
    "Note that we'll pass in a vector of holdings as a numpy array.  For practice, we'll use the h variable that is initialized to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "tcost = # ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## objective function\n",
    "\n",
    "Combine the common risk, idiosyncratic risk, transaction costs and expected portfolio return into the objective function.  Put this inside a function.\n",
    "\n",
    "Objective function is:  \n",
    "factor risk + idiosyncratic risk - expected portfolio return + transaction costs  \n",
    "$f(\\mathbf{h}) = \\frac{1}{2}\\kappa \\mathbf{h}_t^T\\mathbf{Q}^T\\mathbf{Q}\\mathbf{h}_t + \\frac{1}{2} \\kappa \\mathbf{h}_t^T \\mathbf{S} \\mathbf{h}_t - \\mathbf{\\alpha}^T \\mathbf{h}_t + (\\mathbf{h}_{t} - \\mathbf{h}_{t-1})^T \\mathbf{\\Lambda} (\\mathbf{h}_{t} - \\mathbf{h}_{t-1})$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Risk Aversion $\\kappa$\n",
    "\n",
    "The risk aversion term is set to target a particular gross market value (GMV), or to target a desired volatility.  In our case, we tried a few values of the risk aversion term, ran the backtest, and calculated the GMV.  Ideally, a quant who is just starting out may have a targeted GMV of 50 million.  A risk aversion term of $10^{-6}$ gets the GMV to be in the tens of millions.  A higher risk aversion term would decrease the GMV, and a lower risk aversion term would increase the GMV, and also the risk.  Note that this isn't necessarily a linear mapping, so in practice, you'll try different values and check the results.\n",
    "\n",
    "Also, in practice, you'd normally keep the risk aversion term constant, unless your fund is accepting more investor cash, or handling redemptions.  In those instances, the fund size itself changes, so the targeted GMV also changes.  Therefore, we'd adjust the risk aversion term to adjust for the desired GMV. \n",
    "\n",
    "Also, note that we would keep this risk aversion term constant, and not adjust it on a daily basis.  Adjusting the risk aversion term too often would result in unecessary trading that isn't informed by the alphas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "An important point is to think about what matrices can be multiplied independently of the vector of asset holdings, because those can be done once outside of the objective function.  The rest of the objective function that depends on the holdings vector will be evaluated inside the objective function multiple times by the optimizer, as it searches for the optimal holdings.  \n",
    "\n",
    "\n",
    "For instance, \n",
    "\n",
    "$\\mathbf{h}^T\\mathbf{BFB}^T\\mathbf{h}$ became  \n",
    "$\\mathbf{h}^T\\mathbf{BGGB}^T\\mathbf{h}$, where $\\mathbf{GG} = \\mathbf{F}$.  \n",
    "\n",
    "And then, if we let $\\mathbf{Q}^T=\\mathbf{BG}$ and $\\mathbf{Q} = \\mathbf{GB}^T$:  \n",
    "$\\mathbf{h}^T\\mathbf{Q}^T\\mathbf{Qh}$\n",
    "\n",
    "Let $\\mathbf{R} = \\mathbf{Q h}$ and $\\mathbf{R}^T = \\mathbf{h}^T \\mathbf{Q}^T$:  \n",
    "\n",
    "The risk term becomes:  \n",
    "$\\mathbf{R}^T\\mathbf{R}$, where $\\mathbf{R}^T=\\mathbf{h}^T\\mathbf{Q}$ and $\\mathbf{R}=\\mathbf{Q}^T\\mathbf{h}$\n",
    "\n",
    "* Can we pre-compute Q outside of the objective function?  \n",
    "* Can we pre-compute R outside of the objective function?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer\n",
    "Q doesn't depend on h, the holdings vector, so it can be pre-computed once outside of the objective function.\n",
    "\n",
    "R is created using h, the holdings vector.  This should be computed each time the objective function is called, not pre-computed beforehand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Risk Aversion parameter\n",
    "\n",
    "The risk aversion term is set to target a particular gross market value (GMV), or to target a desired volatility.  \n",
    "\n",
    "The gross market value is the dollar value of the absolute value of the long and short positions.\n",
    "\n",
    "$ GMV = \\sum_i^N(|h_{i,t}|)$\n",
    "\n",
    "When we think about what it means to take more risk when investing, taking bigger bets with more money is a way to take on more risk.  So the risk aversion term controls how much risk we take by controlling the dollar amount of our positions, which is the gross market value.\n",
    "\n",
    "In our case, we tried a few values of the risk aversion term, ran the backtest, and calculated the GMV.  Ideally, a quant who is just starting out may have a targeted book size of 50 million.  In other words, they try to keep their GMV around 50 million.  \n",
    "\n",
    "A risk aversion term of $10^{-6}$ gets the GMV to be in the tens of millions.  A higher risk aversion term would decrease the GMV, and a lower risk aversion term would increase the GMV, and also the risk.  Note that this isn't necessarily a linear mapping, so in practice, you'll try different values and check the results.\n",
    "\n",
    "Also, in practice, you'd normally keep the risk aversion term constant, unless your fund is accepting more investor cash, or handling redemptions.  In those instances, the fund size itself changes, so the targeted GMV also changes.  Therefore, we'd adjust the risk aversion term to adjust for the desired GMV. \n",
    "\n",
    "Also, note that we would keep this risk aversion term constant, and not adjust it on a daily basis.  Adjusting the risk aversion term too often would result in unnecessary trading that isn't informed by the alphas.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Risk aversion\n",
    "risk_aversion=1.0e-6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: define objective function\n",
    "\n",
    "Combine the common risk, idiosyncratic risk, transaction costs and expected portfolio return into the objective function.  Put this inside a function.\n",
    "\n",
    "Objective function is:  \n",
    "factor risk + idiosyncratic risk - expected portfolio return + transaction costs  \n",
    "$f(\\mathbf{h}) = \\frac{1}{2}\\kappa \\mathbf{h}_t^T\\mathbf{Q}^T\\mathbf{Q}\\mathbf{h}_t + \\frac{1}{2} \\kappa \\mathbf{h}_t^T \\mathbf{S} \\mathbf{h}_t - \\mathbf{\\alpha}^T \\mathbf{h}_t + (\\mathbf{h}_{t} - \\mathbf{h}_{t-1})^T \\mathbf{\\Lambda} (\\mathbf{h}_{t} - \\mathbf{h}_{t-1})$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(h): \n",
    "    # TODO: define the objective function, where h is the vector of asset holdings\n",
    "    f = # ...\n",
    "    return(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient\n",
    "\n",
    "Before, when we used cvxpy, we didn't have to calculate the gradient, because the library did that for us.\n",
    "\n",
    "Objective function is:  \n",
    "factor risk + idiosyncratic risk - expected portfolio return + transaction costs  \n",
    "$f(\\mathbf{h}) = \\frac{1}{2}\\kappa \\mathbf{h}^T\\mathbf{Q}^T\\mathbf{Qh} + \\frac{1}{2} \\kappa \\mathbf{h}^T \\mathbf{S h} - \\mathbf{\\alpha^T h} + (\\mathbf{h}_{t} - \\mathbf{h}_{t-1})^T \\Lambda (\\mathbf{h}_{t} - \\mathbf{h}_{t-1})$\n",
    "\n",
    "\n",
    "Let's think about the shape of the resulting gradient. The reason we're interested in calculating the derivative is so that we can tell the optimizer in which direction, and how much, it should shift the portfolio holdings in order to improve the objective function (minimize variance, minimize transaction cost, and maximize expected portfolio return).  So we want to calculate a derivative for each of the N assets (about 2000+ in our defined universe).  So the resulting gradient will be a row vector of length N.\n",
    "\n",
    "The gradient, or derivative of the objective function, with respect to the portfolio holdings h, is:  \n",
    "\n",
    "$f'(\\mathbf{h}) = \\frac{1}{2}\\kappa (2\\mathbf{Q}^T\\mathbf{Qh}) + \\frac{1}{2}\\kappa (2\\mathbf{Sh}) - \\mathbf{\\alpha} + 2(\\mathbf{h}_{t} - \\mathbf{h}_{t-1}) \\mathbf{\\Lambda}$\n",
    "\n",
    "We can check that each of these terms is a row vector with one value for each asset (1 by N row vector)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Calculate the gradient of the common risk term:\n",
    "\n",
    "$\\kappa (\\mathbf{Q}^T\\mathbf{Qh})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: gradient of common risk term\n",
    "tmp = risk_aversion * np.matmul(\"\"\"<your code here>\"\"\", np.matmul(\"\"\"<your code here>\"\"\",\"\"\"<your code here>\"\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify that the calculation returns one value for each asset in the stock universe (about 2000+ )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Calculate gradient of idiosyncratic risk term\n",
    "\n",
    "$\\kappa (\\mathbf{Sh})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: idiosyncratic risk gradient\n",
    "tmp = \"\"\"<your code here>\"\"\" * \"\"\"<your code here>\"\"\" * \"\"\"<your code here>\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Calculate the gradient of the expected return\n",
    "\n",
    "$- \\mathbf{\\alpha} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: expected return gradient\n",
    "tmp = -\"\"\"<your code here>\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz\n",
    "\n",
    "Calculate the gradient of the transaction cost.\n",
    "\n",
    "$ 2(\\mathbf{h}_{t} - \\mathbf{h}_{t-1}) \\mathbf{\\Lambda}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# transaction cost\n",
    "tmp = 2 * \"\"\"<your code here>\"\"\" * \"\"\"<your code here>\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: Define gradient function\n",
    "\n",
    "Put this all together to define the gradient function.  The optimizer will use this to make small adjustments to the portfolio holdings.\n",
    "\n",
    "#### gradient (slightly cleaned up)\n",
    "\n",
    "We'll simplify the expression a bit by pulling the common $\\kappa$ out of the common risk and specific risk.  Also, the 1/2 and 2 cancel for both risk terms.\n",
    "\n",
    "$f'(\\mathbf{h}) = \\frac{1}{2}\\kappa (2\\mathbf{Q}^T\\mathbf{Qh}) + \\frac{1}{2}\\kappa (2\\mathbf{h}^T\\mathbf{S}) - \\mathbf{\\alpha} + 2(\\mathbf{h}_{t} - \\mathbf{h}_{t-1})\\cdot \\Lambda$\n",
    "\n",
    "becomes\n",
    "\n",
    "$f'(\\mathbf{h}) = \\kappa (\\mathbf{Q}^T\\mathbf{Qh} + \\mathbf{Sh}) - \\mathbf{\\alpha} + 2(\\mathbf{h}_{t} - \\mathbf{h}_{t-1}) \\mathbf{\\Lambda}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "def grad(x):\n",
    "    # TODO\n",
    "    g = # ...\n",
    "    return(np.asarray(g))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimizer\n",
    "\n",
    "Choose an optimizer.  You can read about these optimizers:\n",
    "\n",
    "* L-BFGS  \n",
    "* Powell\n",
    "* Nelder-Mead\n",
    "* Conjugate Gradient\n",
    "\n",
    "In this [page about math optimization](http://scipy-lectures.org/advanced/mathematical_optimization/)\n",
    "\n",
    "Also read the [scipy.optimize documentation](https://docs.scipy.org/doc/scipy/reference/optimize.html)\n",
    "\n",
    "Pass in the objective function, prior day's portfolio holdings, and the gradient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "optimizer_result = scipy.optimize.fmin_l_bfgs_b(\"\"\"<your code here>\"\"\", \"\"\"<your code here>\"\"\", fprime=\"\"\"<your code here>\"\"\")\n",
    "\n",
    "h1 = optimizer_result[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_portfolio = pd.DataFrame(data = {\"Barrid\" : universe['Barrid'], \"h.opt\" : h1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_portfolio.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## risk exposures\n",
    "factor exposures times the portfolio holdings for each asset, gives the portfolio's exposure to the factors (portfolio's risk exposure).\n",
    "\n",
    "$\\mathbf{B}^T\\mathbf{h}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: risk exposures\n",
    "risk_exposures = np.matmul(\"\"\"<your code here>\"\"\", \"\"\"<your code here>\"\"\")\n",
    "\n",
    "# put this into a pandas series\n",
    "pd.Series(risk_exposures, index = colnames(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quiz: alpha exposures\n",
    "\n",
    "The portfolio's exposures to the alpha factors is equal to the matrix of alpha exposures times the portfolio holdings.  We'll use the holdings returned by the optimizer.\n",
    "\n",
    "$\\textbf{B}_{\\alpha}^T\\mathbf{h}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution: portfolio's alpha exposure\n",
    "alpha_exposures = np.matmul(\"\"\"<your code here>\"\"\", \"\"\"<your code here>\"\"\")\n",
    "\n",
    "# put into a pandas series\n",
    "pd.Series(alpha_exposures, index = colnames(B_alpha))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hints for the project\n",
    "\n",
    "You'll be putting this optimization code into functions so that you can call the optimizer in a loop, as the backtester walks through each day in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution notebook\n",
    "\n",
    "The solution notebook is [here]('./optimization_with_tcosts_solution.ipnyb')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}