How to sell stocks wisely - the code for Almgren-Chriss optimal execution

Imagine you have a big bag of marbles that you want to sell at your school fair. But if you try to sell them all at once, your friends might not have enough pocket money to buy them, and you might have to sell them for a very low price. This is what grown-ups call "market impact."

Now, imagine if you sold just a few marbles each day of the fair. That way, you don't have too many marbles at once, and your friends can buy them without spending all their money. This is what the Almgren-Chriss model helps with — it's like a smart plan that tells you how many marbles to sell each day so you can get the most money by the end of the fair without making it hard for your friends to buy them.

The Framework

The Almgren-Chriss framework, developed by Robert Almgren and Neil Chriss, is a widely used model in the field of optimal execution for portfolio transactions. The aim of the model is to minimize a combination of volatility risk and transaction costs arising from permanent and temporary market impact. In this post, we'll explore how to implement the Almgren-Chriss model in Python for single and multiple assets.

IF YOU WANT ALL CODE, AND COMPLETE POST PLEASE SUBSCRIBE AT https://quantjourney.substack.com/

The Almgren-Chriss model is based on a few key assumptions (see also below with equation):

1. Temporary market impact: The execution of a trade creates a temporary impact on the price of the asset, which decays over time according to a specific function.
2. Permanent market impact: A portion of the market impact caused by the trade is permanent and persists even after the trade is completed.
3. Risk aversion: The trader is risk-averse and seeks to minimize the expected cost of execution, accounting for both the trading cost and the risk associated with the unexecuted portion of the order.

class AlmgrenChriss_Single:
    """
    Implementation of a simplified Almgren-Chriss optimal execution model for single asset transactions.

    This class provides methods to calculate the optimal trading trajectory for liquidating a position
    in a single asset over a specified time horizon. It aims to minimize the expected cost of execution
    taking into account market impact and risk aversion.

    The model simplifies market impact into temporary and permanent components and assumes risk aversion
    is linearly related to the square of the execution speed.

    Parameters are adjusted to include a simplified market impact model without explicit reference to
    trading days, hours, or specific asset details.
    """

    def __init__(self, params):
        """
        Initialize the Almgren-Chriss optimal execution model with specified parameters.

        Args:
            params (dict): Dictionary containing model parameters:
                - lambda (float): The temporary market impact parameter.
                - sigma (float): The volatility of the asset.
                - epsilon (float): The fixed cost of trading.
                - eta (float): The permanent market impact parameter.
                - gamma (float): The risk aversion parameter.
                - tau (float): The time interval between trades.
        """
        self._lambda = params['lambda']
        self._sigma = params['sigma']
        self._epsilon = params['epsilon']
        self._eta = params['eta']
        self._gamma = params['gamma']
        self._tau = params['tau']

        self._eta_tilda = self._eta - 0.5*self._gamma*self._tau

        assert self._eta_tilda > 0, "Model assumption violated: Adjusted permanent impact must be positive."

				# Parameter noted as (1) below
        self._kappa_tilda_squared = (self._lambda * self._sigma ** 2) / self._eta_tilda
        self._kappa = np.arccosh(0.5 * (self._kappa_tilda_squared * self._tau ** 2) + 1) / self._tau

    def trajectory(self, X, T):
        """
        Calculate the optimal liquidation trajectory for the given position over the specified time horizon.

        Args:
            X (int): The total number of shares to be liquidated.
            T (int): The total time horizon for liquidation in discrete intervals.

        Returns:
            np.array: An array of integers representing the number of shares to hold after each time interval.
        """

				# Parameter noted as (2) below
        holds = [int(np.sinh(self._kappa * (T - t)) / np.sinh(self._kappa * T) * X) for t in range(T)]
        holds.append(0)  # Ensure complete liquidation by the end of the horizon
        return np.array(holds)

    def strategy(self, X, T):
        """
        Calculate the optimal liquidation trade list based on the Almgren-Chriss model.

        Args:
            X (int): The total number of shares to be liquidated.
            T (int): The total time horizon for liquidation in discrete intervals.

        Returns:
            np.array: An array representing the number of shares to execute at each time interval.
        """
        return np.diff(self.trajectory(X, T))        

with the code we can see the optimal outcome over time, e.g.

or with different lambda:

or 3D for both gamma and eta.

