Updated code and Notebook to work on Python 3

[adalseno]

Hi, I have updated the code and created a Jupyter notebook to work on Python 3 (the original code was in Python 2).
You can use them to update your repo.

#!/usr/bin/env python
# On 20121129
# Critical Line Algorithm
# by MLdP <[email protected]>
import numpy as np


# ---------------------------------------------------------------
# ---------------------------------------------------------------
class CLA:
    def __init__(self, mean, covar, lB, uB):
        # Initialize the class
        self.mean = mean
        self.covar = covar
        self.lB = lB
        self.uB = uB
        self.w = []  # solution
        self.l = []  # lambdas
        self.g = []  # gammas
        self.f = []  # free weights

    # ---------------------------------------------------------------
    def solve(self):
        # Compute the turning points,free sets and weights
        f, w = self.initAlgo()
        self.w.append(np.copy(w))  # store solution
        self.l.append(-np.inf)
        self.g.append(None)
        self.f.append(f[:])
        while True:
            # 1) case a): Bound one free weight
            l_in = -np.inf
            if len(f) > 1:
                covarF, covarFB, meanF, wB = self.getMatrices(f)
                covarF_inv = np.linalg.inv(covarF)
                j = 0
                for i in f:
                    l, bi = self.computeLambda(
                        covarF_inv, covarFB, meanF, wB, j, [self.lB[i], self.uB[i]]
                    )
                    if l > l_in:
                        l_in, i_in, bi_in = l, i, bi
                    j += 1
            # 2) case b): Free one bounded weight
            l_out = -np.inf
            if len(f) < self.mean.shape[0]:
                b = self.getB(f)
                for i in b:
                    covarF, covarFB, meanF, wB = self.getMatrices(f + [i])
                    covarF_inv = np.linalg.inv(covarF)
                    l, bi = self.computeLambda(
                        covarF_inv,
                        covarFB,
                        meanF,
                        wB,
                        meanF.shape[0] - 1,
                        self.w[-1][i],
                    )
                    if (self.l[-1] == -np.inf or l < self.l[-1]) and l > l_out:
                        l_out, i_out = l, i
            # 3) decide lambda
            if (l_in == -np.inf or l_in < 0) and (l_out == -np.inf or l_out < 0):
                break
            if l_in > l_out:
                # If it's the first elment and it's Nan (-np.inf) replace it
                if len(self.l) == 1 and self.l[0] == -np.inf:
                    self.l[0] = l_in
                else:
                    self.l.append(l_in)
                f.remove(i_in)
                w[i_in] = bi_in  # set value at the correct boundary
            else:
                # If it's the first elment and it's Nan (-np.inf) replace it
                if len(self.l) == 1 and self.l[0] == -np.inf:
                    self.l[0] = l_out
                else:
                    self.l.append(l_out)
                f.append(i_out)
            # 4) compute solution vector
            covarF, covarFB, meanF, wB = self.getMatrices(f)
            covarF_inv = np.linalg.inv(covarF)
            wF, g = self.computeW(covarF_inv, covarFB, meanF, wB)
            for i in range(len(f)):
                w[f[i]] = wF[i]
            self.w.append(np.copy(w))  # store solution
            self.g.append(g)
            self.f.append(f[:])
            if len(f) == self.mean.shape[0]:
                # 5) minimum variance solution
                wF, g = self.computeW(covarF_inv, covarFB, np.zeros(meanF.shape), wB)
                for i in range(len(f)):
                    w[f[i]] = wF[i]
                self.w.append(np.copy(w))  # store solution
                self.g.append(g)
                self.f.append(f[:])
        # Remove first element from weights (it's a duplicate)
        self.w = self.w[1:]
        # If we miss the last lambda insert it
        if len(self.l) < len(self.w):
            self.l += [0.0]

    # ---------------------------------------------------------------
    def initAlgo(self):
        # Initialize the algo
        # 1) Form structured array
        a = np.zeros((self.mean.shape[0]), dtype=[("id", int), ("mu", float)])
        b = [self.mean[i][0] for i in range(self.mean.shape[0])]  # dump array into list
        a[:] = list(zip(range(self.mean.shape[0]), b))  # fill structured array
        # 2) Sort structured array
        b = np.sort(a, order="mu")
        # 3) First free weight
        i, w = b.shape[0], np.copy(self.lB)
        while sum(w) < 1:
            i -= 1
            w[b[i][0]] = self.uB[b[i][0]]
        w[b[i][0]] += 1 - sum(w)
        return [b[i][0]], w

    # ---------------------------------------------------------------
    def computeBi(self, c, bi):
        if c > 0:
            bi = bi[1]
        if c < 0:
            bi = bi[0]
        return bi

    # ---------------------------------------------------------------
    def computeW(self, covarF_inv, covarFB, meanF, wB):
        # 1) compute gamma
        onesF = np.ones(meanF.shape)
        g1 = np.dot(np.dot(onesF.T, covarF_inv), meanF)
        g2 = np.dot(np.dot(onesF.T, covarF_inv), onesF)
        if wB is None:
            g, w1 = float(-self.l[-1] * g1 / g2 + 1 / g2), 0
        else:
            onesB = np.ones(wB.shape)
            g3 = np.dot(onesB.T, wB)
            g4 = np.dot(covarF_inv, covarFB)
            w1 = np.dot(g4, wB)
            g4 = np.dot(onesF.T, w1)
            g = float(-self.l[-1] * g1 / g2 + (1 - g3 + g4) / g2)
        # 2) compute weights
        w2 = np.dot(covarF_inv, onesF)
        w3 = np.dot(covarF_inv, meanF)
        return -w1 + g * w2 + self.l[-1] * w3, g

    # ---------------------------------------------------------------
    def computeLambda(self, covarF_inv, covarFB, meanF, wB, i, bi):
        # 1) C
        onesF = np.ones(meanF.shape)
        c1 = np.dot(np.dot(onesF.T, covarF_inv), onesF)
        c2 = np.dot(covarF_inv, meanF)
        c3 = np.dot(np.dot(onesF.T, covarF_inv), meanF)
        c4 = np.dot(covarF_inv, onesF)
        c = -c1 * c2[i] + c3 * c4[i]
        if c == 0:
            return
        # 2) bi
        if type(bi) == list:
            bi = self.computeBi(c, bi)
        # 3) Lambda
        if wB is None:
            # All free assets
            return float((c4[i] - c1 * bi) / c), bi
        else:
            onesB = np.ones(wB.shape)
            l1 = np.dot(onesB.T, wB)
            l2 = np.dot(covarF_inv, covarFB)
            l3 = np.dot(l2, wB)
            l2 = np.dot(onesF.T, l3)
            return float(((1 - l1 + l2) * c4[i] - c1 * (bi + l3[i])) / c), bi

    # ---------------------------------------------------------------
    def getMatrices(self, f):
        # Slice covarF,covarFB,covarB,meanF,meanB,wF,wB
        covarF = self.reduceMatrix(self.covar, f, f)
        meanF = self.reduceMatrix(self.mean, f, [0])
        b = self.getB(f)
        covarFB = self.reduceMatrix(self.covar, f, b)
        wB = self.reduceMatrix(self.w[-1], b, [0])
        return covarF, covarFB, meanF, wB

    # ---------------------------------------------------------------
    def getB(self, f):
        return self.diffLists(np.arange(self.mean.shape[0]), f)

    # ---------------------------------------------------------------
    def diffLists(self, list1, list2):
        return list(set(list1) - set(list2))

    # ---------------------------------------------------------------
    def reduceMatrix(self, matrix, listX, listY):
        # Reduce a matrix to the provided list of rows and columns
        if len(listX) == 0 or len(listY) == 0:
            return
        matrix_ = matrix[:, listY[0] : listY[0] + 1]
        for i in listY[1:]:
            a = matrix[:, i : i + 1]
            matrix_ = np.append(matrix_, a, 1)
        matrix__ = matrix_[listX[0] : listX[0] + 1, :]
        for i in listX[1:]:
            a = matrix_[i : i + 1, :]
            matrix__ = np.append(matrix__, a, 0)
        return matrix__

    # ---------------------------------------------------------------
    def getMinVar(self):
        # Get the minimum variance solution
        var = []
        for w in self.w:
            a = np.dot(np.dot(w.T, self.covar), w)
            var.append(a)
        return min(var) ** 0.5, self.w[var.index(min(var))]

    # ---------------------------------------------------------------
    def getMaxSR(self):
        # Get the max Sharpe ratio portfolio
        # 1) Compute the local max SR portfolio between any two neighbor turning points
        w_sr, sr = [], []
        for i in range(len(self.w) - 1):
            w0 = np.copy(self.w[i])
            w1 = np.copy(self.w[i + 1])
            kargs = {"minimum": False, "args": (w0, w1)}
            a, b = self.goldenSection(self.evalSR, 0, 1, **kargs)
            w_sr.append(a * w0 + (1 - a) * w1)
            sr.append(b)
        return max(sr), w_sr[sr.index(max(sr))]

    # ---------------------------------------------------------------
    def evalSR(self, a, w0, w1):
        # Evaluate SR of the portfolio within the convex combination
        w = a * w0 + (1 - a) * w1
        b = np.dot(w.ravel(), self.mean)[0]
        c = np.sqrt(np.dot(np.dot(w.ravel(), self.covar), w.ravel()))
        return b / c

    # ---------------------------------------------------------------
    def goldenSection(self, obj, a, b, **kargs):
        # Golden section method. Maximum if kargs['minimum']==False is passed
        # from math import log,ceil
        tol, sign, args = 1.0e-9, 1, None
        if "minimum" in kargs and kargs["minimum"] == False:
            sign = -1
        if "args" in kargs:
            args = kargs["args"]
        numIter = int(np.ceil(-2.078087 * np.log(tol / np.abs(b - a))))
        r = 0.618033989
        c = 1.0 - r
        # Initialize
        x1 = r * a + c * b
        x2 = c * a + r * b
        f1 = sign * obj(x1, *args)
        f2 = sign * obj(x2, *args)
        # Loop
        for i in range(numIter):
            if f1 > f2:
                a = x1
                x1 = x2
                f1 = f2
                x2 = c * a + r * b
                f2 = sign * obj(x2, *args)
            else:
                b = x2
                x2 = x1
                f2 = f1
                x1 = r * a + c * b
                f1 = sign * obj(x1, *args)
        if f1 < f2:
            return x1, sign * f1
        else:
            return x2, sign * f2

    # ---------------------------------------------------------------
    def efFrontier(self, points):
        # Get the efficient frontier
        mu, sigma, weights = [], [], []
        a = np.linspace(0, 1, points // len(self.w))[
            :-1
        ]  # remove the 1, to avoid duplications
        b = np.arange(len(self.w) - 1)
        for i in b:
            w0, w1 = self.w[i], self.w[i + 1]
            if i == b[-1]:
                a = np.linspace(
                    0, 1, points // len(self.w)
                )  # include the 1 in the last iteration
            for j in a:
                w = w1 * j + (1 - j) * w0
                weights.append(np.copy(w))
                mu.append(np.dot(w.T, self.mean)[0, 0])
                sigma.append(np.dot(np.dot(w.T, self.covar), w)[0, 0] ** 0.5)
        return mu, sigma, weights


