Workshop 3, Financial Modeling and Programming

Portfolio optimization and backtesting

Abstract

This is an INDIVIDUAL workshop. In the Algorithms course you simulated thousands of random portfolios and looked for the best one in the cloud. Here we find the optimal portfolios exactly: the Global Minimum Variance portfolio and the tangency (maximum Sharpe ratio) portfolio with matrix algebra, the efficient frontier and the Capital Market Line, and portfolios without short sales with a numerical optimizer. Then we learn why optimized portfolios are fragile (estimation risk) and we test our portfolio with data the model never saw: an out-of-sample backtest against the market and against an equally weighted portfolio.

1 General directions

Same rules as before: notebook W3-FinProg-YourFirstName-YourLastName, Predict → Run → Investigate → Modify → Make, pseudocode before the LLM, verification, and submission of the .ipynb file. Grade: notebook 60% + Canvas Quiz W3 40%.

import os, io, requests
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import minimize

def import_csv(file, base_url="https://www.apradie.com/datos/", **kwargs):
    """Imports a csv file from base_url; downloads it only once and keeps a local copy."""
    if not os.path.exists(file):
        headers = {"User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 "
                                 "(KHTML, like Gecko) Chrome/91.0.4472.124 Safari/537.36"}
        response = requests.get(base_url + file, headers=headers)
        response.raise_for_status()
        with open(file, "wb") as f:
            f.write(response.content)
    return pd.read_csv(file, **kwargs)

prices = import_csv("fz2024_prices_monthly.csv", index_col="date", parse_dates=True)
rf = import_csv("fz2024_rf_monthly.csv", index_col="date", parse_dates=True)
R = prices.pct_change(fill_method=None).iloc[1:]          # simple returns
r = np.log(prices).diff().iloc[1:]                        # cc returns

2 Inputs of portfolio optimization

We keep the conventions of the Algorithms course:

• Expected returns: the annual expected simple return of each asset, from the mean of its monthly cc returns: E[R_i]=e^{12\bar r_i}-1.
• Risk: the annual variance-covariance matrix of cc returns, \Sigma = 12\times Cov(r).
• Portfolio returns over time: weights times simple returns, R[ASSETS] @ w.

We use 5 assets as a running example (the same ones of the Algorithms course) and the 36-month window before our decision date (March 31, 2023):

ASSETS = ["AAPL", "JPM", "WMT", "XOM", "GLD"]
START, END = "2020-04-30", "2023-03-31"

def opt_inputs(assets, start, end):
    """Annual expected simple returns (from cc means) and annual covariance matrix of cc returns."""
    rr = r.loc[start:end, assets].dropna()
    mu = np.exp(12 * rr.mean()) - 1
    cov = 12 * rr.cov()
    return mu, cov

mu, COV = opt_inputs(ASSETS, START, END)
RF_ANNUAL = rf.loc[END, "tb3ms"] / 100           # T-Bill rate at the decision date
print(mu.round(4)); print(f"Annual risk-free rate at {END}: {RF_ANNUAL:.2%}")

AAPL    0.3830
JPM     0.1666
WMT     0.1082
XOM     0.5065
GLD     0.0736
dtype: float64
Annual risk-free rate at 2023-03-31: 4.64%

The expected return and the volatility of a portfolio with weights w are (in matrix form):

E[R_P]=w^{T}\mu \qquad\qquad \sigma_P=\sqrt{w^{T}\Sigma w}

def port_perf(w, mu, cov, rf=0.0):
    """Expected return, volatility and Sharpe ratio of a portfolio with weights w."""
    w = np.asarray(w)
    ret = w @ mu.values
    vol = np.sqrt(w @ cov.values @ w)
    return ret, vol, (ret - rf) / vol

w_eq = np.repeat(1 / len(ASSETS), len(ASSETS))
print("Equally weighted: return = %.4f, volatility = %.4f, Sharpe = %.4f" % port_perf(w_eq, mu, COV, RF_ANNUAL))

Equally weighted: return = 0.2476, volatility = 0.1649, Sharpe = 1.2204

2.1 From random portfolios to exact optimization

In the Algorithms course you simulated random portfolios. Let’s do it again, to see the cloud of feasible portfolios:

rng = np.random.default_rng(seed=42)
W = rng.random((5000, len(ASSETS)))
W = W / W.sum(axis=1, keepdims=True)                         # each row adds up to 1
cloud_ret = W @ mu.values
cloud_vol = np.sqrt(np.einsum("ij,jk,ik->i", W, COV.values, W))   # w' COV w for every row
best = np.argmax((cloud_ret - RF_ANNUAL) / cloud_vol)
print("Best random portfolio (Sharpe):", dict(zip(ASSETS, W[best].round(3))))

Best random portfolio (Sharpe): {'AAPL': np.float64(0.382), 'JPM': np.float64(0.022), 'WMT': np.float64(0.028), 'XOM': np.float64(0.491), 'GLD': np.float64(0.077)}

np.einsum calculates w^T\Sigma w for the 5,000 rows at once (it replaces the loop or list comprehension of the Algorithms course). The best random portfolio is only an approximation: with 5,000 tries we never hit the exact optimum, and with 20 assets random search becomes hopeless. We need optimization.

3 The Global Minimum Variance (GMV) portfolio

The GMV portfolio has the lowest possible volatility. It solves:

\min_{w}\; w^{T}\Sigma w \qquad \text{subject to} \qquad w^{T}\mathbf{1}=1

If short sales are allowed (weights can be negative), this problem has an exact solution (using the method of Lagrange multipliers):

w_{GMV}=\frac{\Sigma^{-1}\mathbf{1}}{\mathbf{1}^{T}\Sigma^{-1}\mathbf{1}}

Read it as an algorithm: (1) invert the covariance matrix, (2) multiply it by a vector of ones (this gives a “raw” weight for each asset: higher for assets with low variance and low covariances), (3) divide by the sum of the raw weights so that they add up to 1.

ones = np.ones(len(ASSETS))
Sinv = np.linalg.inv(COV.values)
raw = Sinv @ ones
w_gmv = raw / (ones @ raw)
ret_gmv, vol_gmv, sh_gmv = port_perf(w_gmv, mu, COV, RF_ANNUAL)
print(pd.Series(w_gmv, index=ASSETS).round(4))
print(f"GMV: return = {ret_gmv:.4f}, volatility = {vol_gmv:.4f}; minimum random volatility = {cloud_vol.min():.4f}")

AAPL   -0.0706
JPM     0.0913
WMT     0.2477
XOM     0.1349
GLD     0.5967
dtype: float64
GMV: return = 0.1272, volatility = 0.1229; minimum random volatility = 0.1249

The GMV volatility is lower than the volatility of every random portfolio. Notice that the GMV portfolio does not use the expected returns \mu: only the covariance matrix.

4 The tangency (maximum Sharpe ratio) portfolio

The Sharpe ratio is the premium return per unit of risk: \frac{E[R_P]-R_f}{\sigma_P}. The portfolio with the highest Sharpe ratio is called the tangency (or optimal) portfolio. With short sales allowed:

w_{T}=\frac{\Sigma^{-1}(\mu-R_f\mathbf{1})}{\mathbf{1}^{T}\Sigma^{-1}(\mu-R_f\mathbf{1})}

It is the same algorithm as the GMV, but the vector of ones is replaced by the vector of premium returns: assets with a higher premium return get a higher raw weight.

raw_t = Sinv @ (mu.values - RF_ANNUAL)
w_tan = raw_t / (ones @ raw_t)
ret_tan, vol_tan, sh_tan = port_perf(w_tan, mu, COV, RF_ANNUAL)
print(pd.Series(w_tan, index=ASSETS).round(4))
print(f"Tangency: return = {ret_tan:.4f}, volatility = {vol_tan:.4f}, Sharpe = {sh_tan:.4f}")

AAPL    0.5827
JPM    -0.5032
WMT    -0.1962
XOM     0.9266
GLD     0.1901
dtype: float64
Tangency: return = 0.6015, volatility = 0.3221, Sharpe = 1.7232

5 The efficient frontier and the Capital Market Line

A beautiful result of portfolio theory (the two-fund theorem): every portfolio on the efficient frontier is a combination of any two frontier portfolios. Then, we can draw the whole frontier by combining the GMV and the tangency portfolios: w(\lambda)=\lambda\, w_T+(1-\lambda)\,w_{GMV}.

Adding the risk-free asset, the best possible portfolios lie on a straight line from R_f through the tangency portfolio: the Capital Market Line (CML):