Model for Multiple-Security Portfolios

We also created class for multiple-security portfolios which has m securities, and thus is a bit more complex (below only trajectory function):

def trajectory(self, X, T, general=False):
        # Optimal Liquidation Trajectory
        trajectory = []
        
        if not general:
            # Diagonal Model 
            if not is_diagonal(self._H_tilda): raise ValueError
            
            z0 = self._U.T@sqrtm(self._H_tilda)@X
            for t in range(T+1):
                # Parameter noted as (3) below
                z = np.sinh(self._kappas*(T - t))/np.sinh(self._kappas*T)*z0
                # Parameter noted as (4) below
                x = np.floor(np.linalg.inv(sqrtm(self._H_tilda))@self._U@z)
                trajectory.append(x)
        else:
            # General Model (only supported for 2 dimensions - did not have time to generalize)
            if self._dims != 2: raise ValueError
            
            # Transformation
            # Parameter noted as (5) below
            y0 = sqrtm(self._H_tilda)@X
            
            # Build Linear System
            rhs = np.zeros(2 * (T + 1))  # Initialize the right-hand side as zero for the linear equation system

            # Initial Conditions: Set the first two conditions based on the initial transformed positions y0
            rhs[0] = y0[0]  # The initial condition for the first asset
            rhs[1] = y0[1]  # The initial condition for the second asset

            # Initialize system matrix with initial conditions
            init1 = np.zeros(2 * (T + 1))  # Coefficients for the initial condition of the first asset
            init1[0] = 1  # Setting the first element to 1 as it's the coefficient of y_0 for the first asset
            init2 = np.zeros(2 * (T + 1))  # Coefficients for the initial condition of the second asset
            init2[1] = 1  # Setting the second element to 1 as it's the coefficient of y_0 for the second asset

            system = [init1, init2]  # Start the system matrix with initial condition rows

            # Define system dynamics based on the discretized differential equation
            for k in range(0, T - 1):
                # Coefficients derived from discretizing the differential equation governing the system
                a = (1 / self._tau ** 2)
                b = 1 / (2 * self._tau)
                c = -2 / (self._tau ** 2)
                l = self._lambda
                A = self._A
                B = self._B
            
                # Coefficients for the system equations at time k for both assets
                equation1_coeff = [a - b * B[0, 0], -b * B[0, 1], c - l * A[0, 0], -l * A[0, 1], a + b * B[0, 0], b * B[0, 1]]
                equation2_coeff = [-b * B[1, 0], a - b * B[1, 1], -l * A[1, 0], c - l * A[1, 1], b * B[1, 0], a + b * B[1, 1]]

                # Append equations for time k into the system matrix with correct placement
                system.append(np.pad(equation1_coeff, (2 * k, 2 * (T - k - 2))))
                system.append(np.pad(equation2_coeff, (2 * k, 2 * (T - k - 2))))

            # Final Conditions: Ensure that no shares are held at the end of the time horizon
            final1 = np.zeros(2 * (T + 1))
            final1[-2] = 1  # y_{1,N} = 0 for the first asset
            system.append(final1)

            final2 = np.zeros(2 * (T + 1))
            final2[-1] = 1  # y_{2,N} = 0 for the second asset
            system.append(final2)

            # Solve Linear System: Solve for the trajectory y_k at all time steps k
            solve = np.linalg.solve(np.array(system), rhs)
            y = solve.reshape(T + 1, 2)  # Reshape the solution to match the trajectory shape

            # Undo Transformation: Convert the transformed positions back to the original position space
            for yk in y:
                trajectory.append(np.linalg.inv(sqrtm(self._H_tilda)) @ yk)

        return np.array(trajectory).T        

giving some plots as:

and the plots for different \lambda (risk-averse impact) - Interestingly, for values of \lambda around the order of magnitude 10^{?8}, the optimal strategy involves selling shares early and may even involve short selling. This is due to the correlation between asset prices. As noted in the paper, some components of the velocity could be non-monotonic in time. Although the model may not be entirely accurate due to the changing sign of the cost associated with the bid-ask spread, this effect is minimal.

and for different values of Gamma (permanent market impact):

Summary

The Almgren-Chriss model enables the selection of an optimal liquidation strategy based on the standard mean-variance optimization of the total trading cost. However, for risk-averse traders, the size of their trades diminishes over time.