Long/Short Global Macro Strategies with Target Beta Using the Three-Factor Model:

ABSTRACT— This project involves the construction of a Long/Short Global Macro strategy based on Fama-French Three-Factor Model with a target beta and the evaluation of its sensitivity to variation of Beta and the length of the estimation terms for the covariance matrix and the expected returns under different market scenarios. Several comparisons are drawn between different target betas as well as different term structures.

KEYWORDS— Global Macroeconomic Trends, Long/Short, French, Fama, Beta, Three-Factor Model, Trading Strategies, Portfolio Management

GOOGLE COLAB OVERVIEW— This is the complete design and implementation docs in a Google Colab notebook,

I. INTRODUCTION

A. Fama-French Three-Factor Model
Historically, Fama-French Three-Factor Model is regarded as a development of CAPM which explains a relationship between expected returns and risk factors. Sharpe, Lintner, and Black developed an asset pricing model referred as CAPM which illustrates expected returns on the securities are a positive linear function of market beta. The CAPM, developed theoretically, had an empirical success, and became the standard model for describing the cross-sectional structure of expected returns on equity. However, some researchers revealed there were some anomalies that cannot be explained by the CAPM. Banz found the market equity had a relationship between expected returns. Stocks with smaller market equity had higher rates of return, usually referred to as the small cap effect. Bhandari found leverage helps explain the cross-section of average stock returns. Stattman, Rosenberg, Reid and Lanstein found that average returns on U.S stocks were positively related to the ratio of a firm’s book value of common equity to its market value. In response to this criticism, Eugene Farmer, and Kenneth French, in a paper published in 1992, empirically showed that the four representative anomaly factors discovered at the time in the U.S. stock market: market equity, book-to-market ratio, leverage, and E/P (the inverse of P/E ratio), were aggregated into market equity and book-to-market ratio. They advanced their study and proposed Fama-French Three-Factor model which describes a cross section of average stock returns with three factors, market risk premium, market equity, and book-to-market ratio. Under the Fama-French Three-Factor model, the random return of security is given by the following formula,

With \(E[\epsilon_{i}]=0\), \(E[r_{i}]=0\), \(E[r{m}]=R_{m}\), \(E[r_{i}^{SMB}]=R^{SMB}\) and \(E[r_{i}^{HML}]=R^{HML}\), the expected return on a given security can be written as, Where \(r_{i}\) is the random return on a given security, \(r_{f}\) is the risk-free rate, \(r_{M}\) is the market return, and \(\beta_{i}^{MKT}\), \(\beta_{i}^{SMB}\), \(\beta_{i}^{HML}\) are the sensitivity measure for the risk premium of the market portfolio, the risk factor of market equity, and the risk factor of the book-to-market ratio respectively.

B. Markowitz Portfolio
Markowitz portfolio theory also known as modern portfolio theory is a theory on how risk-averse investors can construct portfolio to maximize expected return based on a given level of market risk. The theory can also be used to construct a portfolio that minimize risk for a given level of expected return,
where \(\omega\) is the vector of weights of the portfolio securities, \(\Sigma\) is the covariance matrix , \(\rho\) is the expected return, \(\rho_{T}\) is the target return.

C. Linear Regression
An approach for predicting a quantitative response \(Y\) and the basis of multiple predictor variable \(X_{j}\) that assume an approximately linear relationship between \(X_{j}\) and \(Y\). For a model with \(p\) predictors, the linear regression takes the form, where \(X_{j}\) is the \(j^{th}\) predictor and \(\beta_{j}\) qualifies the relationship between that predictor and the response. Given estimates for \(\beta\)'s, it can generate predictors using the model, where we estimate these parameters by minimizing the residual sum of squares,

II. INVESTMENT UNIVERSE AND BACKTESTING

A. Data
We used the 12 ETFs belowfrom March 1, 2007 to June 30, 2020:

1. CurrencyShares Euro Trust (FXE)
2. iShares MSCI Japan Index (EWJ)
3. SPDR GOLD Trust (GLD)
4. PowerShares NASDAQ-100 Trust (QQQ)
5. SPDR S&P500 (SPY)
6. iShares Lehman Short Treasury Bond (SHV)
7. PowerShares DB Agriculture Fund (DBA)
8. United States Oil Fund LP (USO)
9. SPDR S&P Biontech (XBI)
10. iShares S&P Latin America 40 Index (ILF)
11. iShares MSCI Pacific ex-Japan Index Fund (EPP)
12. SPDR DJ Euro Stoxx 50 (FEZ)

The S&P500 (SPY) was chosen to be the analysis benchmark. Lastly, the data used to construct the Fama-French Three-Factor Model is quoted from Ken French’s website for the factors’ historical values.


B. Investment Horizon
The investment horizon was divided into the following sub-periods,

1. Pre-Subprime Crisis : March 22, 2007 – March 3, 2008
2. During Subprime Crisis : March 3, 2008 – September 10, 2010
3. Post-Subprime Crisis : September 10, 2010 – January1, 2015
4. Pre-COVID-19 Pandemic : January 1, 2015 – March 9, 2020
5. During COVID-19 Pandemic : March 9, 2020 – October 30, 2020

C. Backtesting
Individual backtests were executed for each sub-period to compare strategies. We compared with different perspectives.

1. Impact on Beta Target: Compared our strategy in terms of target beta in the same sub-period. Changed target beta \(\beta_{T}^{m}\), and compare performance.
2. Impact of various term structure given Beta: Compare portfolios’ performance with different term structures and fixed beta. \(S_j^{i}\) represents a term structure with \(i\) days lookback period to estimate the expected return and \(j\) days lookback period to estimate the covariance matrix.

We also ran a backtest across the full investment horizon from March 1st, 2007 to November 30th, 2020.

III. INVESTMENT STRATEGY

A. Objective Function
We consider the following investment strategy,

where,