# ---------------------------------------------------------------
# ---------------------------------------------------------------

And the notebook:

CLA in Python 3

Original code in Python 2, together with the example dataset, is here: https://github.com/mdengler/cla

import CLA
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import Markdown, display

# Python >= 3.9 or you might have issues with type-hints

Helper functions

def compute_return_risk(
    df: pd.DataFrame, mean: np.ndarray, covar: np.ndarray
) -> tuple[pd.Series, pd.Series]:
    """Compute Return and Risk (volatility) series for the CLW weights

    Args:
        df (pd.DataFrame): Dataframe with the weights
        mean (np.ndarray): Array of mean returns
        covar (np.ndarray): Array of covariances

    Returns:
        tuple[pd.Series, pd.Series]: series of returns and risk
    """
    p_ret = []
    risk = []
    for i in df.index:
        p_ret.append(np.dot(df.loc[i, :].values, mean)[0])
        risk.append(
            np.sqrt(np.dot(df.loc[i, :].values, np.dot(covar, df.loc[i, :].values)))
        )
    return pd.Series(p_ret, index=df.index, name="Return"), pd.Series(
        risk, index=df.index, name="Risk"
    )

def create_weight_table(
    cla: CLA, mean: np.ndarray, covar: np.ndarray, var_names: list
) -> pd.DataFrame:
    """Create the table with Return, Risk, Lambda and Weights

    Args:
        cla (CLA): The CLA object with the solution
        mean (np.ndarray): Array with the return means
        covar (np.ndarray): Array of covariances
        var_names (list): List of names for the columns

    Returns:
        pd.DataFrame: The weights table
    """
    weights = pd.DataFrame(
        np.hstack(cla.w).T, columns=var_names, index=np.arange(1, len(cla.w) + 1)
    )
    port_return, port_risk = compute_return_risk(weights, mean, covar)
    port_lambda = pd.Series(cla.l, index=weights.index, name="Lambda")
    return pd.concat([port_return, port_risk, port_lambda, weights], axis=1)

