<jupyter_start><jupyter_text>Portfolio optimisation with Python> "In this article, we provide a bird's eye perspective on how to use Python and a dedicated library called PyPortfolioOpt to put Modern Portfolio Optimisation techniques into practice. In particular, we look at the approach developed in 1952 by <NAME>."- toc: false- badges: true- comments: false- author: <NAME>.- categories: [portfolio-optimisation]- tags: [tutorial, markowitz, optimisation, portfolio, construction]- image: images/portfolio-optimisation.png- series: "Finance and Markets with Python and Jupyter: a step-by-step tutorial" ![](photos/portfolio-optimisation-banner.png) > Tip: This article is written as a Jupyter Notebook. It has been published using [Fastpages](https://portfolio-geek.com/blogging/publishing/workflow/2022/03/06/write-articles-with-jupyter-fastpages.html). The Jupyter notebook is available on GitHub and if you want to, you can run it directly using the provided Binder link displayed at the top of the article.{% series_list%} Introduction**<NAME>** is one of the (if not "the") fathers of modern portfolio construction and his seminal paper Portfolio Selection{% fn 1 %} has driven an entire research stream. The **intuition** behind this paper is that one can combine the information gathered on expected returns, risks and diversification of various assets or asset classes with a view to optimise the risk-return profile of a given portfolio.In practice, we typically find a limited stability of the portfolios generated using this approach, especially due the high sensitivity of the portfolio to the expected returns; nonetheless, the insights it provides are very useful and it's a perfect start to see the impact of risk and diversification.For this article, we will mostly rely on a fantastic Python library, PyPortfolioOpt {% fn 2%} which will do the optimisation heavy lifting for us.  Setup Loading LibrariesThe main library to load here is indeed **PyPortfolioOpt**{% fn 2%}, and we will rely extensively on it in this article. We need the following tools as well:- **Pandas**{% fn 3%}: load, filter, sort and pretty much all data wrangling operations;- **Numpy**{% fn 4%}: provides most matrix and advanced numerical operations; this library is the calculation backbone for pandas;- **Matplotlib**{% fn 5%}: the *de facto* reference library to draw scientific charts;- **Pandas_datareader**{% fn 6%}: a very handy library to access many different online databases, including Yahoo Finance.<jupyter_code>#collapse-hide
import numpy as np 
import pandas_datareader.data as web_reader
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import math
from pypfopt.efficient_frontier import EfficientFrontier
from pypfopt import risk_models
from pypfopt import expected_returns<jupyter_output><empty_output><jupyter_text>Misc ParametersIn the code below, we set some variables to adjust the jupyter loo<jupyter_code>#collapse-hide
# Note that this change the decimals places inside Jupyter, but not on the website
pd.options.display.float_format = '{:,.1f}'.format<jupyter_output><empty_output><jupyter_text>Defining our Investable UniverseFirst and foremost, we need to define our **Investable Universe** i.e. the set of asset classes that we will allow in our portfolio. Many investors would typically think about single stocks, but on my side, I am typically looking I am looking here a Long Term Investing, and I would be keen to Speaking about asset classes, we need some ETFs to analyse!> "Important: Please bear in mind this article if purely for pedagogical purpose and should by no mean be understood as a recommendation or advice. Investing brings risk and in particular risk of loss of capital. I have no intention to recommend anything! I will follow among others this recent article in US News{% fn 7%} and ETF.com{% fn 8%} as sources of inspiration to identify relevant ETFs." **Equity**Let's start with various Equity sub-asset classes:- **IVV**: iShares Core **S&P 500** ETF representing US Equity Large Cap- **SCHA**: Schwab **U.S. Small-Cap** ETF representing US Equity Small Cap- **IJH**: iShares Core **S&P Mid-Cap** ETF representing US Equity Mid Cap- **SCHD**: Schwab **U.S. Dividend Equity** ETF representing US Equity Dividend (ie US stocks which are deemed to pay higher dividends)- **VTI**: Vanguard **Total Stock Market** ETF representing US Equity, with all market cap included- **VXUS**: Vanguard **Total International Stock** ETF representing World ex-US equities- **EEM**: iShares **MSCI Emerging Markets** ETF representing the Emerging Market Equities **Fixed Income**Let's add Fixed Income, which would typically aim to reduce the overall portfolio's volatility:- **AGG**: iShares Core **U.S. Aggregate Bond** ETFrepresenting the entire US Bond market- **GOVT**: iShares **U.S. Treasury Bond** ETF representing the performance of US Government Bonds- **VCLT**: Vanguard **Long-Term Corporate Bond** ETF representing the Investment Grade USD denominated bonds. **Alternative Assets**Let's add 2 additional asset classes, Gold and Commodities:- **GLD**: SPDR **Gold** Trust representing the price of **Gold- **PDBC**: Invesco Optimum Yield Diversified **Commodity** Strategy No K-1 ETF** which will represent the performance of the Broad Commodities asset classEverything in the above is **denominated in US Dollars**, this will make our life easier in what follows, ie we will not need any currency conversion, which is always a bit painfull in the process.  Parameters and data gatheringI recommend you to read this **[article about data gathering](https://portfolio-geek.com/data/yahoo%20finance/pandas/2022/02/17/gathering_etf_and_stock_historical_data.html)**.Like we did in this article, we will utilise **pandas_datareader**{% fn 6%} to get historical time series. As mentioned above, we will look at the performance of ETFs, which we will consider as as relevant proxies for their respective asset classes.You can of course utilise whatever asset class you want, and take single stocks, funds as historical data points. On my side, I am quite familiar with indices ETFs, and I will go with the selection above. Last but not least, we need to define the **start_date** and **end_date** for gathering the historical time series. For this study, we will gather 7 years of data.<jupyter_code>#collapse-hide

start_date = '2015-03-01'
end_date = '2022-03-01'

# Define Investable Universe
investable_universe_tickers = ['IVV', 'SCHA', 'IJH', 'SCHD', 'VTI', 'VXUS', 'EEM', 'VCLT', 'AGG', 'GOVT','PDBC','GLD']

# Get Historical Data
df = web_reader.DataReader(investable_universe_tickers, 'yahoo', start=start_date, end=end_date)
df = df['Adj Close']

# This is required to round the blog's table into 2 decimals, Jupyter's formatting does not apply on the published website
df = df.round(decimals=1)<jupyter_output><empty_output><jupyter_text>Data gathering: **results**The request above delivered a **pandas data_frame**, and here is a snapshot of the last 5 rows:<jupyter_code>#collapse-hide
df.tail(3)<jupyter_output><empty_output><jupyter_text>NormalisationA table like the above is not very useful: in the absence a of particle knowledge of the ETFs' values, we have no way of knowing if a value is "high" or "low", hence apart from telling us that the value is a number we have now way to let's try to make this table a bit more insightful.The request above delivered a data_frame, and here is a snapshot of the last 5 rows:<jupyter_code>#collapse-hide
# Same table, but this time, normalised
(df/df.iloc[0, ]*100).round(decimals=1).tail(3)<jupyter_output><empty_output><jupyter_text>This is **much** better.Now we can at least see that:- **US Equities** (e.g. IVV) had a **fantastic ride since 2015**;- Government bonds (e.g. GOVT) under-performed;- Commodities (PDBC) and Gold (GLD) had several rough years as well, with recent massive pick-up in the current geopolitical context. Visualise asset classes' returnsBeing able to sanity check the data is very important, and it's often more efficient with a quick chart.<jupyter_code>#collapse-hide

# Normalise to 100
normalised_data = (df/df.iloc[0, ]*100)

# A bit of data wrangling
transposed_data = normalised_data.tail(1).reset_index().transpose()
transposed_data = transposed_data.iloc[1:len(transposed_data)]
transposed_data = transposed_data.rename(columns={0: "Last Value"})

# We want a bar chart sorted by decreaseing values
transposed_data = transposed_data.sort_values(by="Last Value", ascending=False)

# Theming Seaborn results
sns.set_theme()

# Draw
transposed_data.plot.bar(figsize=(10, 6))
pass<jupyter_output><empty_output><jupyter_text>Draw the **wealth curve**The table above is useful, but when it comes to grasping and long term risks and returns, a chart is worth a thousand words.A very common issue when charting multiple time series is the very different stock levels, and this can make the chart hard to read.This is why we will once again normalise the data.<jupyter_code>#collapse-hide

# Theming Seaborn results
sns.set_theme()

# Plot the time series
plt.figure(figsize=(12,6))

# Legends and Axis titles
p = sns.lineplot(data=normalised_data)
p.set_ylabel(f"Close Price, Basis 100 in {start_date}")
pass<jupyter_output><empty_output><jupyter_text>Estimate risk and returnIn Markowitz 1952{% fn 1 %}, the optimal portfolio is obtained as a function of *expected* returns and *expected* risks of the portfolio. This takes a strong assumption that we have a "crystal ball", or at least access to a predictive model, which at this stage is well beyond this article.For this first example, we will rely on historical parameters estimation, which precisely means that we expect the past to be a good prediction of what might happen in the future.   Volatility There are many ways to estimate the volatility, and I will only cover the simplest approach. We have access to daily close prices of the ETFs, hence we could look at the standard deviation of **daily** returns. In theory this would utilise most of the data we have access to, which is a good thing. In practice, and especially when looking at asset classes which can be observed with an **8 to 12 hours** time difference which might lead the "sample volatility" to become inconsistent between 2 asset classes. This is even more important for estimating the **correlation**.Practitioners often use **weekly returns** to alleviate this issue, and we will do the same here, and as such calculate the returns using a *resampled* time series. To calculate the **annualized volatility** requires an additional factor in this case the square root of 52. Why this? Volatility is essential to option traders, and when pricing options, practitioners often model asset prices as **Wiener processes**(number of weeks in a year) {% fn 9 %}. The variance of a Wiener process is proportional to the time, and the volatility is the square root of the the variance, which gives us that to convert standard deviation of weekly returns into an *annualized* figure, we need to multiply our results by the square root of the number of weeks in a year (more on Wiener processes here {% fn 9 %}).<jupyter_code>weekly_returns = df.resample("W").last().pct_change()
((weekly_returns.std()*math.sqrt(52)*100).sort_values()).plot.bar();<jupyter_output><empty_output><jupyter_text>As one could have guessed, on the left of the chart above, we can find the low volatility asset classes (e.g. Government Bonds, Investment Grade Bonds), whereas on the right we have high volatility asset classes (Mid Cap and Small Cap Equities).  Correlation<jupyter_code>sns.heatmap(weekly_returns.corr());
#collapse-hide
from pypfopt import risk_models
from pypfopt import plotting

# Calculate expected returns and sample covariance
mu = expected_returns.mean_historical_return(df)
sample_cov = risk_models.sample_cov(df, frequency=252)

S = risk_models.CovarianceShrinkage(df).ledoit_wolf()<jupyter_output><empty_output><jupyter_text>The main results: Efficient Frontier & Maximum Sharpe Portfolio Calculate and draw the efficient frontierWith expected returns and risk estimated, we are ready to utilise PyPortfolioOpt's optimiser to draw the **efficient frontier**.The idea behind the efficient frontier{% fn 1 %} is relatively simple: - For **each level of investor's risk**, there is an optimal portfolio which is expected to deliver the **highest level of return**; or conversely- For **each level of investor's return**, there is an optimal portfolio which is expected to deliver the **lowest level of return**.<jupyter_code>#collapse-hide
ef = EfficientFrontier(mu, sample_cov)

# We create 2 efficient frontiers
fig, ax = plt.subplots()
plotting.plot_efficient_frontier(ef, ax=ax, show_assets=True)
ax.set_title("Asset Classes & Efficient Frontier")

plt.show()<jupyter_output><empty_output><jupyter_text>Add the Maximum Sharpe portfolio<jupyter_code>#collapse-hide

# Find and plot the tangency portfolio
fig2, ax = plt.subplots()
ef2 = EfficientFrontier(mu, S) 
plotting.plot_efficient_frontier(ef2, ax=ax, show_assets=True)

ef3 = EfficientFrontier(mu, S) 

ef3.max_sharpe()
ret_tangent, std_tangent, _ = ef3.portfolio_performance()
ax.scatter(std_tangent, ret_tangent, marker="*", s=100, c="r", label="Max Sharpe")
ax.set_title("Asset Classes, Efficient Frontier & Max Sharpe portfolio")
plt.show()<jupyter_output><empty_output><jupyter_text>Calculate random portfolios<jupyter_code>#collapse-hide

# Plot random portfolios
fig2, ax = plt.subplots()
n_samples = 10000
w = np.random.dirichlet(np.ones(len(mu)), n_samples)
rets = w.dot(mu)
stds = np.sqrt((w.T * (S @ w.T)).sum(axis=0))
sharpes = rets / stds
ax.scatter(stds, rets, marker=".", c=sharpes, cmap="viridis_r")

# Format
ax.set_title("Random portfolios, based on the same asset classes")
#ax.legend()
plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08461649158152253, 'fraction_numerical': 0.0092099582673766, 'mean_word_length': 5.640707118967988, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 55, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '188843', 'license': 'Apache-2.0', 'n_tokens_mistral': 4103, 'n_tokens_neox': 3839, 'n_words': 1775, 'path': '_notebooks/2022-03-07-optimise-portfolios-with-python.ipynb', 'repo_name': 'vdenoise/portfolio-quant-tech', 'size': '13919'}