• \(\omega_{i}\): weight allocated to each security \(S_{i}\),
• \(\omega\) is a vector of weights,
• \(\rho\): vector of expected returns of security \(S_{i},
• \(\Sigma\): The covariance matrix between securities reuturns derived from the factor model,
• \(\omega_{p}\): composition of a reference Portfolio, the previous portfolio when rebalancing the portfolio,
• \(\beta_{i}^{m}=\frac{cov(r_{i},r_{M})}{\sigma^{2}(r_{M})}\): the Beta of security \(S_{i}\) as defined in the CAPM model,
• \(\beta_{T}^{m}\): the portfolio's target beta.

B. Term Structures
We investigated the following term structures,
\begin{equation} \begin{cases} S_{60}^{60} \\ S_{120}^{60} \\ S_{60}^{90} \\ S_{120}^{90} \\ S_{60}^{120} \\ S_{120}^{120} \end{cases} \end{equation}
C. Target Beta
We investigated combinations of the term structures and the following target \(\beta\)'s,
\begin{equation} \beta_{T}^{m}=\begin{cases} -1.0 \\ -0.5 \\ 0.5 \\ 1.0 \\ 1.5 \\ 2.0 \end{cases} \end{equation} D. Performance Metrics
The trading strategy pipeline performs computations of the following performance metrics,

• Cumulated Return
• Annual Arithmetic Mean / Geometric Mean Return
• Annual Min Return
• Max 10-days Drawdown
• Sharpe Ratio

E. Risk Metrics
The trading strategy pipeline also performs computations of the following risk metrics,

• Volatility
• Daily VaR
• Annual VaR
• Modified VaR
• Annual CVaR
• Skewness
• Kurtosis

IV. RESULTS & DISCUSSION

A. Return on Investment
The returns and log-returns of the portfolio securities are summarized below,

B. PnL
The evolution of cumulated daily Profit and Loss as summing with an initial investment of $100 on the first allocation date is as follows,







C. Distribution of Daily Returns
The following plots show the distribution of daily returns for different term structures and \(\beta\)'s (Click on the right/left side of each slideshow to go to the next/previous plot),






D. Strategy Performance and Risk Metrics

i. Before the Subprime Crisis

ii. During the Subrime Crisis

iii. After the Subprime Crisis

iv. Before COVID-19

v. During COVID-19

vi. Full Investment Horizon

ACKNOWLEDGEMENTS
I would like to thank Professor Ndiaye for his invaluable guidance through his teachings in FE 630: Portfolio Theory & Applications, and his continued support in identifying the pain points and potential developments in the project proposal, as well as the developers of Google's Collaboratory for providing a streamlined platform for the testing and execution of the models.

AUTHORS

• Theo Dimitrasopoulos– MSc. Financial Engineering, Stevens Institute of Technology 2021; Bachelor of Science in Engineering, Structural Engineering, Princeton University 2017; tdimitr1(at)stevens(dot)edu
• Yuki Homma– MSc. Financial Engineering, Stevens Institute of Technology 2020; yhomma(at)stevens(dot)edu

REFERENCES

• [1] Eugene F. Fama, Kenneth R. French, “The Cross-Section of Expected Stock Return” ,The Journal of Finance, vol. 47, No.2 , 1992, pp. 427–465.
• [2] Eugene F. Fama, Kenneth R. French, “Common risk factors in the returns on stocks and bonds” ,Journal of Financial Economics, vol. 33, No.1 , 1993, pp. 3–56.
• [3] Eugene F. Fama, Kenneth R. French, “Size and Book-to-Market Factors in Earnings and Return” ,The Journal of Finance, vol. 50, No.1 , 1995, pp. 131–155.
• [4] Eugene F. Fama, Kenneth R. French, “Multifactor Explanation of Asset Pricing Anomalies” ,The Journal of Finance, vol. 51, No.1 , 1996, pp. 55–84.
• [5] Francis, Kim, “Modern Portfolio Thoery”, Wiley Finance

SOURCE CODE

# -*- coding: utf-8 -*-
"""global_macro.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1fmXkzxM2Q6rcfWPwg6OAUw3pf5HHxMsO

# **Long/Short Global Macro Strategies with Target \(\beta\) Using the Three-Factor Model**
Authors: Theo Dimitrasopoulos✝, Yuki Homma✝

Advisor: Papa Momar Ndiaye✝

✝ Department of Financial Engineering; Stevens Institute of Technology Babbio School of Business

## **Authors**

**Final Project**

FE630: Modern Portfolio Theory & Applications

### Yuki Homma

Contribution: Presentation Preparation, Manuscript Generation

Stevens MS in Financial Engineering '21

### Theo Dimitrasopoulos

Contribution: Quant Analysis, Programming

Stevens MS in Financial Engineering '21 | Princeton '17

theonovak.com | +1 (609) 933-2990

## **Introduction**

### Investment Strategy

We build a Long-Short Global Macro Strategy with a Beta target using a factor-based model and evaluate its sensitivity to variations of Beta.

To optimize the portfolios, we deploy the following strategies:
1. Maximize the return of the portfolio subject to a constraint of target Beta, where Beta is the single-factor market risk measure. This allows us to evaluate the sensitivity of the portfolios to variations of Beta. The portfolio is re-optimized (weight recalibration) every week for the investment horizon between March 2007 to the end of October 2020.
2. Minimum variance with a target return.

### Optimization Problem:

The strategy aims to maximize return with a certain Target Beta under constraints.

It is defined as,

\begin{cases}
\max\limits_{{\omega ∈ ℝ^{n}}}\rho^{T}\omega-\lambda(\omega-\omega_{p})^{T}\Sigma(\omega-\omega_{p})\\
\sum_{i=1}^{n} \beta_{i}^{m}\omega_{i}=\beta_{T}^{m}\\
\sum_{i=1}^{n} \omega_{i}=1, -2\leq\omega_{i}\leq2
\end{cases}

$\Sigma$ is the the covariance matrix between the securities returns (computed from
the Factor Model), $\omega_{p}$ is the composition of a reference Portfolio (the previous Portfolio when rebalancing the portfolio and $\omega_{p}$ has all its components equal to $1/n$ for the first allocation) and $\lambda$ is a small regularization parameter to limit the turnover;

$\beta_{i}^{m}=\frac{cov(r_{i},r_{M}}{\sigma^{2}(r_{M})}$ is the Beta of security $S_{i}$ as defined in the CAPM Model so that $\beta_{P}^{m}=\sum_{i=1}^{n}\beta_{i}^{m}\omega_{i}$ is the Beta of the Portfolio;

$\beta_{T}^{m}$ is the Portfolio's Target Beta, for example $\beta_{T}^{m}=-1$, $\beta_{T}^{m}=-0.5$, $\beta_{T}^{m}=0$, $\beta_{T}^{m}=0.5$, $\beta_{T}^{m}=1.5$.

### Equivalent Optimization Problem:

We can reformulate the optimization problem above to make the programming process more straightforward:

$(\omega-\omega_{p})^{T}\Sigma(\omega-\omega_{p})\rightarrow$

$=(\omega-\omega_{p})^{T}\Sigma\omega-(\omega-\omega_{p} )^{T}\Sigma\omega_{p}$

$=\omega^{T} \Sigma\omega-2(\omega^{T} \Sigma\omega_{p})+\omega_{p}^{T}\Sigma \omega_{p}$

We simplify,
- $d=\rho-2\lambda\Sigma\omega_{p}$
- $P=\lambda\Sigma$

Finally,

$\max\limits_{{\omega ∈ ℝ^{n}}}(\rho-2\lambda\Sigma\omega_{p} )^{T} \omega-\lambda\omega^{T}\Sigma\omega+\lambda\omega_{p}^{T}\Sigma\omega_{p}=\max\limits_{{\omega ∈ ℝ^{n}}}d^{T}\omega-\omega^{T}P\omega$

---

The following formulation is equivalent,

\begin{cases}
\max\limits_{{\omega ∈ ℝ^{n}}}d^{T}\omega-\omega^{T}P\omega\\
\sum_{i=1}^{n} \beta_{i}^{m}\omega_{i}=\beta_{T}^{m}\\
\sum_{i=1}^{n} \omega_{i}=1, -2\leq\omega_{i}\leq2
\end{cases}
- $\Sigma$ is the the covariance matrix between the returns of the portfolio assets;
- $\omega_{p}$ is the composition of a reference Portfolio:
  - When rebalancing the portfolio, $\omega_{p}$ is the previous portfolio
  - $\omega_{p}$ has all its components equal to $1/n$ for the first allocation
- $\lambda$ is a regularization parameter to limit the turnover
- $\beta_{i}^{m}=\frac{cov(r_{i},r_{M}}{\sigma^{2}(r_{M})}$ is the Beta of security $S_{i}$ as defined in the CAPM Model s.t. $\beta_{P}^{m}=\sum_{i=1}^{n}\beta_{i}^{m}\omega_{i}$ is the portfolio Beta
- $\beta_{T}^{m}$ is the Portfolio's Target Beta.

### Fama French Three-Factor Model

A three-factor model proposed by Fama and French in 1993 includes not only market excess return, but capitalization size and a book-to-market ratio as influencing factors.

The random return of a given security is given by the formulas (equivalent),

\begin{equation}
\boxed{r = r_{f}+\beta_{1}(r_{m}-r_{f})+\beta{2}(SMB)+\beta_{3}(HML)+\epsilon}
\end{equation}

\begin{equation}
\boxed{R_{i}-r_{f}=\alpha_{i}+\beta{i}^{MKT}(R_{M}-r_{f})+\beta_{i}^{SMB}R_{SMB}+\beta_{i}^{HML}R_{HML}}
\end{equation}

- rSMB represents small size variables minus big one
- rHML represents high minus low in book value to equity to book value to the market.

### ETF Data

The following ETFs represent the investment Universe of our portfolios. They range from the S&P500 to ETFs representing all continents such as Europe, Asia and Africa and asset types such as bonds, stocks, and commodities.

1. CurrencyShares Euro Trust (FXE)
2. iShares MSCI Japan Index (EWJ)
3. SPDR GOLD Trust (GLD)
4. Powershares NASDAQ-100 Trust (QQQ)
5. SPDR S&P500 (SPY)$^*$
6. iShares Lehman Short Treasury Bond (SHV)
7. PowerShares DB Agriculture Fund (DBA)
8. United States Oil Fund LP (USO)
9. SPDR S&P Biotech (XBI)
10. iShares S&P Latin America 40 Index (ILF)
11. iShares MSCI Pacific ex-Japan Index Fund (EPP)
12. SPDR DJ Euro Stoxx 50 (FEZ)

From this universe, we have created portfolios by utilizing the Three-Factor Fama-French model. The investment portfolio that we created is compared to the following benchmark portfolios:

1.	The Market Portfolio (S&P500)

The dataset includes daily price data between March 1st, 2007 to October 31th, 2020. We choose this investment horizon to match the Fama-French Factor data currently available.

We have used three different look-back periods, which we have defined as: A. Short Term – 60 Days B. Medium Term – 120 Days C. Long Term – 200 Days To calculate the risk-return parameters of then portfolio we have used the target Beta as -1, -0.5, 0, 0.5, 1 and 1.5. The rebalance period is kept to one week as specified in the project.

*$^*$The SPY market portfolio is the chosen benchmark*

## **Setup**

### Environment
"""

# -*- coding: utf-8 -*-

# ENVIRONMENT CHECK:
import sys, os, inspect, site, pprint
# Check whether in Colab:
IN_COLAB = 'google.colab' in sys.modules
if IN_COLAB == True:
  print('YES, this is a Google Colaboratory environment.')
else:
  print('NO, this is not a Google Colaboratory environment.')
print(' ')

# Python installation files:
stdlib = os.path.dirname(inspect.getfile(os))
python_version = !python --version
print('Python Standard Library is located in:\n' + stdlib)
print(' ')
print('This environment is using {}'.format(str(python_version[0])))
print(' ')
print('Local system packages are located in:')
pprint.pprint(site.getsitepackages())
print(' ')
print('Local user packages are located in:\n' + site.getusersitepackages())

# Installed packages:
#!pip list -v
#!pip list --user -v

# Mount Google Drive:
if IN_COLAB:
  from google.colab import drive
  drive.mount('/content/drive', force_remount=True)

# Define Paths:
if IN_COLAB:
  graphs_dir = '/content/drive/MyDrive/Colab Notebooks/FE630_Final/report/graphics/'
  data_dir = '/content/drive/MyDrive/Colab Notebooks/FE630_Final/src/data/'
  source_dir = '/content/drive/MyDrive/Colab Notebooks/FE630_Final/src/'
else:
  graphs_dir = 'C:/Users/theon/GDrive/Colab Notebooks/FE630_Final/report/graphics/'
  data_dir = 'C:/Users/theon/GDrive/Colab Notebooks/FE630_Final/src/data/'
  source_dir = '/content/drive/MyDrive/Colab Notebooks/FE630_Final/src/'

"""### Packages

#### Uninstall Packages:
"""

# UNINSTALL PACKAGES:
#!pip uninstall pandas -y
#!pip uninstall numpy -y
#!pip uninstall cvxopt -y
#!pip uninstall matplotlib -y
#!pip uninstall pandas-datareader -y
#!pip uninstall zipline -y
#!pip uninstall pyfolio -y
#!pip uninstall alphalens -y
#!pip uninstall empyrical -y
#!pip uninstall mlfinlab -y

"""#### Install Packages:"""

# INSTALL PACKAGES:
#!pip install pandas
#!pip install numpy
#!pip install cvxopt
#!pip install matplotlib
#!pip install pandas-datareader
#!pip install zipline
#!pip install pyfolio
#!pip install alphalens
#!pip install empyrical
#!pip install mlfinlab

"""#### Inspect Packages"""

#!pip list -v
#!pip list --user -v

"""#### Import Packages:"""

# Core:
import sys
import io
import os
from datetime import datetime

# Numerical:
import numpy as np
import pandas as pd
import scipy
import scipy.stats as ss
from scipy.stats import kurtosis,skew,norm
from scipy.optimize import minimize, least_squares
import sklearn
from sklearn import linear_model
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KernelDensity
from sklearn.preprocessing import normalize
import statsmodels.api as smf

# Data Processing:
import pandas_datareader.data as web
import pickle
import urllib.request
import zipfile

# Plotting:
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.mlab as mlab

# Optimization:
from cvxopt import matrix, solvers

# Bundles:
#import zipline
#import pyfolio as pf
#import alphalens
#import empyrical
#import mlfinlab

"""#### Settings"""

# Commented out IPython magic to ensure Python compatibility.
# General:
import warnings
warnings.filterwarnings('ignore')
# %matplotlib inline
get_ipython().run_line_magic('matplotlib', 'inline')

# Other:
#%load_ext zipline
#%reload_ext zipline
#!zipline ingest

"""#### Custom Packages:"""

sys.path.append('/content/drive/MyDrive/Colab Notebooks/FE630_Final/src/packages/alphalens/')
sys.path.append('/content/drive/MyDrive/Colab Notebooks/FE630_Final/src/packages/empyrical/')
sys.path.append('/content/drive/MyDrive/Colab Notebooks/FE630_Final/src/packages/mlfinlab/')
sys.path.append('/content/drive/MyDrive/Colab Notebooks/FE630_Final/src/packages/pyfolio/')
sys.path.append('/content/drive/MyDrive/Colab Notebooks/FE630_Final/src/packages/zipline/')
path = pd.DataFrame(sys.path)
path.T

"""## **Definitions**

### General Utilities
"""

# Generate random weights that sum up to 1:
def weights_randn(n):
    k = np.random.rand(n)
    return k / sum(k)

def check_array(arr):
    if len(np.array(arr).shape)==1:
        days = len(np.array(arr))
        cols = 1
    elif len(np.array(arr).shape)==2:
        days = np.array(arr).shape[0]
        cols = np.array(arr).shape[1]
    else:
        raise TypeError('Input should be 1-D np.array or pd.Series or a 2-D np.array.')
    return cols,days

def var_w(rho, lamb, Q, wp, beta_im_ ,beta_T):
  def constrain1(w):
    return np.dot(beta_im_,w)-beta_T

  def constrain2(w):
    return np.sum(w)-1

  cons = [{'type':'eq', 'fun': constrain1},
          {'type':'eq', 'fun': constrain2}]
  bnds = scipy.optimize.Bounds(-2.0, 2.0, keep_feasible = True)

  def f(w):
    return -rho.dot(w) + lamb*(w-wp).dot(Q.dot(w-wp))

  w0 = np.array([1/12]*12)
  res = minimize(f, w0, method='SLSQP', bounds=bnds, constraints=cons,
                  tol=1e-9)
  return res.x

"""### Data Retrieval/Processing"""

# Save data:
def save_data(df, file_name, dir_name=data_dir):
    if not os.path.exists(dir_name):
        os.mkdir(dir_name)
    # Save results to a csv file
    df.to_csv(dir_name + file_name + '.csv', index=True)
    print('Successfully saved {}.csv. in {}'.format(file_name, dir_name + file_name + '.csv'))

# Download and prepare Fama French data:
def fama_french(frequency, no_factors):
  if frequency == 'annual':
    date_format = '  %Y'
    if no_factors == 3:
      ff_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_CSV.zip'
      filename = 'F-F_Research_Data_Factors'
    elif no_factors == 5:
      ff_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_5_Factors_2x3_CSV.zip'
      filename = 'F-F_Research_Data_5_Factors_2x3'
    else:
      print('Please choose 3 or 5 for the 3- and 5-Factor Model respectively.')
  elif frequency == 'monthly':
    date_format = '%Y%m'
    if no_factors == 3:
      ff_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_CSV.zip'
      filename = 'F-F_Research_Data_Factors'
    elif no_factors == 5:
      ff_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_5_Factors_2x3_CSV.zip'
      filename = 'F-F_Research_Data_5_Factors_2x3'
    else:
      print('Please choose 3 or 5 for the 3- and 5-Factor Model respectively.')
  elif frequency == 'weekly':
    date_format = '%Y%m%d'
    if no_factors == 3:
      ff_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_weekly_CSV.zip'
      filename = 'F-F_Research_Data_Factors_weekly'
    elif no_factors == 5:
      print ('No weekly data available for the 5-Factor Model.')
    else:
      print('Please choose 3 or 5 for the 3- and 5-Factor Model respectively.')
  elif frequency == 'daily':
    date_format = '%Y%m%d'
    if no_factors == 3:
      ff_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_daily_CSV.zip'
      filename = 'F-F_Research_Data_Factors_daily'
    elif no_factors == 5:
      ff_url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_5_Factors_2x3_daily_CSV.zip'
      filename = 'F-F_Research_Data_5_Factors_2x3_daily'
    else:
      print('Please choose 3 or 5 for the 3- and 5-Factor Model respectively.')
  else:
    print('Please choose between annual, monthly, weekly or daily for the frequency.')

  urllib.request.urlretrieve(ff_url, data_dir + filename + '.zip')
  zip = zipfile.ZipFile(data_dir + filename + '.zip', 'r')
  with zipfile.ZipFile(data_dir + filename + '.zip', 'r') as zip_ref:
    zip_ref.extractall(data_dir)
    zip.close()

  try:
    ff_factors = pd.read_csv(data_dir + filename + '.CSV', skiprows = 3, index_col = 0)
  except ValueError:
    ff_factors = pd.read_csv(data_dir + filename + '.CSV', skiprows = 3, index_col = 0)
  ff_row = ff_factors.isnull().any(1).nonzero()[0][0]
  try:
    ff_factors = pd.read_csv(data_dir + filename + '.CSV', skiprows = 3, index_col = 0)
  except ValueError:
    ff_factors = pd.read_csv(data_dir + filename + '.csv', skiprows = 3, index_col = 0)
  ff_factors = ff_factors.iloc[:-1]
  if frequency == 'annual':
    ff_factors = ff_factors.iloc[1134:,]
  elif frequency == 'monthly':
    ff_factors = ff_factors.iloc[0:1131,]
  else:
    pass
  ff_factors = ff_factors.dropna()
  ff_factors.index = pd.to_datetime(ff_factors.index, format=date_format)
  ff_factors.index = ff_factors.index + pd.offsets.MonthEnd()
  return ff_factors

"""### Risk/Performance Metrics"""

def PnL(arr,P = 1000000):
    cols,days = check_array(arr)
    data = np.array(arr).reshape(days, cols)
    ret = []
    s = (np.array([1.0 for _ in range(cols)]))*P
    for i in range(days):
        s += data[i,:]*s
        ret.append(s.copy())
    return np.array(ret)

def geom_mean(arr):
    cols,days = check_array(arr)
    data = np.array(arr).reshape(days, cols)
    return np.power(np.prod(1+data,axis=0),1/days)-1

def MaxDrawdown(arr, n=10):
    cols,days = check_array(arr)
    data = np.array(arr).reshape(days, cols)
    D_ = []
    d_ = []
    for day in range(n,days):
        returns = pd.DataFrame(1+data[(day-n):day,:]).cumprod(axis = 0)
        D = returns.cummax(axis=0)-returns
        d = np.array(D)/(np.array(D+returns))
        D_.append(np.max(np.array(D),axis=0))
        d_.append(np.max(np.array(d),axis = 0))
    return np.max(np.array(D_),axis=0),np.max(np.array(d_),axis=0)

def Volatility(arr,yearly=False):
    cols,days = check_array(arr)
    data = np.array(arr).reshape(days, cols)
    if yearly:
        return np.sqrt(np.var(data,axis=0))
    else:
        return np.sqrt((252/days)*np.sum((data-np.mean(data,axis=0))**2,axis=0))

def Sharpe(arr,rf,yearly = False):
    cols,days = check_array(arr)
    c,row = check_array(rf)
    if not days == row:
        raise RuntimeError('length of columns of inputs do not match (%s, %s).'% (days,row))
    data = np.array(arr).reshape(days, cols)
    r = np.array(rf).reshape(days,1)*250
    ER = np.power(np.product(1+data,axis=0),250/days)-np.mean(r,axis=0)-1
    return ER/Volatility(data)

def Kurt(arr):
    cols,days = check_array(arr)
    data = np.array(arr).reshape(days, cols)
    return ss.kurtosis(data,axis=0)

def Skew(arr):
    cols,days = check_array(arr)
    data = np.array(arr).reshape(days, cols)
    return ss.skew(data,axis=0)

def VaR(arr,q):
    cols,days = check_array(arr)
    data = np.array(arr).reshape(days, cols)
    tmp = np.sort(data,axis=0)
    n = int(np.around((1-q)*days))
    return -tmp[max(0,n-1),:]

def CVaR(arr,q):
    cols,days = check_array(arr)
    data = np.array(arr).reshape(days, cols)
    tmp = np.sort(data,axis=0)
    # print(tmp)
    n = int(np.around((1 - q) * days))
    return np.mean(-tmp[0:max(0, n - 1),:],axis=0)

def Summary(arr,RF, q=0.99):
    result = arr
    cols,days = check_array(result)
    print('Last PnL after %s: ' % days,PnL(result)[-1,:])
    # Geometric mean
    print('Geometric mean', geom_mean(result))
    # min
    print('Daily min', np.min(result, axis=0))
    # max drawdown
    print('max drawdown: ', MaxDrawdown(result))
    # Vol
    print('Volatility', Volatility(result))
    print('Sharp ratio: ', Sharpe(result, RF))
    print('Mean sharp: ', np.mean(Sharpe(result, RF), axis=0))
    print('Kurt: ', Kurt(result))
    print('Skewness: ', Skew(result))
    print('%s VaR %s days: ' % (q,days), VaR(result,q))
    print('%s CVaR %s days: ' % (q, days), CVaR(result, q))

"""### Backtesting"""

def backtesting(ret_etf, ff_factors, return_period, variance_period, lamb, beta_tm):
  port_returns = []
  omegas = []
  omega_p = np.array([1/12] *12)
  look_back =  max(return_period,
                   variance_period)
  next_chang_date = look_back - 1
  for i in range(len(ret_etf)):
    omegas.append(omega_p)
    today_return = np.asarray(ret_etf.iloc[i,:])
    pr = np.dot(omega_p,today_return)
    port_returns.append(pr)
    if i == next_chang_date:
      omega_p = omega(
          ret_r = ret_etf.iloc[i+1-return_period:i+1],
          factor_r =ff_factors.iloc[i+1-return_period:i+1],
          return_v = ret_etf.iloc[i+1-variance_period:i+1],
          factor_v = ff_factors.iloc[i+1-variance_period:i+1],
          lamb_ = lamb,
          beta_tm_ = beta_tm,
          wp_ = omega_p)
      next_chang_date += 5

    else:
        continue

  return port_returns,omegas

"""### Analytics"""

def analytics(X,rf,confidenceLevel,position):
  cum_ret_day=np.cumprod((X+1))
  cum_ret_annual = (np.power(cum_ret_day.iloc[-1,0],1/len(X)))**250
  arith_mean_ret_annual=np.mean(X)*250
  geom_mean_ret_annual=(np.power(cum_ret_day.iloc[-1,0],1/len(X))-1)*250
  min_ret_annual = np.min(X)*250
  p_v =np.cumprod((X+1))*100
  p_v_extend = pd.DataFrame(np.append([p_v.iloc[0,0]]*9,p_v))
  rolling_window_max = p_v_extend.rolling(window=10).max()
  ten_day_drawdown = float(np.min(p_v_extend/rolling_window_max-1)[0])
  vol_annual=np.std(X)*np.sqrt(250)
  ratio_annual=(arith_mean_ret_annual-rf)/vol_annual
  kurt_annual=kurtosis(X*250)
  skew_annual=skew(X*250)
  kurt_day=kurtosis(X)
  skew_day=skew(X)
  z=norm.ppf(1-confidenceLevel)
  t=z+((1/6)*(z**2-1)*skew_day)+((1/24)*(z**3-3*z))*kurt_day-((1/36)*(2*z**3-5*z)*(skew_day**2))
  mVaR= position*(np.mean(X)+t*np.std(X))*np.sqrt(250)
  alpha=norm.ppf(1-confidenceLevel, np.mean(X), np.std(X))
  VaR= position*(alpha)
  VaR_annual=VaR*np.sqrt(250)
  CVaR = position*np.mean(X[X≤np.quantile(X,1-confidenceLevel)])[0]*np.sqrt(250)
  df=pd.DataFrame([
                   cum_ret_annual,
                   arith_mean_ret_annual[0],
                   geom_mean_ret_annual,min_ret_annual[0],
                   ten_day_drawdown,vol_annual[0],
                   ratio_annual[0],
                   kurt_annual[0],
                   skew_annual[0],
                   mVaR[0],
                   VaR[0],
                   VaR_annual[0],
                   CVaR],
                  index=['Cumulative Returns (Annual)',
                         'Arithmetic Mean Returns (Annual)',
                         'Geometric Mean Returns (Annual)',
                         'Minimum Return (Annual)',
                         'Max 10-day Drawdown',
                         'Volatility',
                         'Sharpe Ratio (Annual)',
                         'Kurtosis (Annual)',
                         'Skew (Annual)',
                         'mVaR (Annual)',
                         'VaR (Daily)',
                         'VaR (Annual)',
                         'CVaR (Annual)'],
                  columns=['result'])
  return df

def omega(ret_r, factor_r, return_v, factor_v, lamb_, beta_tm_, wp_):
  rf = np.asarray(factor_r['RF'])
  rM_rf = np.asarray(factor_r['Mkt-RF'])
  rSMB = np.asarray(factor_r['SMB'])
  rHML = np.asarray(factor_r['HML'])
  SPY = np.asarray(ret_r['SPY'])
  ri = np.asarray(ret_r)
  var_market = np.var(SPY,ddof=1)
  beta_im = np.array([0.0]*12)
  for i in range (12):
    temp = np.cov(ri[:,i],SPY,ddof=1)
    beta_im[i] = temp[0,1] / var_market
  Ri = ri - rf.reshape(-1,1)
  f = np.array([rM_rf, rSMB, rHML])
  F = f.T
  lr = linear_model.LinearRegression().fit(F, Ri)
  alpha = lr.intercept_
  B = lr.coef_
  ft = f[:,-1]
  rho_r = alpha + B.dot(ft) + rf[-1]

  rf_v = np.asarray(factor_v['RF'])
  rM_rf_v = np.asarray(factor_v['Mkt-RF'])
  rSMB_v = np.asarray(factor_v['SMB'])
  rHML_v = np.asarray(factor_v['HML'])
  SPY_v = np.asarray(return_v['SPY'])
  ri_v = np.asarray(return_v)
  var_market_v = np.var(SPY_v,ddof=1)
  beta_im_v = np.array([0.0]*12)
  for i in range (12):
    temp_v = np.cov(ri_v[:,i],SPY_v,ddof=1)
    beta_im_v[i] = temp_v[0,1] / var_market_v
  Ri_v = ri_v - rf_v.reshape(-1,1)
  f_v = np.array([rM_rf_v, rSMB_v, rHML_v])
  F_v = f_v.T
  lr_v = linear_model.LinearRegression().fit(F_v, Ri_v)
  alpha_v = lr_v.intercept_
  B_v = lr_v.coef_
  eph_v = Ri_v.T - (alpha_v.reshape(-1,1) + B_v.dot(f_v))
  eph2_v = np.cov(eph_v,ddof=1)
  eph2_diag_v = np.diag(eph2_v)
  D_v = np.diag(eph2_diag_v)
  omega_f_v = np.cov(f_v,ddof=1)
  cov_Rt_v = B_v.dot(omega_f_v).dot(B_v.T) + D_v
  result = var_w(rho_r, lamb_, cov_Rt_v, wp_, beta_im_v ,beta_tm_)
  return result

"""## **Data Processing**

### Containers
"""

# Data containers:
p_u = pd.DataFrame()
p_aapl = pd.DataFrame()
p_spy = pd.DataFrame()

# Ticker containers:
u_tix = ['FXE', 'EWJ', 'GLD', 'QQQ', 'SPY', 'SHV', 'GAF', 'DBA', 'USO', 'XBI', 'ILF', 'EPP', 'FEZ']
aapl_tix = ['AAPL']
spy_tix = ['SPY']

"""### Load Data

#### Fama French Factors

A three-factor model proposed by Fama and French(1993), includes not only market excess return, but a capitalization size and book to market ratio will also be added in as influencing factors.

The random return of a given security is given by the formulas (equivalent),

\begin{equation}
\boxed{r = r_{f}+\beta_{1}(r_{m}-r_{f})+\beta{2}(SMB)+\beta_{3}(HML)+\epsilon}
\end{equation}

\begin{equation}
\boxed{R_{i}-r_{f}=\alpha_{i}+\beta{i}^{MKT}(R_{M}-r_{f})+\beta_{i}^{SMB}R_{SMB}+\beta_{i}^{HML}R_{HML}}
\end{equation}

- rSMB represents small size variables minus big one
- rHML represents high minus low in book value to equity to book value to the market.
"""

# Using definition above:
# Fama/French Three-Factor Model:
'''
ff_3_daily = fama_french('daily', 3)
print('Fama/French Three-Factor Model Daily Data\n' + str(ff_3_daily.tail(10)))

ff_3_weekly = fama_french('weekly', 3)
print('Fama/French Three-Factor Model Weekly Data\n' + str(ff_3_weekly.tail(10)))

ff_3_monthly = fama_french('monthly',3)
print('Fama/French Three-Factor Model Monthly Data\n' + str(ff_3_monthly.tail(10)))

ff_3_annual = fama_french('annual', 3)
print('Fama/French Three-Factor Model Annual Data\n' + str(ff_3_annual.tail(10)))
'''

# Manual loading:
# Fama/French Three-Factor Model:
ff_3_annual = pd.read_csv(data_dir + 'F-F_Research_Data_Factors.CSV', skiprows = 3, index_col = 0)
ff_3_annual = ff_3_annual.iloc[:-1]
ff_3_annual = ff_3_annual.iloc[1134:,]
ff_3_annual.index = ff_3_annual.index.map(lambda h: '  '.join(h).replace(' ', ''))
ff_3_annual.index = pd.to_datetime(ff_3_annual.index, format='%Y')
ff_3_annual = ff_3_annual.dropna()
#ff_3_annual = ff_3_annual/100
print('Fama/French Three-Factor Model Annual Data\n' + str(ff_3_annual.head(10)))

ff_3_monthly = pd.read_csv(data_dir + 'F-F_Research_Data_Factors.CSV', skiprows = 3, index_col = 0)
ff_3_monthly = ff_3_monthly.iloc[0:1131,]
ff_3_monthly = ff_3_monthly.dropna()
ff_3_monthly.index = pd.to_datetime(ff_3_monthly.index, format= '%Y%m')
#ff_3_monthly = ff_3_monthly/100
print('Fama/French Three-Factor Model Monthly Data\n' + str(ff_3_monthly.head(10)))

ff_3_weekly = pd.read_csv(data_dir + 'F-F_Research_Data_Factors_weekly.csv', skiprows = 3, index_col = 0)
ff_3_weekly = ff_3_weekly.dropna()
ff_3_weekly.index = pd.to_datetime(ff_3_weekly.index, format= '%Y%m%d')
ff_3_weekly = ff_3_weekly/100
print('Fama/French Three-Factor Model Weekly Data\n' + str(ff_3_weekly.head(10)) + '\n')

ff_3_daily = pd.read_csv(data_dir + 'F-F_Research_Data_Factors_daily.csv', skiprows = 3, index_col = 0)
ff_3_daily = ff_3_daily.dropna()
ff_3_daily.index = pd.to_datetime(ff_3_daily.index, format= '%Y%m%d')
ff_3_daily = ff_3_daily/100
print('Fama/French Three-Factor Model Daily Data\n' + str(ff_3_daily.head(10)) + '\n')

# Fama/French 5-Factor Model:
ff_5_annual = pd.read_csv(data_dir + 'F-F_Research_Data_5_Factors_2x3.csv', skiprows = 3, index_col = 0)
ff_5_annual = ff_5_annual.iloc[690:,]
ff_5_annual.index = ff_5_annual.index.map(lambda h: '  '.join(h).replace(' ', ''))
ff_5_annual.index = pd.to_datetime(ff_5_annual.index, format='%Y')
ff_5_annual = ff_5_annual.dropna()
#ff_5_annual = ff_5_annual/100
print('Fama/French 5-Factor Model Annual Data\n' + str(ff_5_annual.head(10)) + '\n')

ff_5_monthly = pd.read_csv(data_dir + 'F-F_Research_Data_5_Factors_2x3.csv', skiprows = 3, index_col = 0)
ff_5_monthly = ff_5_monthly.iloc[:688,]
ff_5_monthly = ff_5_monthly.dropna()
ff_5_monthly.index = pd.to_datetime(ff_5_monthly.index, format='%Y%m')
#ff_5_monthly = ff_5_monthly/100
print('Fama/French 5-Factor Model Monthly Data\n' + str(ff_5_monthly.tail(10)) + '\n')

ff_5_daily = pd.read_csv(data_dir + 'F-F_Research_Data_5_Factors_2x3_daily.csv', skiprows = 3, index_col = 0)
ff_5_daily = ff_5_daily.dropna()
ff_5_daily.index = pd.to_datetime(ff_5_daily.index, format='%Y%m%d')
ff_5_daily = ff_5_daily/100
print('Fama/French 5-Factor Model Daily Data\n' + str(ff_5_daily.head(10)) + '\n')

# Last date of time series data must match that of the Fama/French data:
last_datapoint = str(ff_3_daily.index[-1].strftime('%m/%d/%Y'))
print('Last Date for Fama/French data: ' + last_datapoint)

"""#### Historical Time Series

The following ETFs represent the investment Universe of our portfolios. They range from the S&P500 to ETFs representing all continents such as Europe, Asia and Africa and asset types such as bonds, stocks, and commodities.

1. CurrencyShares Euro Trust (FXE)
2. iShares MSCI Japan Index (EWJ)
3. SPDR GOLD Trust (GLD)
4. Powershares NASDAQ-100 Trust (QQQ)
5. SPDR S&P500 (SPY) **(THE MARKET PORTFOLIO S&P500 IS THE BENCHMARK)**
6. iShares Lehman Short Treasury Bond (SHV)
7. PowerShares DB Agriculture Fund (DBA)
8. United States Oil Fund LP (USO)
9. SPDR S&P Biotech (XBI)
10. iShares S&P Latin America 40 Index (ILF)
11. iShares MSCI Pacific ex-Japan Index Fund (EPP)
12. SPDR DJ Euro Stoxx 50 (FEZ)

From this universe, we have created portfolios by utilizing the Three-Factor Fama-French model. The investment portfolio that we created is compared to the following benchmark portfolios:

1.	The Market Portfolio (S&P500)

The dataset includes daily price data between March 1st, 2007 to October 31th, 2020. We choose this investment horizon to match the Fama-French Factor data available.

We have used three different look-back periods, which we have defined as: A. Short Term – 60 Days B. Medium Term – 120 Days C. Long Term – 200 Days To calculate the risk-return parameters of then portfolio we have used the target Beta as -1, -0.5, 0, 0.5, 1 and 1.5. The rebalance period is kept as one week as specified in the project.
"""

# Retrieve ETF Data:
start_date = '07/24/2007'
end_date = '10/30/2020'

for i in u_tix:
    tmp = web.DataReader(i, 'yahoo', start_date, end_date)
    p_u[i] = tmp['Adj Close']
for i in aapl_tix:
    tmp = web.DataReader(i, 'yahoo', start_date, end_date)
    p_aapl[i] = tmp['Adj Close']
for i in spy_tix:
    tmp = web.DataReader(i, 'yahoo', start_date, end_date)
    p_spy[i] = tmp['Adj Close']

"""### Preprocess Data"""

# Clean data:
p_u.isnull().sum().sum()
p_aapl.isnull().sum().sum()
p_spy.isnull().sum().sum()

p_u = p_u.dropna()
p_aapl = p_aapl.dropna()
p_spy = p_spy.dropna()

# Sample data:
sys.stdout.write('\nInvestment Universe U Time Series ({} - {}):\n'.format(start_date, end_date) + str(p_u.tail(10)))
print('\n' + str(p_u.shape))
sys.stdout.write('\nAAPL Time Series ({} - {}):\n'.format(start_date, end_date) + str(p_aapl.head(10)))
print('\n' + str(p_aapl.shape))
sys.stdout.write('\nSPY Time Series ({} - {}):\n'.format(start_date, end_date) + str(p_spy.head(10)))
print('\n' + str(p_spy.shape))

"""### Transform Data"""

# Useful variables:
w = weights_randn(len(u_tix))
print('Random weights:\n' + str(w))
days = 252

# ETF Returns:
R_u = p_u.pct_change()
R_u = R_u.dropna()
print('\nETF Portfolio Prices:\n' + str(R_u.head(10)))
rho_u = np.mean(R_u, axis=0)
print('\nETF Portfolio Mean Returns:\n' + str(rho_u))

# SPY Returns:
R_spy = p_spy.pct_change()
R_spy = R_spy.dropna()
print('\nSPY Prices:\n' + str(R_spy.head(10)))
rho_spy = np.mean(R_spy, axis=0)
print('\nSPY Mean Return:\n' + str(rho_spy))

"""### Save Data"""

# Securities:
save_data(p_u, 'p_u')
save_data(p_aapl, 'p_aapl')
save_data(p_spy, 'p_spy')
save_data(R_spy,'R_spy')
save_data(R_u,'R_u')

# Fama-French Processed Datasets (Archive):
save_data(ff_3_annual, 'ff_3_annual')
save_data(ff_5_annual, 'ff_5_annual')
save_data(ff_3_monthly, 'ff_3_monthly')
save_data(ff_5_monthly, 'ff_5_monthly')
save_data(ff_3_weekly, 'ff_3_weekly')
save_data(ff_3_daily, 'ff_3_daily')
save_data(ff_5_daily, 'ff_5_daily')

"""### Visualize Data"""

# Visualize ETF Log-Returns:
R_spy = np.log(p_spy / p_spy.shift(1))
returns_u, axs = plt.subplots(4,3,figsize=(15, 7.5))
returns_u.suptitle('Log-Returns of Portfolio Securities', fontweight='bold', fontsize=15)
axs[0,0].plot(R_u['FXE'], 'black', linewidth=0.5, alpha=0.9)
axs[0,0].set_title('FXE')
axs[0,1].plot(R_u['EWJ'], 'black', linewidth=0.5, alpha=0.9)
axs[0,1].set_title('EWJ')
axs[0,2].plot(R_u['GLD'], 'black', linewidth=0.5, alpha=0.9)
axs[0,2].set_title('GLD')
axs[1,0].plot(R_u['QQQ'], 'black', linewidth=0.5, alpha=0.9)
axs[1,0].set_title('QQQ')
axs[1,1].plot(R_u['SPY'], 'black', linewidth=0.5, alpha=0.9)
axs[1,1].set_title('SPY')
axs[1,2].plot(R_u['SHV'], 'black', linewidth=0.5, alpha=0.9)
axs[1,2].set_title('SHV')
axs[2,0].plot(R_u['DBA'], 'black', linewidth=0.5, alpha=0.9)
axs[2,0].set_title('DBA')
axs[2,1].plot(R_u['USO'], 'black', linewidth=0.5, alpha=0.9)
axs[2,1].set_title('USO')
axs[2,2].plot(R_u['XBI'], 'black', linewidth=0.5, alpha=0.9)
axs[2,2].set_title('XBI')
axs[3,0].plot(R_u['ILF'], 'black', linewidth=0.5, alpha=0.9)
axs[3,0].set_title('ILF')
axs[3,1].plot(R_u['EPP'], 'black', linewidth=0.5, alpha=0.9)
axs[3,1].set_title('EPP')
axs[3,2].plot(R_u['FEZ'], 'black', linewidth=0.5, alpha=0.9)
axs[3,2].set_title('FEZ')
plt.tight_layout()
returns_u.subplots_adjust(top=0.9)
plt.savefig(graphs_dir + 'returns_u_log.png', bbox_inches='tight')

# Visualize ETF Price Time Series:
R_u = p_u
returns_u, axs = plt.subplots(4,3,figsize=(15, 7.5))
returns_u.suptitle('Historical Time Series of Portfolio Securities', fontweight='bold', fontsize=15)
axs[0,0].plot(R_u['FXE'], 'black', linewidth=0.5, alpha=0.9)
axs[0,0].set_title('FXE')
axs[0,1].plot(R_u['EWJ'], 'black', linewidth=0.5, alpha=0.9)
axs[0,1].set_title('EWJ')
axs[0,2].plot(R_u['GLD'], 'black', linewidth=0.5, alpha=0.9)
axs[0,2].set_title('GLD')
axs[1,0].plot(R_u['QQQ'], 'black', linewidth=0.5, alpha=0.9)
axs[1,0].set_title('QQQ')
axs[1,1].plot(R_u['SPY'], 'black', linewidth=0.5, alpha=0.9)
axs[1,1].set_title('SPY')
axs[1,2].plot(R_u['SHV'], 'black', linewidth=0.5, alpha=0.9)
axs[1,2].set_title('SHV')
axs[2,0].plot(R_u['DBA'], 'black', linewidth=0.5, alpha=0.9)
axs[2,0].set_title('DBA')
axs[2,1].plot(R_u['USO'], 'black', linewidth=0.5, alpha=0.9)
axs[2,1].set_title('USO')
axs[2,2].plot(R_u['XBI'], 'black', linewidth=0.5, alpha=0.9)
axs[2,2].set_title('XBI')
axs[3,0].plot(R_u['ILF'], 'black', linewidth=0.5, alpha=0.9)
axs[3,0].set_title('ILF')
axs[3,1].plot(R_u['EPP'], 'black', linewidth=0.5, alpha=0.9)
axs[3,1].set_title('EPP')
axs[3,2].plot(R_u['FEZ'], 'black', linewidth=0.5, alpha=0.9)
axs[3,2].set_title('FEZ')
plt.tight_layout()
returns_u.subplots_adjust(top=0.9)
plt.savefig(graphs_dir + 'prices_u_raw.png', bbox_inches='tight')

# Visualize ETF Price Time Series:
fig = plt.figure(figsize=(15, 7.5))
ts_u = fig.add_subplot(111)
ts_u.plot(R_u['FXE'], linewidth=0.5, alpha=0.9, label='FXE')
ts_u.plot(R_u['EWJ'], linewidth=0.5, alpha=0.9, label='EWJ')
ts_u.plot(R_u['GLD'], linewidth=0.5, alpha=0.9, label='GLD')
ts_u.plot(R_u['QQQ'], linewidth=0.5, alpha=0.9, label='QQQ')
ts_u.plot(R_u['SPY'], linewidth=0.5, alpha=0.9, label='SPY')
ts_u.plot(R_u['SHV'], linewidth=0.5, alpha=0.9, label='SHV')
ts_u.plot(R_u['DBA'], linewidth=0.5, alpha=0.9, label='DBA')
ts_u.plot(R_u['USO'], linewidth=0.5, alpha=0.9, label='USO')
ts_u.plot(R_u['XBI'], linewidth=0.5, alpha=0.9, label='XBI')
ts_u.plot(R_u['ILF'], linewidth=0.5, alpha=0.9, label='ILF')
ts_u.plot(R_u['EPP'], linewidth=0.5, alpha=0.9, label='EPP')
ts_u.plot(R_u['FEZ'], linewidth=0.5, alpha=0.9, label='FEZ')
ts_u.set_xlabel('Year', fontweight='bold', fontsize=12)
ts_u.set_ylabel('Price', fontweight='bold', fontsize=12)
ts_u.set_title('Historical Time Series of Portfolio Securities', fontweight='bold', fontsize=15)
ts_u.legend(loc='upper right', fontsize=10)
plt.savefig(graphs_dir + 'rho_u.png', bbox_inches='tight')

"""## **Analysis**

### Theory

#### Optimization Problem:

The strategy aims to maximize return with a certain Target Beta under constraints.

It is defined as,

\begin{cases}
\max\limits_{{\omega ∈ ℝ^{n}}}\rho^{T}\omega-\lambda(\omega-\omega_{p})^{T}\Sigma(\omega-\omega_{p})\\
\sum_{i=1}^{n} \beta_{i}^{m}\omega_{i}=\beta_{T}^{m}\\
\sum_{i=1}^{n} \omega_{i}=1, -2\leq\omega_{i}\leq2
\end{cases}

$\Sigma$ is the the covariance matrix between the securities returns (computed from
the Factor Model), $\omega_{p}$ is the composition of a reference Portfolio (the previous Portfolio when rebalancing the portfolio and $\omega_{p}$ has all its components equal to $1/n$ for the first allocation) and $\lambda$ is a small regularization parameter to limit the turnover;

$\beta_{i}^{m}=\frac{cov(r_{i},r_{M}}{\sigma^{2}(r_{M})}$ is the Beta of security $S_{i}$ as defined in the CAPM Model so that $\beta_{P}^{m}=\sum_{i=1}^{n}\beta_{i}^{m}\omega_{i}$ is the Beta of the Portfolio;

$\beta_{T}^{m}$ is the Portfolio's Target Beta, for example $\beta_{T}^{m}=-1$, $\beta_{T}^{m}=-0.5$, $\beta_{T}^{m}=0$, $\beta_{T}^{m}=0.5$, $\beta_{T}^{m}=1.5$.

#### Equivalent Optimization Problem:

We can reformulate the optimization problem above to make the programming process more straightforward:

$(\omega-\omega_{p})^{T}\Sigma(\omega-\omega_{p})\rightarrow$

$=(\omega-\omega_{p})^{T}\Sigma\omega-(\omega-\omega_{p} )^{T}\Sigma\omega_{p}$

$=\omega^{T} \Sigma\omega-2(\omega^{T} \Sigma\omega_{p})+\omega_{p}^{T}\Sigma \omega_{p}$

We simplify,
- $d=\rho-2\lambda\Sigma\omega_{p}$
- $P=\lambda\Sigma$

Finally,

$\max\limits_{{\omega ∈ ℝ^{n}}}(\rho-2\lambda\Sigma\omega_{p} )^{T} \omega-\lambda\omega^{T}\Sigma\omega+\lambda\omega_{p}^{T}\Sigma\omega_{p}=\max\limits_{{\omega ∈ ℝ^{n}}}d^{T}\omega-\omega^{T}P\omega$

---

The following formulation is equivalent,

\begin{cases}
\max\limits_{{\omega ∈ ℝ^{n}}}d^{T}\omega-\omega^{T}P\omega\\
\sum_{i=1}^{n} \beta_{i}^{m}\omega_{i}=\beta_{T}^{m}\\
\sum_{i=1}^{n} \omega_{i}=1, -2\leq\omega_{i}\leq2
\end{cases}
- $\Sigma$ is the the covariance matrix between the returns of the portfolio assets;
- $\omega_{p}$ is the composition of a reference Portfolio:
  - When rebalancing the portfolio, $\omega_{p}$ is the previous portfolio
  - $\omega_{p}$ has all its components equal to $1/n$ for the first allocation
- $\lambda$ is a regularization parameter to limit the turnover
- $\beta_{i}^{m}=\frac{cov(r_{i},r_{M}}{\sigma^{2}(r_{M})}$ is the Beta of security $S_{i}$ as defined in the CAPM Model s.t. $\beta_{P}^{m}=\sum_{i=1}^{n}\beta_{i}^{m}\omega_{i}$ is the portfolio Beta
- $\beta_{T}^{m}$ is the Portfolio's Target Beta.

### Algebra
"""

# Create hybrid dataset:
R_u_ff = pd.merge(R_u,ff_3_daily,how='inner',left_index=True,right_index=True)
R_spy_ff = pd.merge(R_spy,ff_3_daily,how='inner',left_index=True,right_index=True)

# Rename Market Excess Column Index:
R_u_ff.rename(columns={'Mkt-RF':'Mkt_RF'}, inplace=True)
R_u_ff['Portfolio_Excess'] = R_u_ff.sum(axis=1) - R_u_ff['RF']
print(R_u_ff.head(10))

# Quick save:
save_data(R_u_ff, 'R_u_ff')
save_data(R_spy_ff, 'R_spy_ff')

# Estimate Security Betas:
betas = []
for i in range(0,len(u_tix)):
  reg_mult = smf.formula.ols(formula = 'R_u_ff.iloc[:, i] - RF ~ Mkt_RF - RF + SMB + HML', data = R_u_ff).fit()
  betas.append(list(reg_mult.params))

betas = pd.DataFrame(betas, index=u_tix)
betas.columns = ['Intercept', 'Mkt_RF', 'SMB', 'HML']
print(betas)

# Quick save:
save_data(betas, 'betas')

# Calculate Annualized Average Expected Returns under FF Three-Factor Model:
rho_daily = []

for i in range(0,len(u_tix)):
  step_0 = (R_spy_ff.sum(axis=1) - R_spy_ff['RF']).mul((betas.iloc[i,0] + betas.iloc[i,1]))
  step_1 = R_spy_ff['SMB'].mul(betas.iloc[i,2])
  step_2 = R_spy_ff['HML'].mul(betas.iloc[i,3])
  step_4 = step_0 + step_1 + step_2
  rho_daily.append(step_4)
rho_daily = pd.DataFrame(rho_daily)
rho_daily = rho_daily.T
rho_daily.columns = u_tix
print('Daily Average Expected Returns:\n' + str(rho_daily.head(10)))

rho_annual = rho_daily * 252
print('Annualized Average Expected Returns:\n' + str(rho_annual.head(10)))

# Quick Save:
save_data(rho_daily, 'rho_daily')
save_data(rho_annual, 'rho_annual')

# Calculate other variables:
ones = np.ones(len(u_tix))
mu_u = np.mean(rho_annual, axis=0)
print('Mean Average Expected Returns (Annual): \n' + str(mu_u))
mu_u_daily = np.mean(rho_daily, axis=0)
print('\nMean Average Expected Returns (Daily): \n' + str(mu_u_daily))
Sigma_u = np.cov(rho_annual, rowvar=False)
print('\nCovariance Matrix: \n' + str(Sigma_u))

P = 2 * (Sigma_u + 0.01 * np.identity(len(mu_u)))
#print('\nP Matrix: \n' + str(P))
omega_u = np.repeat(1/len(mu_u), len(mu_u))
A_eq = np.repeat(1,len(mu_u))
A_mat = pd.DataFrame(np.identity(len(mu_u))).merge(pd.DataFrame(-np.identity(len(mu_u))))

"""### Visualizations"""

# Visualize Daily Average Expected Returns:
r = np.transpose(np.linspace(0, 1, len(rho_annual)))
exp_returns_day, axs = plt.subplots(4,3,figsize=(15, 7.5))
exp_returns_day.suptitle('Daily Average Expected Returns of Portfolio Securities', fontweight='bold', fontsize=15)
axs[0,0].plot(rho_daily['FXE'], 'black', linewidth=0.5, alpha=0.9)
axs[0,0].set_title('FXE')
axs[0,1].plot(rho_daily['EWJ'], 'black', linewidth=0.5, alpha=0.9)
axs[0,1].set_title('EWJ')
axs[0,2].plot(rho_daily['GLD'], 'black', linewidth=0.5, alpha=0.9)
axs[0,2].set_title('GLD')
axs[1,0].plot(rho_daily['QQQ'], 'black', linewidth=0.5, alpha=0.9)
axs[1,0].set_title('QQQ')
axs[1,1].plot(rho_daily['SPY'], 'black', linewidth=0.5, alpha=0.9)
axs[1,1].set_title('SPY')
axs[1,2].plot(rho_daily['SHV'], 'black', linewidth=0.5, alpha=0.9)
axs[1,2].set_title('SHV')
axs[2,0].plot(rho_daily['DBA'], 'black', linewidth=0.5, alpha=0.9)
axs[2,0].set_title('DBA')
axs[2,1].plot(rho_daily['USO'], 'black', linewidth=0.5, alpha=0.9)
axs[2,1].set_title('USO')
axs[2,2].plot(rho_daily['XBI'], 'black', linewidth=0.5, alpha=0.9)
axs[2,2].set_title('XBI')
axs[3,0].plot(rho_daily['ILF'], 'black', linewidth=0.5, alpha=0.9)
axs[3,0].set_title('ILF')
axs[3,1].plot(rho_daily['EPP'], 'black', linewidth=0.5, alpha=0.9)
axs[3,1].set_title('EPP')
axs[3,2].plot(rho_daily['FEZ'], 'black', linewidth=0.5, alpha=0.9)
axs[3,2].set_title('FEZ')
plt.tight_layout()
exp_returns_day.subplots_adjust(top=0.9)
plt.savefig(graphs_dir + 'exp_returns_daily.png', bbox_inches='tight')

# Visualize Annualized Average Expected Returns:
r = np.transpose(np.linspace(0, 1, len(rho_annual)))
exp_returns_yr, axs = plt.subplots(4,3,figsize=(15, 7.5))
exp_returns_yr.suptitle('Annualized Average Expected Returns of Portfolio Securities', fontweight='bold', fontsize=15)
axs[0,0].plot(r,rho_annual['FXE'], 'black', linewidth=0.5, alpha=0.9)
axs[0,0].set_title('FXE')
axs[0,1].plot(r,rho_annual['EWJ'], 'black', linewidth=0.5, alpha=0.9)
axs[0,1].set_title('EWJ')
axs[0,2].plot(r,rho_annual['GLD'], 'black', linewidth=0.5, alpha=0.9)
axs[0,2].set_title('GLD')
axs[1,0].plot(r,rho_annual['QQQ'], 'black', linewidth=0.5, alpha=0.9)
axs[1,0].set_title('QQQ')
axs[1,1].plot(r,rho_annual['SPY'], 'black', linewidth=0.5, alpha=0.9)
axs[1,1].set_title('SPY')
axs[1,2].plot(r,rho_annual['SHV'], 'black', linewidth=0.5, alpha=0.9)
axs[1,2].set_title('SHV')
axs[2,0].plot(r,rho_annual['DBA'], 'black', linewidth=0.5, alpha=0.9)
axs[2,0].set_title('DBA')
axs[2,1].plot(r,rho_annual['USO'], 'black', linewidth=0.5, alpha=0.9)
axs[2,1].set_title('USO')
axs[2,2].plot(r,rho_annual['XBI'], 'black', linewidth=0.5, alpha=0.9)
axs[2,2].set_title('XBI')
axs[3,0].plot(r,rho_annual['ILF'], 'black', linewidth=0.5, alpha=0.9)
axs[3,0].set_title('ILF')
axs[3,1].plot(r,rho_annual['EPP'], 'black', linewidth=0.5, alpha=0.9)
axs[3,1].set_title('EPP')
axs[3,2].plot(r,rho_annual['FEZ'], 'black', linewidth=0.5, alpha=0.9)
axs[3,2].set_title('FEZ')
plt.tight_layout()
exp_returns_yr.subplots_adjust(top=0.9)
plt.savefig(graphs_dir + 'exp_returns_annual.png', bbox_inches='tight')

# Visualize Daily Average Expected Returns (Superimposed):
#fig = plt.figure(figsize=(15, 7.5))
#rho_u_day = fig.add_subplot(111)
#rho_u_day.plot(rho_daily['FXE'], linewidth=0.5, alpha=0.9, label='FXE')
#rho_u_day.plot(rho_daily['EWJ'], linewidth=0.5, alpha=0.9, label='EWJ')
#rho_u_day.plot(rho_daily['GLD'], linewidth=0.5, alpha=0.9, label='GLD')
#rho_u_day.plot(rho_daily['QQQ'], linewidth=0.5, alpha=0.9, label='QQQ')
#rho_u_day.plot(rho_daily['SPY'], linewidth=0.5, alpha=0.9, label='SPY')
#rho_u_day.plot(rho_daily['SHV'], linewidth=0.5, alpha=0.9, label='SHV')
#rho_u_day.plot(rho_daily['DBA'], linewidth=0.5, alpha=0.9, label='DBA')
#rho_u_day.plot(rho_daily['USO'], linewidth=0.5, alpha=0.9, label='USO')
#rho_u_day.plot(rho_daily['XBI'], linewidth=0.5, alpha=0.9, label='XBI')
#rho_u_day.plot(rho_daily['ILF'], linewidth=0.5, alpha=0.9, label='ILF')
#rho_u_day.plot(rho_daily['EPP'], linewidth=0.5, alpha=0.9, label='EPP')
#rho_u_day.plot(rho_daily['FEZ'], linewidth=0.5, alpha=0.9, label='FEZ')
#rho_u_day.set_xlabel('Year', fontweight='bold', fontsize=12)
#rho_u_day.set_ylabel('Return', fontweight='bold', fontsize=12)
#rho_u_day.set_title('Daily Average Expected Returns of Portfolio Securities', fontweight='bold', fontsize=15)
#rho_u_day.legend(loc='upper right', fontsize=10)
#plt.savefig(graphs_dir + 'exp_returns_daily_all.png', bbox_inches='tight')

# Visualize Annualized Average Expected Returns (Superimposed):
#fig = plt.figure(figsize=(15, 7.5))
#rho_u_yr = fig.add_subplot(111)
#rho_u_yr.plot(rho_annual['FXE'], linewidth=0.5, alpha=0.9, label='FXE')
#rho_u_yr.plot(rho_annual['EWJ'], linewidth=0.5, alpha=0.9, label='EWJ')
#rho_u_yr.plot(rho_annual['GLD'], linewidth=0.5, alpha=0.9, label='GLD')
#rho_u_yr.plot(rho_annual['QQQ'], linewidth=0.5, alpha=0.9, label='QQQ')
#rho_u_yr.plot(rho_annual['SPY'], linewidth=0.5, alpha=0.9, label='SPY')
#rho_u_yr.plot(rho_annual['SHV'], linewidth=0.5, alpha=0.9, label='SHV')
#rho_u_yr.plot(rho_annual['DBA'], linewidth=0.5, alpha=0.9, label='DBA')
#rho_u_yr.plot(rho_annual['USO'], linewidth=0.5, alpha=0.9, label='USO')
#rho_u_yr.plot(rho_annual['XBI'], linewidth=0.5, alpha=0.9, label='XBI')
#rho_u_yr.plot(rho_annual['ILF'], linewidth=0.5, alpha=0.9, label='ILF')
#rho_u_yr.plot(rho_annual['EPP'], linewidth=0.5, alpha=0.9, label='EPP')
#rho_u_yr.plot(rho_annual['FEZ'], linewidth=0.5, alpha=0.9, label='FEZ')
#rho_u_yr.set_xlabel('Year', fontweight='bold', fontsize=12)
#rho_u_yr.set_ylabel('Return', fontweight='bold', fontsize=12)
#rho_u_yr.set_title('Annualized Average Expected Returns of Portfolio Securities', fontweight='bold', fontsize=15)
#rho_u_yr.legend(loc='upper right', fontsize=10)
#plt.savefig(graphs_dir + 'exp_returns_annual_all.png', bbox_inches='tight')

"""## **Deployment**

### Load Data
"""

etf_ff = pd.read_csv(data_dir + '2020.12.22_ETF_FF.csv',index_col=0).fillna(0)
ff_3_daily = etf_ff.loc[:,'Mkt-RF':r'RF']
etf = etf_ff.loc[:,:'FEZ']
etf = etf.drop(['GAF'], axis=1)
R_etf = (etf/etf.shift(1)-1).replace(np.nan, 0)
R_etf = R_etf.loc['3/22/2007':'10/30/2020',:]
ff_3_daily = ff_3_daily.loc['3/22/2007':'10/30/2020',:].replace(np.nan, 0)

# ETF Data Inspection:
R_etf

# Fama French Factor Data Inspection:
ff_3_daily

"""### Before the Subprime Crisis
Period: March 22, 2007 - March 3, 2008

#### Tearsheet
"""

# Pre-Subprime Crisis:
pre_subprime_R_u = R_etf.loc[:'3/3/2008',:'FEZ']
pre_subprime_ff_factors = ff_3_daily.loc[:'3/3/2008','Mkt-RF':'RF']
pre_subprime_lookbacks = [[60,60], [60,120], [90,60], [90,120], [120,60], [120,120]]
pre_subprime_betas = [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0]
pre_subprime_exec = pd.DataFrame([])
pre_subprime_final = pd.DataFrame([])
omegas = []
for lb in pre_subprime_lookbacks:
    for bt in pre_subprime_betas:
        res = backtesting(pre_subprime_R_u,
                        pre_subprime_ff_factors,
                        return_period = lb[0],
                        variance_period = lb[1],
                        lamb = 10,
                        beta_tm = bt)
        omegas.append(res[1])
        res = pd.DataFrame(res[0],index = pd.to_datetime(pre_subprime_R_u.index))
        res_perf = analytics(X = res,rf = 0.06, confidenceLevel = 0.95, position = 100)
        pre_subprime_final = pd.concat([pre_subprime_final,res],axis = 1)
        pre_subprime_exec = pd.concat([pre_subprime_exec,res_perf],axis = 1)

pre_subprime_final = pd.concat([pre_subprime_final,pre_subprime_R_u['SPY']],axis = 1)

pre_subprime_spy_performance = analytics(X = pd.DataFrame(pre_subprime_R_u.loc[:,'SPY']),rf = 0.06, confidenceLevel = 0.95, position = 100)
pre_subprime_exec = pd.concat([pre_subprime_exec,pre_subprime_spy_performance],axis = 1)
pre_subprime_exec.columns = [['$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$',
                                  '$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$',
                                  '$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$',
                                  '$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$',
                                  '$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$',
                                  '$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','SPY'],
                                 ['\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0','']]
pre_subprime_final.columns = pre_subprime_exec.columns
save_data(pre_subprime_exec, 'pre_subprime_exec')

"""#### Plots

##### Notes

An investment strategy is abbreviated as $S_{Cov}^{E[r]}(\beta)$. In this implementation, *post_subprime_lookbacks* contains pairs $[E[r], Cov]$. *post_subprime_betas* contains the various target $\beta$.

##### Total Value
"""

# Total Value:
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111)
for i in range(36):
  ax.plot(100*(np.cumprod(pre_subprime_final.iloc[:,i]+1)),label = pre_subprime_final.columns[i][0]+', '+pre_subprime_final.columns[i][1])
  ax.legend(loc='best', ncol=4, fontsize=10)