def display_results(ret_mean: np.ndarray, cla: CLA) -> None:
    """Plot the CLA results and dispaly the maximum Sharpe portfolio
    and the minimum variance one

    Args:
        ret_mean (np.ndarray): Array or return means
        cla (CLA): CLA object with the solution
    """
    plot_results()
    print()
    # 5) Get Maximum Sharpe ratio portfolio
    sr, w_sr = cla.getMaxSR()
    display(
        Markdown(
            f"### Portfolio volatility: {np.sqrt(np.dot(np.dot(w_sr.ravel(),cla.covar),w_sr.ravel())):.2%}, Sharpe ratio: {sr:.2f}"
        )
    )
    display(Markdown("### Weights (rounded to 4 decimal places):"))
    display(
        pd.Series(
            np.round(w_sr.ravel(), 4),
            name="sr_weights",
            index=np.arange(1, len(w_sr) + 1),
        )
    )
    print()
    # 6) Get Minimum Variance portfolio
    mv, w_mv = cla.getMinVar()
    mu_v = np.dot(w_mv.ravel(), ret_mean)
    sr_v = np.round(mu_v / mv, 2).ravel()[0]
    display(
        Markdown(
            f"### Portfolio minimum volatility: {mv.ravel()[0]:.2%}, Sharpe ratio: {sr_v:.2f}"
        )
    )
    display(Markdown("### Weights (rounded to 4 decimal places):"))
    display(
        pd.Series(
            np.round(w_mv.ravel(), 4),
            name="mv_weights",
            index=np.arange(1, len(w_sr) + 1),
        )
    )
    print()
    return