E[R_P]=R_f+Sharpe_T\times\sigma_P

An investor puts a fraction y in the tangency portfolio and 1-y in T-Bills: y<1 for conservative investors; y>1 means borrowing at the risk-free rate (leverage).

lams = np.linspace(-0.5, 2.0, 200)
front = [port_perf(l * w_tan + (1 - l) * w_gmv, mu, COV) for l in lams]
f_ret, f_vol = np.array([f[0] for f in front]), np.array([f[1] for f in front])
x = np.linspace(0, f_vol.max(), 50)

plt.figure(figsize=(9, 6))
plt.scatter(cloud_vol, cloud_ret, s=4, alpha=0.3, label="random portfolios (long only)")
plt.plot(f_vol, f_ret, color="black", label="efficient frontier (short sales allowed)")
plt.plot(x, RF_ANNUAL + sh_tan * x, color="green", ls="--", label="Capital Market Line")
plt.scatter(vol_gmv, ret_gmv, marker="D", s=90, color="red", label="GMV", zorder=5)
plt.scatter(vol_tan, ret_tan, marker="*", s=250, color="orange", label="Tangency", zorder=5)
asset_vol = np.sqrt(np.diag(COV.values))
plt.scatter(asset_vol, mu.values, color="purple", zorder=5)
for i, a in enumerate(ASSETS):
    plt.annotate(a, (asset_vol[i], mu.values[i]), fontsize=8)
plt.xlabel("annual volatility"); plt.ylabel("annual expected return")
plt.legend(fontsize=8); plt.grid(alpha=0.3); plt.show()

6 Without short sales: a numerical optimizer

Most investors (and most funds) cannot short-sell. With the constraint w_i\geq0 there is no formula, so we use a numerical optimizer: scipy.optimize.minimize. It works like the bisection and gradient descent algorithms of the Algorithms course: start from an initial guess and improve it step by step until the objective cannot be improved. We must give it:

• the objective function to minimize (the variance, or the negative Sharpe ratio, because the function minimizes);
• an initial guess (equal weights);
• bounds for each weight: (0, 1);
• constraints: the weights add up to 1.

def optimize_portfolio(mu, cov, objective="gmv", rf=0.0, long_only=True):
    """GMV or max-Sharpe portfolio with scipy (SLSQP)."""
    n = len(mu)
    if objective == "gmv":
        fun = lambda w: w @ cov.values @ w
    else:
        fun = lambda w: -port_perf(w, mu, cov, rf)[2]
    bounds = [(0, 1)] * n if long_only else [(-2, 2)] * n
    cons = [{"type": "eq", "fun": lambda w: w.sum() - 1}]
    res = minimize(fun, np.repeat(1 / n, n), method="SLSQP", bounds=bounds, constraints=cons,
                   options={"ftol": 1e-12, "maxiter": 1000})
    assert res.success, res.message
    return pd.Series(res.x, index=mu.index)

# Check: without the long-only restriction, the optimizer must find the closed-form GMV
w_check = optimize_portfolio(mu, COV, "gmv", long_only=False)
assert np.allclose(w_check.values, w_gmv, atol=1e-4)
w_gmv_lo = optimize_portfolio(mu, COV, "gmv")
w_gmv_lo.round(4)

AAPL    0.0000
JPM     0.0858
WMT     0.2048
XOM     0.1286
GLD     0.5808
dtype: float64

The assert is our verification: when the optimizer solves the same problem as the formula, both must give the same answer.

NotePyPortfolioOpt (optional)

The library PyPortfolioOpt does the same calculations (!pip install PyPortfolioOpt in Colab). If you use it, be careful: max_sharpe() uses a risk-free rate of 2% unless you give it risk_free_rate=, and portfolio_performance() reports the portfolio of the last optimization you ran on that object. Always verify its results with your own code.

from pypfopt import EfficientFrontier
ef = EfficientFrontier(mu, COV, weight_bounds=(0, 1))
ef.max_sharpe(risk_free_rate=RF_ANNUAL)
print(ef.clean_weights()); ef.portfolio_performance(risk_free_rate=RF_ANNUAL, verbose=True)

7 CHALLENGE 1: Optimal portfolios for your 5 stocks

In Workshop 2 (Challenge 2, step 5) you selected 5 stocks: from the 20 stocks of the ROA screen, the 10 with the highest alpha, and from them the 5 with the lowest beta.