plt.xlabel('t', fontweight='bold', fontsize=15)
plt.ylabel('Value', fontweight='bold', fontsize=15)
plt.title('Value of Investment Strategies During the Subprime Crisis', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '0_pre_subprime_strategy_val.png', bbox_inches='tight')

"""##### $S_{60}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_subprime_final.iloc[:,i]
  col_name = pre_subprime_final.columns[i]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(loc='left', fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.9).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 50, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{60}$ Returns Before the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '1_pre_subprime_ret_distS6060.png', bbox_inches='tight')

"""##### $S_{120}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_subprime_final.iloc[:,i+6]
  col_name = pre_subprime_final.columns[i+6]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(loc='left', fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.9).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 50, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('r', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{60}$ Returns Before the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '2_pre_subprime_ret_distS60120.png', bbox_inches='tight')

"""##### $S_{60}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_subprime_final.iloc[:,i+12]
  col_name = pre_subprime_final.columns[i+12]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(loc='left', fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.9).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 50, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{90}$ Returns Before the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '3_pre_subprime_ret_distS9060.png', bbox_inches='tight')

"""##### $S_{120}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_subprime_final.iloc[:,i+18]
  col_name = pre_subprime_final.columns[i+18]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(loc='left', fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.9).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 50, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{90}$ Returns Before the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '4_pre_subprime_ret_distS90120.png', bbox_inches='tight')