def plot_results() -> None:
    """Plots the results of CLA"""
    mu, sigma, weights = cla.efFrontier(100)
    mu = np.array(mu)
    sigma = np.array(sigma)
    x_low = max(sigma.min() - 0.05, 0.0)
    fig, ax = plt.subplots(1, 2, figsize=(20, 7))  # one row, one column, first plot
    ax[0].plot(sigma, mu, color="blue")
    ax[0].set_xlabel("Risk")
    ax[0].set_ylabel("Expected Excess Return", rotation=90)
    # ax[0].set_xticklabels(ax[0].get_xticks(), rotation='vertical')
    ax[0].set_xlim(x_low, 1.0)
    ax[0].set_title("CLA-derived Efficient Frontier")
    ax[1].plot(sigma, np.round(np.array(mu) / np.array(sigma), 2), color="blue")
    ax[1].set_xlabel("Risk")
    ax[1].set_ylabel("Sharpe ratio", rotation=90)
    # ax[1].set_xticklabels(ax[0].get_xticks(), rotation='vertical')
    ax[1].set_xlim(x_low, 1.0)
    ax[1].set_title("CLA-derived Sharpe Ratio function")
    plt.suptitle("CLA results", fontsize=16)
    plt.show()

    return

Import dataset

df = pd.read_csv("CLA_Data.csv")
df

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style>

	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10