1. Using the window April 2020 – March 2023, calculate mu and COV for your 5 stocks with opt_inputs, and the correlation matrix of their monthly cc returns. Which pair is the most correlated?
2. Calculate the GMV and the tangency portfolios with short sales allowed (closed-form formulas), and their expected return, volatility and Sharpe ratio (use the T-Bill rate of March 2023).
3. Calculate the GMV and the tangency portfolios without short sales (optimize_portfolio). Compare the weights with the previous ones. Which version would you recommend to a client, and why?
4. Interpret the long-only GMV portfolio: weights, expected return and volatility, in words that a client understands.

# Write your code here

8 CHALLENGE 2: Efficient frontier and CML of your 5 stocks

1. Draw, in the same plot: 5,000 random long-only portfolios of your 5 stocks, the efficient frontier with short sales (two-fund theorem), the CML, the 5 stocks, the two GMV portfolios and the two tangency portfolios.
2. A client wants an expected volatility of 15% per year. Using the CML, what fraction y must be invested in the tangency portfolio (short sales allowed) and what is the expected return? Is the client borrowing or lending at the risk-free rate?

# Write your code here

9 Estimation risk: why optimized portfolios are fragile

The inputs \mu and \Sigma are estimates from 36 months of data, and the expected returns are especially noisy: the standard error of a mean annual return is about \sigma/\sqrt{T}, which for a stock with 30% volatility and 3 years of data is about 17 percentage points! The optimizer takes the estimates as if they were true, and it puts large weights on the assets whose expected returns are overestimated. This is why many practitioners prefer the GMV (which does not use \mu) or simply equal weights.

Let’s see it with our 5 example assets: we add a small random error (standard deviation of 3 percentage points, much smaller than the real estimation error) to the expected returns, 20 times, and recompute the tangency portfolio with and without short sales. The table shows the lowest and the highest weight that each asset received in the 20 experiments:

rng = np.random.default_rng(seed=1)
w_short, w_long = [], []
for k in range(20):
    mu_noisy = mu + rng.normal(0, 0.03, len(mu))
    raw_k = Sinv @ (mu_noisy.values - RF_ANNUAL)
    w_short.append(pd.Series(raw_k / (ones @ raw_k), index=ASSETS))          # short sales allowed
    w_long.append(optimize_portfolio(mu_noisy, COV, "sharpe", rf=RF_ANNUAL))   # long only
w_short, w_long = pd.DataFrame(w_short), pd.DataFrame(w_long)
pd.DataFrame({"short allowed: min": w_short.min(), "short allowed: max": w_short.max(),
              "long only: min": w_long.min(), "long only: max": w_long.max()}).round(2)

	short allowed: min	short allowed: max	long only: min	long only: max
AAPL	0.38	1.06	0.24	0.44
JPM	-1.09	-0.29	0.00	0.00
WMT	-0.69	0.23	0.00	0.00
XOM	0.68	1.72	0.45	0.61
GLD	-0.92	0.41	0.00	0.24

Errors of only 3 percentage points in the expected returns move the “optimal” weights a lot, especially when short sales are allowed. The long-only constraint makes the portfolio more stable, but the weights still change by many percentage points. The GMV does not change at all in this experiment, because it does not use \mu.

10 Backtesting

A backtest answers: if I had invested with this rule in the past, what would have happened? To be honest, a backtest must be out of sample: the portfolio is estimated with data before the decision date, and it is evaluated with data after the decision date. Evaluating a portfolio in the same window used to estimate it is meaningless: of course the optimized portfolio looks great with the data that was used to optimize it.

Two ways to hold a portfolio:

• Buy and hold: buy the initial weights and do nothing. The weights drift: the assets that go up gain weight.
• Monthly rebalancing: every month, trade back to the initial weights. The monthly portfolio return is simply R[assets] @ w.

def backtest(w, assets, start, end, rebalance=False):
    """Monthly simple returns of a portfolio held from start to end."""
    Rw = R.loc[start:end, assets].fillna(0)            # a missing return is treated as 0
    w = np.asarray(w)
    if rebalance:
        return Rw @ w
    growth = (1 + Rw).cumprod()                         # value of $1 in each asset
    value = growth @ w                                  # value of the portfolio
    return value.pct_change().fillna(value.iloc[0] - 1)