"""##### $S_{60}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_subprime_final.iloc[:,i+24]
  col_name = pre_subprime_final.columns[i+24]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(loc='left', fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.9).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 50, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{120}$ Returns Before the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '5_pre_subprime_ret_distS12060.png', bbox_inches='tight')

"""##### $S_{120}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_subprime_final.iloc[:,i+30]
  col_name = pre_subprime_final.columns[i+30]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(loc='left', fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.9).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 50, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{120}$ Returns Before the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '6_pre_subprime_ret_distS120120.png', bbox_inches='tight')

"""### During the Subprime Crisis
Period: March 3, 2008 - September 1, 2010

#### Tearsheet
"""

# During  the Subprime Crisis:
on_subprime_R_u = R_etf.loc['3/3/2008':'9/1/2010',:'FEZ']
on_subprime_ff_factors = ff_3_daily.loc['3/3/2008':'9/1/2010','Mkt-RF':'RF']
on_subprime_lookbacks = [[60,60], [60,120], [90,60], [90,120], [120,60], [120,120]]
on_subprime_betas = [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0]
on_subprime_exec = pd.DataFrame([])
on_subprime_final = pd.DataFrame([])
omegas = []
for lb in on_subprime_lookbacks:
  for bt in on_subprime_betas:
    res = backtesting(on_subprime_R_u,
                    on_subprime_ff_factors,
                    return_period = lb[0],
                    variance_period = lb[1],
                    lamb = 10,
                    beta_tm = bt)
    omegas.append(res[1])
    res = pd.DataFrame(res[0],index = pd.to_datetime(on_subprime_R_u.index))
    res_perf = analytics(X = res,rf = 0.06, confidenceLevel = 0.95, position = 100)
    on_subprime_final = pd.concat([on_subprime_final,res],axis = 1)
    on_subprime_exec = pd.concat([on_subprime_exec,res_perf],axis = 1)

