Portfolio Optimisation: Markowitz in Python

Modern Portfolio Theory

Harry Markowitz published "Portfolio Selection" in 1952 and changed finance forever. His insight: diversification is mathematically provable. By combining assets that are not perfectly correlated, you can reduce risk without sacrificing return. The optimal portfolio depends on the expected returns, volatilities, and correlations of all assets.

The efficient frontier is the set of portfolios that offer the maximum expected return for each level of risk. Any portfolio below the frontier is suboptimal: you could get more return for the same risk, or less risk for the same return. Rational investors should only hold portfolios on the frontier.

Loading Data and Calculating Returns

We start with 5 large FTSE 100 stocks across different sectors: Shell (energy), AstraZeneca (pharma), HSBC (banking), Unilever (consumer goods), and Diageo (beverages). Sector diversification gives us meaningful optimisation results.

import numpy as np
import pandas as pd
import yfinance as yf
from scipy.optimize import minimize

# Download FTSE 100 stocks
tickers = ["SHEL.L", "AZN.L", "HSBA.L", "ULVR.L", "DGE.L"]
data = yf.download(tickers, start="2020-01-01", end="2023-12-31")["Close"]

# Calculate daily returns
returns = data.pct_change().dropna()

# Annualised expected returns and covariance
mu = returns.mean() * 252
cov_matrix = returns.cov() * 252

print("Expected Annual Returns:")
for t, r in zip(tickers, mu):
    print(f"  {t}: {r:.2%}")
print(f"
Covariance Matrix Shape: {cov_matrix.shape}")

Random Portfolio Simulation

A simple but effective approach: generate 10,000 random portfolio weight combinations and plot their risk-return profiles. This gives you a visual feel for the efficient frontier before doing formal optimisation. The upper-left boundary of the scatter cloud is the frontier.

def random_portfolios(n_portfolios, mu, cov_matrix, n_assets):
    """Generate random portfolio weights and calculate risk/return."""
    results = np.zeros((n_portfolios, 3))  # return, vol, sharpe
    weights_record = []

    for i in range(n_portfolios):
        # Random weights that sum to 1
        w = np.random.random(n_assets)
        w /= w.sum()
        weights_record.append(w)

        # Portfolio return
        port_return = np.dot(w, mu)

        # Portfolio volatility
        port_vol = np.sqrt(np.dot(w.T, np.dot(cov_matrix, w)))

        # Sharpe ratio (assuming 4% risk-free rate)
        sharpe = (port_return - 0.04) / port_vol

        results[i] = [port_return, port_vol, sharpe]

    return results, weights_record

# Generate 10,000 random portfolios
results, weights = random_portfolios(10000, mu.values, cov_matrix.values, len(tickers))

# Find max Sharpe portfolio
max_sharpe_idx = results[:, 2].argmax()
print(f"Max Sharpe Portfolio:")
print(f"  Return: {results[max_sharpe_idx, 0]:.2%}")
print(f"  Vol:    {results[max_sharpe_idx, 1]:.2%}")
print(f"  Sharpe: {results[max_sharpe_idx, 2]:.2f}")
for t, w in zip(tickers, weights[max_sharpe_idx]):
    print(f"  {t}: {w:.1%}")

Optimal Portfolios with scipy

For precise optimisation, we use scipy.optimize.minimize with the SLSQP algorithm. We find two special portfolios: the maximum Sharpe ratio portfolio (best risk-adjusted return) and the minimum variance portfolio (lowest possible risk).

def max_sharpe_portfolio(mu, cov_matrix, rf=0.04):
    """Find the portfolio that maximises the Sharpe ratio."""
    n = len(mu)

    def neg_sharpe(w):
        port_return = np.dot(w, mu)
        port_vol = np.sqrt(np.dot(w.T, np.dot(cov_matrix, w)))
        return -(port_return - rf) / port_vol

    # Constraints: weights sum to 1
    constraints = ({
        'type': 'eq',
        'fun': lambda w: np.sum(w) - 1
    })

    # Bounds: no shorting (0 to 1 for each weight)
    bounds = tuple((0, 1) for _ in range(n))

    # Initial guess: equal weight
    w0 = np.array([1/n] * n)

    result = minimize(neg_sharpe, w0, method='SLSQP',
                      bounds=bounds, constraints=constraints)

    return result.x

def min_variance_portfolio(cov_matrix):
    """Find the minimum variance portfolio."""
    n = cov_matrix.shape[0]

    def portfolio_vol(w):
        return np.sqrt(np.dot(w.T, np.dot(cov_matrix, w)))

    constraints = ({'type': 'eq', 'fun': lambda w: np.sum(w) - 1})
    bounds = tuple((0, 1) for _ in range(n))
    w0 = np.array([1/n] * n)

    result = minimize(portfolio_vol, w0, method='SLSQP',
                      bounds=bounds, constraints=constraints)

    return result.x

# Optimise
opt_sharpe = max_sharpe_portfolio(mu.values, cov_matrix.values)
opt_minvar = min_variance_portfolio(cov_matrix.values)

print("Max Sharpe Weights:")
for t, w in zip(tickers, opt_sharpe):
    print(f"  {t}: {w:.1%}")

print("
Min Variance Weights:")
for t, w in zip(tickers, opt_minvar):
    print(f"  {t}: {w:.1%}")

Tracing the Efficient Frontier

To plot the exact efficient frontier, we solve the minimum variance problem for a range of target returns. Each point on the frontier is an optimisation problem: minimise volatility subject to achieving at least a given return.

def efficient_frontier(mu, cov_matrix, n_points=50):
    """Trace out the efficient frontier."""
    n = len(mu)
    target_returns = np.linspace(mu.min(), mu.max(), n_points)
    efficient_vols = []

    for target in target_returns:
        def portfolio_vol(w):
            return np.sqrt(np.dot(w.T, np.dot(cov_matrix, w)))

        constraints = [
            {'type': 'eq', 'fun': lambda w: np.sum(w) - 1},
            {'type': 'eq', 'fun': lambda w: np.dot(w, mu) - target}
        ]
        bounds = tuple((0, 1) for _ in range(n))
        w0 = np.array([1/n] * n)

        result = minimize(portfolio_vol, w0, method='SLSQP',
                          bounds=bounds, constraints=constraints)

        if result.success:
            efficient_vols.append(result.fun)
        else:
            efficient_vols.append(None)

    return target_returns, efficient_vols

targets, vols = efficient_frontier(mu.values, cov_matrix.values)
print("Efficient Frontier Points:")
for t, v in zip(targets, vols):
    if v is not None:
        print(f"  Return: {t:.2%}  Vol: {v:.2%}")