0	1.175000	1.190000	0.396000	1.120000	0.346000	0.679000	0.089000	0.730000	0.481000	1.080000
1	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000
3	0.407552	0.031758	0.051839	0.056639	0.033023	0.008278	0.021659	0.013324	0.034348	0.022499
4	0.031758	0.906305	0.031364	0.026873	0.019172	0.009344	0.024950	0.007610	0.028749	0.013369
5	0.051839	0.031364	0.194909	0.044085	0.030068	0.013227	0.035260	0.011549	0.042756	0.020573
6	0.056639	0.026873	0.044085	0.195285	0.027773	0.005267	0.013758	0.007809	0.029142	0.016404
7	0.033023	0.019172	0.030068	0.027773	0.340591	0.007771	0.020678	0.007364	0.025427	0.012841
8	0.008278	0.009344	0.013227	0.005267	0.007771	0.159839	0.021056	0.005187	0.017237	0.007238
9	0.021659	0.024950	0.035260	0.013758	0.020678	0.021056	0.680567	0.013779	0.046270	0.019261
10	0.013324	0.007610	0.011549	0.007809	0.007364	0.005187	0.013779	0.955269	0.010655	0.007610
11	0.034348	0.028749	0.042756	0.029142	0.025427	0.017237	0.046270	0.010655	0.316816	0.018543
12	0.022499	0.013369	0.020573	0.016404	0.012841	0.007238	0.019261	0.007610	0.018543	0.110793

Create variables

var_names = df.columns.to_list()

var_names

['X1', 'X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8', 'X9', 'X10']

mean = df.iloc[0, :].values.reshape(-1, 1)

lB = df.iloc[1, :].values.reshape(-1, 1)
uB = df.iloc[2, :].values.reshape(-1, 1)

covar = df.iloc[3:, :].values

Compute CLA

# 3) Invoke object
cla = CLA.CLA(mean, covar, lB, uB)

cla.solve()

Display results

weight_df = create_weight_table(cla, mean, covar, var_names)

weight_df.style.format("{:.3f}")

<style type="text/css"> </style>

	Return	Risk	Lambda	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10
1	1.190	0.952	58.303	0.000	1.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
2	1.180	0.546	4.174	0.649	0.351	0.000	-0.000	0.000	0.000	0.000	0.000	0.000	0.000
3	1.160	0.417	1.946	0.434	0.231	0.000	0.335	0.000	0.000	0.000	0.000	0.000	-0.000
4	1.111	0.267	0.165	0.127	0.072	0.000	0.281	0.000	0.000	0.000	-0.000	0.000	0.520
5	1.108	0.265	0.147	0.123	0.070	0.000	0.279	0.000	0.000	0.000	0.006	0.000	0.521
6	1.022	0.230	0.056	0.087	0.050	0.000	0.224	0.000	0.174	0.000	0.030	0.000	0.435
7	1.015	0.228	0.052	0.085	0.049	0.000	0.220	-0.000	0.180	0.000	0.031	0.006	0.429
8	0.973	0.220	0.037	0.074	0.044	0.000	0.199	0.026	0.198	0.000	0.033	0.028	0.398
9	0.950	0.216	0.031	0.068	0.041	0.015	0.188	0.034	0.202	0.000	0.034	0.034	0.383
10	0.803	0.205	0.000	0.037	0.027	0.095	0.126	0.077	0.219	0.030	0.036	0.061	0.292

display_results(mean, cla)

Portfolio volatility: 22.74%, Sharpe ratio: 4.45

Weights (rounded to 4 decimal places):

1     0.0840
2     0.0489
3     0.0000
4     0.2183
5     0.0017
6     0.1812
7     0.0000
8     0.0312
9     0.0079
10    0.4269
Name: sr_weights, dtype: float64

Portfolio minimum volatility: 20.52%, Sharpe ratio: 3.91

Weights (rounded to 4 decimal places):

1     0.0370
2     0.0269
3     0.0949
4     0.1258
5     0.0767
6     0.2194
7     0.0300
8     0.0360
9     0.0613
10    0.2920
Name: mv_weights, dtype: float64

Thank you very much for your repo.

Kind regards,
Andrea Dalseno