on_subprime_final = pd.concat([on_subprime_final,on_subprime_R_u['SPY']],axis = 1)

on_subprime_spy_performance = analytics(X = pd.DataFrame(on_subprime_R_u.loc[:,'SPY']),rf = 0.06, confidenceLevel = 0.95, position = 100)
on_subprime_exec = pd.concat([on_subprime_exec,on_subprime_spy_performance],axis = 1)
on_subprime_exec.columns = [['$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$',
                                  '$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$',
                                  '$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$',
                                  '$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$',
                                  '$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$',
                                  '$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','SPY'],
                                 ['\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0','']]

on_subprime_final.columns = on_subprime_exec.columns
save_data(on_subprime_exec, 'on_subprime_exec')

"""#### Plots

##### Notes

An investment strategy is abbreviated as $S_{Cov}^{E[r]}(\beta)$. In this implementation, *post_subprime_lookbacks* contains pairs $[E[r], Cov]$. *post_subprime_betas* contains the various target $\beta$.

##### Total Value
"""

# During the Subprime crisis, plots:
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111)
for i in range(36):
  ax.plot(100*(np.cumprod(on_subprime_final.iloc[:,i]+1)),label = on_subprime_final.columns[i][0]+', '+on_subprime_final.columns[i][1])
  ax.legend(loc='best', ncol=4, fontsize=10)