def summary_stats(port_R, rf_annual=0.0):
    """HPR, annualized return, annualized volatility and Sharpe ratio from monthly simple returns."""
    n = len(port_R)
    hpr = (1 + port_R).prod() - 1
    ann_ret = (1 + hpr) ** (12 / n) - 1
    ann_vol = np.log(1 + port_R).std() * np.sqrt(12)
    return pd.Series({"HPR": hpr, "annual return": ann_ret, "annual volatility": ann_vol,
                      "Sharpe": (ann_ret - rf_annual) / ann_vol, "months": n})

OOS_START, OOS_END = "2023-04-30", "2024-03-31"          # the 12 months after the decision date
bt = pd.DataFrame({
    "GMV long-only (B&H)": backtest(w_gmv_lo, ASSETS, OOS_START, OOS_END),
    "Equal weights (B&H)": backtest(w_eq, ASSETS, OOS_START, OOS_END),
    "SPY": R.loc[OOS_START:OOS_END, "SPY"]})
rf_oos = (1 + rf.loc[OOS_START:OOS_END, "rf"]).prod() ** (12 / len(bt)) - 1
bt.apply(summary_stats, rf_annual=rf_oos).T.round(4)

	HPR	annual return	annual volatility	Sharpe	months
GMV long-only (B&H)	0.1833	0.1833	0.0935	1.3918	12.0
Equal weights (B&H)	0.2178	0.2178	0.1004	1.6397	12.0
SPY	0.2962	0.2962	0.1332	1.8240	12.0

((1 + bt).cumprod()).plot(figsize=(9, 4.5), title="Value of $1 invested on March 31, 2023 (out of sample)")
plt.grid(alpha=0.3); plt.show()

summary_stats calculates the HPR with the product of the growth factors, \prod(1+R_t)-1, and then the annualized (geometric) return (1+HPR)^{12/n}-1. The volatility is calculated with cc returns, as in the Algorithms course.

11 CHALLENGE 3: Out-of-sample backtest of your portfolio

1. Backtest from April 2023 to March 2024 (buy and hold) the following portfolios of your 5 stocks: the long-only GMV, the long-only tangency, and the equally weighted portfolio. Add SPY as the benchmark.
2. Report the HPR, annual return, volatility and Sharpe ratio of each one, and plot the value of $1.
3. Repeat the backtest of the long-only GMV with monthly rebalancing. Is the difference large? Why?
4. Did the optimized portfolios beat the equally weighted portfolio and the market? Was the realized volatility of the GMV close to its expected volatility? Write a short paragraph for your investment committee.
5. Extension: extend the backtest to the last month of the data. Does your conclusion change?

# Write your code here

12 CHALLENGE 4: Bug hunt

4.1 Weights and returns:

w = np.array([0.4, 0.3, 0.2, 0.1, 0.0])           # weights for ["XOM", "AAPL", "WMT", "JPM", "GLD"]
port = R[["AAPL", "GLD", "JPM", "WMT", "XOM"]] @ w

4.2 Portfolio returns with log returns:

port_r = r[ASSETS] @ w_eq
hpr = np.exp(port_r.sum()) - 1

4.3 Inputs:

mu = np.exp(12 * r.loc[START:END, ASSETS].mean()) - 1
COV = r.loc[START:END, ASSETS].cov()

4.4 Backtest:

w = optimize_portfolio(mu, COV, "gmv")
bt = backtest(w, ASSETS, START, END)               # START, END = estimation window

4.5 PyPortfolioOpt:

ef = EfficientFrontier(mu, COV, weight_bounds=(0, 1))
w = ef.min_volatility()
print("Tangency portfolio:", ef.portfolio_performance(risk_free_rate=0.045))

Write your answers here

13 Reflection: beyond the rational investor

Portfolio theory assumes that investors are rational and that the expected returns and covariances are known. Neither is true. Based on what you learned in this workshop (estimation risk, out-of-sample tests) and on your own view of the market (trends, regulations, fears, technologies), propose in a short paragraph how you would modify the optimized portfolio for a real client, and why.

14 DataCamp

Take the chapter “Portfolio Investing” (Chapter 2 of Introduction to Portfolio Risk Management in Python): efficient frontier, GMV and maximum-Sharpe portfolios.

15 W3 submission

Submit your .ipynb file in Canvas and answer the Canvas Quiz W3 before the first class of next week.