plt.xlabel('t', fontweight='bold', fontsize=15)
plt.ylabel('Value', fontweight='bold', fontsize=15)
plt.title('Value of Investment Strategies During the Subprime Crisis', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '7_on_subprime_strategy_val.png', bbox_inches='tight')

"""##### $S_{60}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = on_subprime_final.iloc[:,i]
  col_name = on_subprime_final.columns[i]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x, y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:]) / 2
  cs = [c] * len(x)
  ax.bar(y, x, zs = z, zdir = 'y', color = cs, alpha = 0.7, width = 0.003, label = col_name[0]+', ' + col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 50, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{60}$ Returns During the Subprime Crisis', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '8_on_subprime_ret_distS6060.png', bbox_inches='tight')

"""##### $S_{120}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = on_subprime_final.iloc[:,i+6]
  col_name = on_subprime_final.columns[i+6]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/40,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('r', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{60}$ Returns During the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '9_on_subprime_ret_distS60120.png', bbox_inches='tight')

"""##### $S_{60}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = on_subprime_final.iloc[:,i+12]
  col_name = on_subprime_final.columns[i+12]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/50,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{90}$ Returns During the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '10_on_subprime_ret_distS9060.png', bbox_inches='tight')

"""##### $S_{120}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = on_subprime_final.iloc[:,i+18]
  col_name = on_subprime_final.columns[i+18]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/50,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{90}$ Returns During the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '11_on_subprime_ret_distS90120.png', bbox_inches='tight')

"""##### $S_{60}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = on_subprime_final.iloc[:,i+24]
  col_name = on_subprime_final.columns[i+24]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/50,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{120}$ Returns During the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '12_on_subprime_ret_distS12060.png', bbox_inches='tight')

"""##### $S_{120}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = on_subprime_final.iloc[:,i+30]
  col_name = on_subprime_final.columns[i+30]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/50,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{120}$ Returns During the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '13_on_subprime_ret_distS120120.png', bbox_inches='tight')

"""### After the Subprime Crisis
Period: September 1, 2010 - January 2, 2015

#### Tearsheet
"""

# After the Subprime Crisis:
post_subprime_R_u = R_etf.loc['9/1/2010':'1/2/2015',:'FEZ']
post_subprime_ff_factors = ff_3_daily.loc['9/1/2010':'1/2/2015','Mkt-RF':'RF']
post_subprime_lookbacks = [[60,60], [60,120], [90,60], [90,120], [120,60], [120,120]]
post_subprime_betas = [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0]
post_subprime_exec = pd.DataFrame([])
post_subprime_final = pd.DataFrame([])
omegas = []
for lb in post_subprime_lookbacks:
  for bt in post_subprime_betas:
    res = backtesting(post_subprime_R_u,
                    post_subprime_ff_factors,
                    return_period = lb[0],
                    variance_period = lb[1],
                    lamb = 10,
                    beta_tm = bt)
    omegas.append(res[1])
    res = pd.DataFrame(res[0],index = pd.to_datetime(post_subprime_R_u.index))
    res_perf = analytics(X = res,rf = 0.06, confidenceLevel = 0.95, position = 100)
    post_subprime_final = pd.concat([post_subprime_final,res],axis = 1)
    post_subprime_exec = pd.concat([post_subprime_exec,res_perf],axis = 1)

post_subprime_final = pd.concat([post_subprime_final,post_subprime_R_u['SPY']],axis = 1)
post_subprime_spy_performance = analytics(X = pd.DataFrame(post_subprime_R_u.loc[:,'SPY']),rf = 0.06, confidenceLevel = 0.95, position = 100)
post_subprime_exec = pd.concat([post_subprime_exec,post_subprime_spy_performance],axis = 1)

post_subprime_exec.columns = [['$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$',
                                  '$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$',
                                  '$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$',
                                  '$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$',
                                  '$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$',
                                  '$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','SPY'],
                                 ['\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0','']]

post_subprime_final.columns = post_subprime_exec.columns

# Save Data:
save_data(post_subprime_exec, 'post_subprime_exec')

"""#### Plots

##### Notes

An investment strategy is abbreviated as $S_{Cov}^{E[r]}(\beta)$. In this implementation, *post_subprime_lookbacks* contains pairs $[E[r], Cov]$. *post_subprime_betas* contains the various target $\beta$.

##### Total Value
"""

# After Crisis Plot:
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111)
for i in range(36):
  ax.plot(100*(np.cumprod(post_subprime_final.iloc[:,i]+1)),label = post_subprime_final.columns[i][0]+', '+post_subprime_final.columns[i][1])
  ax.legend(loc='best', ncol=4, fontsize=10)
plt.xlabel('t', fontweight='bold', fontsize=15)
plt.ylabel('Value', fontweight='bold', fontsize=15)
plt.title('Value of Investment Strategies After the Subprime Crisis', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '14_post_subprime_strategy_val.png', bbox_inches='tight')

"""##### $S_{60}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = post_subprime_final.iloc[:,i]
  col_name = post_subprime_final.columns[i]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/30,[z]*len(y),dens/10,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{60}$ Returns After the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '15_post_subprime_ret_distS6060.png', bbox_inches='tight')

"""##### $S_{120}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = post_subprime_final.iloc[:,i+6]
  col_name = post_subprime_final.columns[i+6]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/40,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('r', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{60}$ Returns After the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '16_post_subprime_ret_distS60120.png', bbox_inches='tight')

"""##### $S_{60}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = post_subprime_final.iloc[:,i+12]
  col_name = post_subprime_final.columns[i+12]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/30,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('r', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{90}$ Returns After the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '17_post_subprime_ret_distS9060.png', bbox_inches='tight')

"""##### $S_{120}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = post_subprime_final.iloc[:,i+18]
  col_name = post_subprime_final.columns[i+18]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/30,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{90}$ Returns After the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '18_post_subprime_ret_distS90120.png', bbox_inches='tight')

"""##### $S_{60}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = post_subprime_final.iloc[:,i+24]
  col_name = post_subprime_final.columns[i+24]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/30,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{120}$ Returns After the Subprime Crisis', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '19_post_subprime_ret_distS12060.png', bbox_inches='tight')

"""##### $S_{120}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = post_subprime_final.iloc[:,i+30]
  col_name = post_subprime_final.columns[i+30]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/50,[z]*len(y),dens/9,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{120}$ Returns After the Subprime Crisis', fontweight='bold', fontsize=18)
+
plt.savefig(graphs_dir + '20_post_subprime_ret_distS120120.png', bbox_inches='tight')

"""### Before COVID-19
Period: January 2, 2015 - March 9, 2020

#### Tearsheet
"""

# Pre-COVID:
pre_covid_R_u = R_etf.loc['1/2/2015':'3/9/2020',:'FEZ']
pre_covid_ff_factors = ff_3_daily.loc['1/2/2015':'3/9/2020','Mkt-RF':'RF']
pre_covid_lookbacks = [[60,60], [60,120], [90,60], [90,120], [120,60], [120,120]]
pre_covid_betas = [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0]
pre_covid_exec = pd.DataFrame([])
pre_covid_final = pd.DataFrame([])
omegas = []
for lb in pre_covid_lookbacks:
  for bt in pre_covid_betas:
    res = backtesting(pre_covid_R_u,
                    pre_covid_ff_factors,
                    return_period = lb[0],
                    variance_period = lb[1],
                    lamb = 10,
                    beta_tm = bt)
    omegas.append(res[1])
    res = pd.DataFrame(res[0],index = pd.to_datetime(pre_covid_R_u.index))
    res_perf = analytics(X = res,rf = 0.06, confidenceLevel = 0.95, position = 100)
    pre_covid_final = pd.concat([pre_covid_final,res],axis = 1)
    pre_covid_exec = pd.concat([pre_covid_exec,res_perf],axis = 1)

pre_covid_final = pd.concat([pre_covid_final,pre_covid_R_u['SPY']],axis = 1)
pre_covid_spy_performance = analytics(X = pd.DataFrame(pre_covid_R_u.loc[:,'SPY']),rf = 0.06, confidenceLevel = 0.95, position = 100)
pre_covid_exec = pd.concat([pre_covid_exec,pre_covid_spy_performance],axis = 1)

pre_covid_exec.columns = [['$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$',
                                  '$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$',
                                  '$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$',
                                  '$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$',
                                  '$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$',
                                  '$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','SPY'],
                                 ['\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0','']]

pre_covid_final.columns = pre_covid_exec.columns

# Save Data:
save_data(pre_covid_exec, 'pre_covid_exec')

"""#### Plots

##### Notes

An investment strategy is abbreviated as $S_{Cov}^{E[r]}(\beta)$. In this implementation, *post_subprime_lookbacks* contains pairs $[E[r], Cov]$. *post_subprime_betas* contains the various target $\beta$.

##### Total Value
"""

# Pre-COVID Plot:
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111)
for i in range(36):
  ax.plot(100*(np.cumprod(pre_covid_final.iloc[:,i]+1)),label = pre_covid_final.columns[i][0]+', '+pre_covid_final.columns[i][1])
  ax.legend(loc='best', ncol=4, fontsize=10)
plt.xlabel('t', fontweight='bold', fontsize=15)
plt.ylabel('Value', fontweight='bold', fontsize=15)
plt.title('Value of Investment Strategies Before COVID-19', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '21_pre_covid_strategy_val.png', bbox_inches='tight')

"""##### $S_{60}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_covid_final.iloc[:,i]
  col_name = pre_covid_final.columns[i]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/20,[z]*len(y),dens/10,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{60}$ Returns Before COVID-19', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '22_pre_covid_ret_distS6060.png', bbox_inches='tight')

"""##### $S_{120}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_covid_final.iloc[:,i+6]
  col_name = pre_covid_final.columns[i+6]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/20,[z]*len(y),dens/10,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{60}$ Returns Before COVID-19', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '23_pre_covid_ret_distS60120.png', bbox_inches='tight')

"""##### $S_{60}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_covid_final.iloc[:,i+12]
  col_name = pre_covid_final.columns[i+12]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/20,[z]*len(y),dens/10,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{90}$ Returns Before COVID-19', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '24_pre_covid_ret_distS9060.png', bbox_inches='tight')

"""##### $S_{120}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_covid_final.iloc[:,i+18]
  col_name = pre_covid_final.columns[i+18]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/20,[z]*len(y),dens/10,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{90}$ Returns Before COVID-19', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '25_pre_covid_ret_distS90120.png', bbox_inches='tight')

"""##### $S_{60}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_covid_final.iloc[:,i+24]
  col_name = pre_covid_final.columns[i+24]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/20,[z]*len(y),dens/10,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{120}$ Returns Before COVID-19', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '26_pre_covid_ret_distS12060.png', bbox_inches='tight')

"""##### $S_{120}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = pre_covid_final.iloc[:,i+30]
  col_name = pre_covid_final.columns[i+30]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot/20,[z]*len(y),dens/10,color ='black',linewidth = 2.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{120}$ Returns Before COVID-19', fontweight='bold', fontsize=18)

plt.savefig(graphs_dir + '27_pre_covid_ret_distS120120.png', bbox_inches='tight')

"""### During COVID-19
Period: March 9, 2020 - Present

#### Tearsheet
"""

# During-COVID:
on_covid_R_u = R_etf.loc['3/9/2020':,:'FEZ']
on_covid_ff_factors = ff_3_daily.loc['3/9/2020':,'Mkt-RF':'RF']
on_covid_lookbacks = [[60,60], [60,120], [90,60], [90,120], [120,60], [120,120]]
on_covid_betas = [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0]
on_covid_exec = pd.DataFrame([])
on_covid_final = pd.DataFrame([])
omegas = []
for lb in on_covid_lookbacks:
  for bt in on_covid_betas:
    res = backtesting(on_covid_R_u,
                    on_covid_ff_factors,
                    return_period = lb[0],
                    variance_period = lb[1],
                    lamb = 10,
                    beta_tm = bt)
    omegas.append(res[1])
    res = pd.DataFrame(res[0],index = pd.to_datetime(on_covid_R_u.index))
    res_perf = analytics(X = res,rf = 0.06, confidenceLevel = 0.95, position = 100)
    on_covid_final = pd.concat([on_covid_final,res],axis = 1)
    on_covid_exec = pd.concat([on_covid_exec,res_perf],axis = 1)

on_covid_final = pd.concat([on_covid_final,on_covid_R_u['SPY']],axis = 1)
on_covid_spy_performance = analytics(X = pd.DataFrame(on_covid_R_u.loc[:,'SPY']),rf = 0.06, confidenceLevel = 0.95, position = 100)
on_covid_exec = pd.concat([on_covid_exec,on_covid_spy_performance],axis = 1)

on_covid_exec.columns = [['$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$',
                                  '$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$',
                                  '$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$',
                                  '$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$',
                                  '$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$',
                                  '$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','SPY'],
                                 ['\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0','']]

on_covid_final.columns = on_covid_exec.columns

# Save Data:
save_data(on_covid_exec, 'on_covid_exec')

"""#### Plots

##### Notes

An investment strategy is abbreviated as $S_{Cov}^{E[r]}(\beta)$. In this implementation, *post_subprime_lookbacks* contains pairs $[E[r], Cov]$. *post_subprime_betas* contains the various target $\beta$.

##### Total Value
"""

# During COVID Plot:
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111)
for i in range(36):
  ax.plot(100*(np.cumprod(on_covid_final.iloc[:,i]+1)),label = on_covid_final.columns[i][0]+', '+on_covid_final.columns[i][1])
  ax.legend(loc='best', ncol=4, fontsize=10)
plt.xlabel('t', fontweight='bold', fontsize=15)
plt.ylabel('Value', fontweight='bold', fontsize=15)
plt.title('Value of Investment Strategies During COVID-19', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '28_on_covid_strategy_val.png', bbox_inches='tight')

"""##### $S_{60}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = on_covid_final.iloc[:,i]
  col_name = on_covid_final.columns[i]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{60}$ Returns During COVID-19', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '29_on_covid_ret_distS6060.png', bbox_inches='tight')

"""##### $S_{120}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = on_covid_final.iloc[:,i+6]
  col_name = on_covid_final.columns[i+6]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{60}$ Returns During COVID-19', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '30_on_covid_ret_distS60120.png', bbox_inches='tight')

"""##### $S_{60}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = on_covid_final.iloc[:,i+12]
  col_name = on_covid_final.columns[i+12]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{90}$ Returns During COVID-19', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '31_on_covid_ret_distS9060.png', bbox_inches='tight')

"""##### $S_{120}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = on_covid_final.iloc[:,i+18]
  col_name = on_covid_final.columns[i+18]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{90}$ Returns During COVID-19', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '32_on_covid_ret_distS90120.png', bbox_inches='tight')

"""##### $S_{60}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = on_covid_final.iloc[:,i+24]
  col_name = on_covid_final.columns[i+24]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{120}$ Returns During COVID-19', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '33_on_covid_ret_distS12060.png', bbox_inches='tight')

"""##### $S_{120}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')
for i in range(6):
  dt = on_covid_final.iloc[:,i+30]
  col_name = on_covid_final.columns[i+30]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)
  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{120}$ Returns During COVID-19', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '34_on_covid_ret_distS120120.png', bbox_inches='tight')

"""### Full Investment Horizon
Period: March 22, 2007 - October 30, 2020

#### Tearsheet
"""

# Full Period:
full_horizon_R_u = R_etf.loc[:,:]
full_horizon_ff_factors = ff_3_daily.loc[:,'Mkt-RF':'RF']
full_horizon_lookbacks = [[60,60], [60,120], [90,60], [90,120], [120,60], [120,120]]
full_horizon_betas = [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0]
full_horizon_exec = pd.DataFrame([])
full_horizon_final = pd.DataFrame([])
omegas = []

for lb in full_horizon_lookbacks:
  for bt in full_horizon_betas:
    res = backtesting(full_horizon_R_u,
                    full_horizon_ff_factors,
                    return_period = lb[0],
                    variance_period = lb[1],
                    lamb = 10,
                    beta_tm = bt)
    omegas.append(res[1])
    res = pd.DataFrame(res[0],index = pd.to_datetime(full_horizon_R_u.index))
    res_perf = analytics(X = res,rf = 0.06, confidenceLevel = 0.95, position = 100)
    full_horizon_final = pd.concat([full_horizon_final,res],axis = 1)
    full_horizon_exec = pd.concat([full_horizon_exec,res_perf],axis = 1)

full_horizon_final = pd.concat([full_horizon_final,full_horizon_R_u['SPY']],axis = 1)
full_horizon_spy_performance = analytics(X = pd.DataFrame(full_horizon_R_u.loc[:,'SPY']),rf = 0.06, confidenceLevel = 0.95, position = 100)
full_horizon_exec = pd.concat([full_horizon_exec,full_horizon_spy_performance],axis = 1)
full_horizon_exec.columns = [['$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$','$S^{60}_{60}$',
                                  '$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$','$S^{60}_{120}$',
                                  '$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$','$S^{90}_{60}$',
                                  '$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$','$S^{90}_{120}$',
                                  '$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$','$S^{120}_{60}$',
                                  '$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','$S^{120}_{120}$','SPY'],
                                 ['\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0',
                                  '\(\beta\)=-1.0','\(\beta\)=-0.5','\(\beta\)=0.5','\(\beta\)=1.0','\(\beta\)=1.5','\(\beta\)=2.0','']]

full_horizon_final.columns = full_horizon_exec.columns
full_horizon_exec.to_csv('full_horizon_exec')

# Save Data:
save_data(full_horizon_exec, 'full_horizon_exec')

"""#### Plots

##### Notes

An investment strategy is abbreviated as $S_{Cov}^{E[r]}(\beta)$. In this implementation, *post_subprime_lookbacks* contains pairs $[E[r], Cov]$. *post_subprime_betas* contains the various target $\beta$.

##### Total Value
"""

fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111)
for i in range(36):
  ax.plot(100*(np.cumprod(full_horizon_final.iloc[:,i]+1)),label = full_horizon_final.columns[i][0]+', '+full_horizon_final.columns[i][1])
  ax.legend(loc='best', ncol=4, fontsize=10)
plt.xlabel('t', fontweight='bold', fontsize=15)
plt.ylabel('Value', fontweight='bold', fontsize=15)
plt.title('Value of Investment Strategies Across the Investment Horizon', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '35_full_horizon_strategy_val.png', bbox_inches='tight')

"""##### $S_{60}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = full_horizon_final.iloc[:,i]
  col_name = full_horizon_final.columns[i]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{60}$ Returns Across the Investment Horizon', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '36_full_ret_distS6060.png', bbox_inches='tight')

"""##### $S_{120}^{60}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = full_horizon_final.iloc[:,i+6]
  col_name = full_horizon_final.columns[i+6]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{60}$ Returns Across the Investment Horizon', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '37_full_ret_distS60120.png', bbox_inches='tight')

"""##### $S_{60}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = full_horizon_final.iloc[:,i+12]
  col_name = full_horizon_final.columns[i+12]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{90}$ Returns Across the Investment Horizon', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '38_full_ret_distS9060.png', bbox_inches='tight')

"""##### $S_{120}^{90}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = full_horizon_final.iloc[:,i+18]
  col_name = full_horizon_final.columns[i+18]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{90}$ Returns Across the Investment Horizon', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '39_full_ret_distS90120.png', bbox_inches='tight')

"""##### $S_{60}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = full_horizon_final.iloc[:,i+24]
  col_name = full_horizon_final.columns[i+24]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{60}^{120}$ Returns Across the Investment Horizon', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '40_full_ret_distS12060.png', bbox_inches='tight')

"""##### $S_{120}^{120}(\beta^{m}_{T})$"""

fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(111, projection='3d')

for i in range(6):
  dt = full_horizon_final.iloc[:,i+30]
  col_name = full_horizon_final.columns[i+30]
  c =  ['r', 'g', 'b', 'y', 'm','orange'][i]
  z =  [-1.0, -0.5, 0.5, 1.0, 1.5, 2.0][i]
  x,y = np.histogram(dt,bins = 100)
  x = x/len(dt)
  y = (y[:-1]+y[1:])/2
  cs = [c] * len(x)
  ax.bar(y, x, zs=z, zdir='y', color=cs, alpha=0.7,width = 0.003,label = col_name[0]+', '+col_name[1])
  ax.legend(fontsize=13)

  samples = np.asarray(dt).reshape(-1,1)
  x_plot = np.linspace(-10,10,100).reshape(-1,1)
  kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(samples)
  log_dens = kde.score_samples(x_plot)
  dens = np.exp(log_dens)
  ax.view_init(20, 50)
  ax.plot(x_plot / 20, [z] * len(y), dens / 8, color = 'black', linewidth = 3.0)
ax.set_xlabel('$ρ$', fontweight='bold', fontsize=15)
ax.set_ylabel('$\(\beta\)$', fontweight='bold', fontsize=15)
ax.set_zlabel('$f$', fontweight='bold', fontsize=15)
ax.set_title('$S_{120}^{120}$ Returns Across the Investment Horizon', fontweight='bold', fontsize=18)
plt.savefig(graphs_dir + '41_full_ret_distS120120.png', bbox_inches='tight')

"""## **Appendix**

Included is the progress on a Maximum Return solver and a Minimum Variance solver with $\rho_{p}=15\%$

### Data
"""

# ETF Data:
ticker = ['FXE','EWJ','GLD','QQQ','SPY','SHV','DBA','USO',
          'XBI','ILF','GAF','EPP','FEZ']
start =  datetime(2007, 3, 26)
end = datetime(2020, 10, 30)
data = pd.DataFrame()
for i in ticker:
    data[i] = web.DataReader(i, 'yahoo', start, end)['Close']
data.to_csv('ETFs.csv')

# Load data
ETF = pd.read_csv(data_dir + 'ETFs.csv', index_col=0)[55:]
F = pd.read_csv(data_dir + 'Factors.csv', index_col=0)[56:]
F.index = ETF.index[1:]
# Calculte the simple anualized returns for the ETFs
R = (ETF.pct_change(1)[1:])*250
# Calculate the excess annualized return for the ETFs
ER = pd.DataFrame(R.values-F['RF'].values.reshape(-1,1),
                  index=F.index, columns=ticker)
print(str(ER))
F = F.iloc[:,0:3]
F

"""### Maximum Return"""

# Before the subprime crisis(2007/03/26 - 2008/03/23)
R_bc = R['2007-03-26':'2008-03-23'].values
ER_bc = ER['2007-03-26':'2008-03-23'].values
F_bc = F['2007-03-26':'2008-03-23'].values

# During the subprime crisis(2008/03/24 - 2009/06/30)
R_bc = R['2008-03-24':'2009-06-30'].values
ER_bc = ER['2008-03-24':'2009-06-30'].values
F_bc = F['2008-03-24':'2009-06-30'].values

# After the subprime crisis(2007/03/26 - 2008/03/23)
#R_bc = R['2009-06-30':'2016-10-20'].values
#ER_bc = ER['2009-06-30':'2016-10-20'].values
#F_bc = F['2009-06-30':'2016-10-20'].values

Num_days = len(F_bc)
FR_bc = F_bc[1:].copy()

# Short term model(60 days)
Lambda = 0.001
beta_T = [0.5, 1, 1.5]
R_opt = []

# Conduct the max return strategy:
for j in beta_T:
  Rp = []
  wp = np.ones((13,1))*1/13
  for i in range(len(R_bc)-59):
    r = R_bc[i:(i+60),:]
    er = ER_bc[i:(i+60),:]
    f1 = F_bc[i:(i+60),:]
    rho = r.mean(axis=0).reshape(-1,1)
    cov_f = np.cov(f1, rowvar=False)

    # Run regression to get the beta:
    lm = LinearRegression()
    lm.fit(f1, er)
    coeff3 = lm.coef_
    beta = coeff3[:,0]
    error = er - lm.predict(f1)

    # Calculate the covariance matrix:
    Q = coeff3.dot(cov_f).dot(coeff3.T)+np.diag(error.var(axis=0))

    # Preparation for the optimization:
    P = matrix(2*Lambda*Q, tc='d')
    q = matrix(-2*Lambda*(Q.T).dot(wp)-rho, tc='d')
    A = matrix(np.vstack((beta, [1]*13)), tc='d')
    G = matrix(np.vstack((np.diag([1]*13),np.diag([-1]*13))), tc='d')
    h = matrix([2]*26, tc='d')
    b = matrix([j,1], tc='d')

    # Do the optimization using QP solver:
    opt = solvers.qp(P, q, G, h, A, b, options={'show_progress':False})
    w = opt['x']
    wp = np.array(w).reshape(-1,1)
    Rp = Rp + [wp.T.dot(rho)[0,0]]
  R_opt.append(Rp)

R_opt = pd.DataFrame(np.array(R_opt))
R_opt

# Short term model(60 days)
Lambda = 0.001
beta_T = [0.5, 1, 1.5]
R_opt = []

# Conduct the max return strategy
window = 63
alocate = 5
R_opt = R_bc[window:,4]/250
for j in beta_T:
  Rp = []
  wp = np.ones((13,1))*1/13
  for i in range(window,Num_days):
    future_return = R_bc[i, :].reshape(-1, 1)
    if i%alocate==0:
      r = R_bc[(i-window):i,:]
      er = ER_bc[(i-window):i,:]
      f1 = F_bc[(i-window):i,:]
      rho = r.mean(axis=0).reshape(-1,1)
      cov_f = np.cov(f1, rowvar=False)

      # Run regression to get the beta
      lm = LinearRegression()
      lm.fit(f1, er)
      coeff3 = lm.coef_
      beta = coeff3[:,0]
      error = er - lm.predict(f1)

      # Calculate the covariance matrix
      #Q = coeff3.dot(cov_f).dot(coeff3.T)+np.diag(error.var(axis=0))
      Q = np.diag([1]*13)

      # Preparation for the optimization
      P = matrix(2*Lambda*Q, tc='d')
      q = matrix(-2*Lambda*(Q.T).dot(wp)-rho, tc='d')
      A = matrix(np.vstack((beta, [1]*13)), tc='d')
      G = matrix(np.vstack((np.diag([1]*13),np.diag([-1]*13))), tc='d')
      h = matrix([2]*26, tc='d')
      b = matrix([j,1], tc='d')

      # Do the optimization using QP solver
      opt = solvers.qp(P, q, G, h, A, b, options={'show_progress':False})
      w = opt['x']
      wp = np.array(w).reshape(-1,1)
    Rp = Rp + [wp.T.dot(future_return/250)[0,0]]

R_opt = pd.DataFrame(np.array(R_opt).transpose())
R_opt

"""### Minimum variance with $\rho_{p}=15\%$"""

# Before the subprime crisis(2007/03/26 - 2008/03/23)
R_bc = R['2007-03-26':'2008-03-23'].values
ER_bc = ER['2007-03-26':'2008-03-23'].values
F_bc = F['2007-03-26':'2008-03-23'].values

# During the subprime crisis(2008/03/24 - 2009/06/30)
R_bc = R['2008-03-24':'2009-06-30'].values
ER_bc = ER['2008-03-24':'2009-06-30'].values
F_bc = F['2008-03-24':'2009-06-30'].values

# After the subprime crisis(2007/03/26 - 2008/03/23)
#R_bc = R['2009-06-30':'2016-10-20'].values
#ER_bc = ER['2009-06-30':'2016-10-20'].values
#F_bc = F['2009-06-30':'2016-10-20'].values

Num_days = len(F_bc)
FR_bc = F_bc[1:].copy()

Rp = []
wp = np.ones((13,1))*1/13

for i in range(len(R_bc)-59):
  r = R_bc[i:(i+60),:]
  er = ER_bc[i:(i+60),:]
  f2 = F_bc[i:(i+60),:]
  rho = r.mean(axis=0)
  cov_f = np.cov(f2, rowvar=False)

  # Run regression to get the beta
  lm = LinearRegression()
  lm.fit(f2, er)
  coeff3 = lm.coef_
  beta = coeff3[:,0]
  error = er - lm.predict(f2)

  # Calculate the covariance matrix
  Q = coeff3.dot(cov_f).dot(coeff3.T)+np.diag(error.var(axis=0))

  # Preparation for the optimization
  P = matrix(2*(1+Lambda)*Q, tc='d')
  q = matrix(-2*Lambda*(Q.T).dot(wp), tc='d')
  G = matrix(np.vstack((np.diag([1]*13),np.diag([-1]*13))), tc='d')
  h = matrix([2]*26, tc='d')
  A = matrix(np.vstack((beta, [1]*13)), tc='d')
  b = matrix([0.15,1], tc='d')

  # Do the optimization using QP solver
  opt = solvers.qp(P, q, G, h, A, b, options={'show_progress':False})
  w = opt['x']
  wp = np.array(w).reshape(-1,1)
  Rp = Rp + [wp.T.dot(rho.reshape(-1,1))[0,0]]

R_opt.append(Rp)
R_opt = pd.DataFrame(np.array(R_opt))
R_opt

# Conduct the min variance with 15% target return strategy.
Rp = []
wp = np.ones((13,1))*1/13

# Short term model(60 days)
Lambda = 0.001
beta_T = [0.5, 1, 1.5]
R_opt = []

for i in range(window, Num_days):
  future_return = R_bc[i, :].reshape(-1, 1)
  if i % alocate == 0:
    r = R_bc[(i - window):i, :]
    er = ER_bc[(i - window):i, :]
    f1 = F_bc[(i - window):i, :]
    rho = r.mean(axis=0)
    cov_f = np.cov(f1, rowvar=False)

    # Run regression to get the beta
    lm = LinearRegression()
    lm.fit(f1, er)
    coeff3 = lm.coef_
    beta = coeff3[:,0]
    error = er - lm.predict(f1)

    # Calculate the covariance matrix
    Q = coeff3.dot(cov_f).dot(coeff3.T)+np.diag(error.var(axis=0))
    Q_ = np.diag([1]*13)

    # Preparation for the optimization
    P = matrix((Q+Lambda*Q_), tc='d')
    q = matrix(-2*Lambda*(Q_.T).dot(wp), tc='d')
    G = matrix(np.vstack((np.diag([1]*13),np.diag([-1]*13))), tc='d')
    h = matrix([2]*26, tc='d')
    A = matrix(np.vstack((rho, [1]*13)), tc='d')
    b = matrix([0.15,1], tc='d')

    # Do the optimization using QP solver
    opt = solvers.qp(P, q, G, h, A, b, options={'show_progress':False})
    w = opt['x']
    wp = np.array(w).reshape(-1,1)
  Rp = Rp + [wp.T.dot(future_return/250)[0,0]]
  R_opt.append(Rp)

#plt.plot(range(result.shape[0]),result['beta=0.5'])

result = pd.DataFrame(R_opt)
print(result)

"""### PnL"""

# Compute PnL
pnl = PnL(result)
print(pnl)
for i in range(6):
  plt.plot(pnl[:,i],label=i)
plt.legend(loc='best')
plt.show()

# result = R_bc.copy()/250
days = result.shape[0]

print('Last PnL after %s: ' % days, PnL(result,100)[-1, :])

# Geometric Mean
print('Geometric mean',geom_mean(result)*250)

# Min
print('Daily min',np.min(result,axis=0)*250)

# Max Drawdown
print('max drawdown: ', MaxDrawdown(result))

# Vol:
print('Volatility', Volatility(result))

# Sharpe Ratio:
RF = np.array(R_bc-ER_bc)[:,0].reshape(-1,1)/250
#print('Sharp ratio: ', Sharpe(result,RF))
# print('Mean sharp: ', np.mean(Sharpe(result,RF),axis=0))

# Kurt:
print('Kurt: ', Kurt(result))
print('Skewness: ', Skew(result))
print('%s VaR %s days: ' % (0.99, days), VaR(result, 0.99))
print('%s CVaR %s days: ' % (0.99, days), CVaR(result, 0.99))

#for i in range(result.shape[1]):
    # print(i)
    #plt.plot((1+result[:,i]).cumprod(),label=ticker[i])

#plt.legend(loc='best')
RF = np.array(R_bc - ER_bc)[window:, 0].reshape(-1, 1) / 250
#Summary(R_bc,RF,0.99)

"""### Other"""

pd.DataFrame(R_opt,index=['\(\beta\)=0.5','\(\beta\)=1','\(\beta\)=1.5','minvar']).T

r = R_bc[0:60,:]
er = ER_bc[0:60,:]
f1 = F_bc[0:60,:]
rho = r.mean(axis=0).reshape(-1,1)
cov_f = np.cov(f1, rowvar=False)

# Run regression to get the beta
lm = LinearRegression()
lm.fit(f1, er)
coeff3 = lm.coef_
beta = coeff3[:,0]
error = er - lm.predict(f1)
# Calculate the covariance matrix
Q = coeff3.dot(cov_f).dot(coeff3.T)+np.diag(error.var(axis=0))
# Preparation for the optimization
P = matrix(2*Lambda*Q, tc='d')
q = matrix(-2*Lambda*(Q.T).dot(wp)-rho, tc='d')
A = matrix(np.vstack((beta, [1]*13)), tc='d')
b = matrix([1.5,1], tc='d')
# Do the optimization using QP solver
opt = solvers.qp(P=P, q=q, A=A, b=b, options={'show_progress':False})
w = opt['x']
wp = np.array(w).reshape(-1,1)

wp.T.dot(rho)
wp

pd.DataFrame(R_opt,index=['\(\beta\)=0.5','\(\beta\)=1','\(\beta\)=1.5','minvar']).T
	

DOCUMENTS

The Github repository can be found here.
All requests for copies of the research, please forward to my university email address listed on my homepage.