Portfolio optimization using genetic algorithm

Question

I'm working on a (naïve) algorithm for portfolio optimization using GA. It takes a list of stocks, calculates its expected returns and the covariance between all of them and then it returns the portfolio weights that would produce the highest return of investment given a certain maximum risk the investor is willing to be exposed to.

The expected return of the portfolio is calculated multiplying the weight of the stock (a percentage) by that stock's expected return, while the risk is calculated using the matrix equation w^T * Cov * w, where w is the array of weights and Cov is the covariance matrix.

I'm having trouble improving performance with my algorithm. I've timed the execution of the different parts of the algorithm and found that the crossing is the part that takes the longest (around 5 seconds per iteration).

geneticAlgorithmPortfolioGenerator.py

This is the main class. The portfolio is generated by calling its generate_portfolio() method after instantiating the class with the maximum risk, the stock returns and the covariance matrix. It runs for generations generations unless it fails to improve after 100 generations, in which case it short-cuts the loop and returns. This has been the case every time I've executed it.

class Genetic_Algorithm_Portfolio:
    def __init__(self, max_risk, generations, returns, cov_matrix):
        '''
        Creates a genetic algorithm portfolio generator

        - population: list of Candidate objects
        - max_risk: risk constraint the portfolio must meet
        - generations: number of generations the algorithm will run to create a better solution
        '''
        
        self.max_risk = max_risk
        self.generations = generations
        self.returns = returns
        self.cov_matrix = cov_matrix

    def generate_portfolio(self):
        #print("Generating portfolio")
        # Map solutions so that the Candidate objects also have normalized fitness value and cummulative sum
        t0 = time.time()
        random_generator = Random_Portfolio_Generator(5000, self.max_risk, self.returns, self.cov_matrix)
        t1 = time.time()
        print("Time to generate random set: " + str(int(t1 - t0)) + " seconds")
        self.initial_population = random_generator.generate_solutions()[0]
        if len(self.initial_population) == 0:
            return None
        self.map_solutions()
        current_generation = self.initial_population
        initial_weight = [0]*len(self.initial_population[0].weights)
        best_solution = Candidate(
            initial_weight,
            portfolioUtilities.calculate_return(initial_weight, self.returns),
            portfolioUtilities.calculate_std(initial_weight, self.cov_matrix)
        )

        # With this initial population, run 100 generations
        improvements = 100
        for i in range(self.generations):
            t_g0 = time.time()
            if improvements == 0:
                print("No improvements after 100 generations")
                break

            # Next generation will be those selected by the selection method + the 10 best by expected return from the current generation
            t_s0 = time.time()
            next_gen_selected = self.next_generation(current_generation)
            current_generation.sort(key=lambda x: x.expected_return, reverse=True)
            next_gen_selected = next_gen_selected + current_generation[:10]
            t_s1 = time.time()
            print(f"Next generation selection: {str(int(t_s1 - t_s0))} seconds")

            # Pair parents
            t_p0 = time.time()
            parents = self.pair(next_gen_selected)
            t_p1 = time.time()
            print(f"Parents selection: {str(int(t_p1 - t_p0))} seconds")
            
            # Cross them and create the next gen
            next_gen = []
            t_c0 = time.time()
            for j in range(len(parents)):
                next_gen = next_gen + self.cross(parents[j])

            # Check if you found a better solution
            next_gen.sort(key=lambda x: x.expected_return, reverse=True)
            t_c1 = time.time()
            print(f"Cross: {str(int(t_c1 - t_c0))} seconds")
            if len(next_gen) == 0:
                continue

            next_gen_best_solution = next_gen[0]
            if (best_solution.expected_return < next_gen_best_solution.expected_return):
                best_solution = next_gen_best_solution
                improvements = 100
            else:
                improvements -= 1
            
            # Next generation is now current generation
            current_generation = next_gen
            t_g1 = time.time()
            print(f">>> Generation {i} time: {str(int(t_g1 - t_g0))} seconds")
        #print("Done!")
        return best_solution

Random solutions

I generate 5000 random solutions simply checking that all weights are above 0 and the risk is under the max risk passed as parameter:

lass Random_Portfolio_Generator:
    def __init__(self, n, max_risk, returns, cov_matrix):
        self.num_solutions = n
        self.max_risk = max_risk
        self.returns = returns
        self.cov_matrix = cov_matrix

    def generate_solutions(self):
        solutions = []
        best_solution = []
        best_return = float("-inf")
        best_std = 0
        seed = 42
        counter = 0
        attempts = 0
        while counter < self.num_solutions and attempts < 1000:
            random.seed(seed)
            seed += 1
            weights = self.generate_random_weights()
            solution = Candidate(
                weights,
                portfolioUtilities.calculate_return(weights, self.returns),
                portfolioUtilities.calculate_std(weights, self.cov_matrix)
            )
            if (solution.std <= self.max_risk):
                attempts = 0
                counter +=1
                solutions.append(solution)
                if (solution.expected_return > best_return):
                    best_return = solution.expected_return
                    best_solution = solution
                    best_std = solution.std
            else:
                attempts += 1
        return solutions, best_solution, best_return, best_std
    
    def generate_random_weights(self):
        n = len(self.returns)
        arr = [0] * n
        indexes = list(range(n))
        sum = 0;
        for i in range(n - 1):
            rand = random.randint(0,1000 - sum)
            index = random.choice(indexes)
            indexes.remove(index)
            sum += rand
            arr[index] = rand
        arr[indexes[0]] = 1000 - sum
        return np.array([x / 1000 for x in arr])

Selection

My selection method is encapsulated in the next_generation() method, which uses roulette selection by cumulative sum of the expected return of each portfolio. I also sort the current generation by expected return and add the 10 best chromosomes to the next generation.

    def next_generation(self, current):
        # Select chromosomes to pair
        cum_sums = [o.cum_sum for o in self.initial_population]
        selected = []
        for x in range(len(cum_sums)//2):
            selected.append(self.roulette(cum_sums, np.random.rand()))

            # In case the roulette picks an index that has already been picked
            while len(set(selected)) != len(selected):
                selected[x] = self.roulette(cum_sums, np.random.rand())
        
        return [self.initial_population[int(selected[x])] for x in range(len(self.initial_population)//2)]

    def roulette(self, cum_sum, chance):
        veriable = list(cum_sum.copy())
        veriable.append(chance)
        veriable = sorted(veriable)
        return veriable.index(chance)

Pairing

I pair the next generation by a "next best" approach:

    def pair(self, generation):
        # Pairs selected chromosomes by fitness
        return [generation[i:i + 2] for i in range(0, len(generation), 2)]

 Crossing

To cross the pairs, I take a middle point d = (x1 - x2) / 3 and add it and subtract it from each parent to create 4 new solutions. I found this method in a paper after finding that the traditional crossing methods ended up always creating invalid solutions (more risk than the max tolerated by the user). I filter those that don't fit that criteria neither:

    def cross(self, parents):
        x1 = parents[0].weights
        x2 = parents[1].weights
        d = (x1 - x2) / 3
        x3 = self.mutate(x1 + d)
        x4 = self.mutate(x1 - d)
        x5 = self.mutate(x2 + d)
        x6 = self.mutate(x2 - d)

        return sorted(self.judge_candidates([
            parents[0],
            parents[1],
            Candidate(
                x3,
                portfolioUtilities.calculate_return(x3, self.returns),
                portfolioUtilities.calculate_std(x3, self.cov_matrix)
            ),
            Candidate(
                x4,
                portfolioUtilities.calculate_return(x4, self.returns),
                portfolioUtilities.calculate_std(x4, self.cov_matrix)
            ),
            Candidate(
                x5,
                portfolioUtilities.calculate_return(x5, self.returns),
                portfolioUtilities.calculate_std(x5, self.cov_matrix)
            ),
            Candidate(
                x6,
                portfolioUtilities.calculate_return(x6, self.returns),
                portfolioUtilities.calculate_std(x6, self.cov_matrix)
            )
        ]), key=lambda x: x.expected_return, reverse=True)[:2]

    def judge_candidates(self, candidates):
        return list(filter(lambda x: x.std <= self.max_risk and all(i >= 0 for i in x.weights), candidates))

Mutation

Finally, I mutate two random genes in the chromosome with probability 2% in random.gauss(0, 0.001).

    def mutate(self, candidate):
        if (random.randint(0, 100) <= 2):
            gene_one = random.randint(0, len(candidate) - 1)
            gene_two = random.randint(0, len(candidate) - 1)
            deviation = random.gauss(0, 0.001)
            candidate[gene_one] += deviation
            candidate[gene_two] -= deviation
        return candidate

Utilities

These are the two utility classes I use in this class.

Candidate.py

import math

class Candidate:
    def __init__(self, weights, expected_return, std):
        self.weights = weights
        self.expected_return = expected_return
        self.std = math.sqrt(std)
        self.print = f"weights: {self.weights}, return = {self.expected_return}, std = {self.std}"

portfolioUtilities.py

import numpy as np

def calculate_return(weights, returns):
    return (np.asarray(weights) * returns).sum()

def calculate_std(weights, cov_matrix):
    return np.dot(np.dot(np.asarray(weights), cov_matrix),np.asarray(weights).transpose())

Any ideas or criticisms about what I may be doing wrong?

Answer 1

Few tips:

• try hashing fitness function, create dictionary of hashed individuals and their fitness values; fitness evaluation is often costly; no reason to evaluate same specimen twice
• don't use roulette selection, you can loose selective pressure, when there's little difference between fitnesses or converge too fast, when there's huge deviation; use tournament selection instead, it doesn't care about fitness distribution and maintains selective pressure well
• GA is rather well suited for parallel computation; exploit that fact; you got nice multiprocessing package in python's standard lib
• try to use hill climber to improve some of the individuals; there was a lot of research about hybrids: simulated annealing, tabu search, FIHC. It depends on your problem and encoding though.

Now most important stuff - vanilla GA is slow, because it searches through space like a blind man. Try to learn about linkage learning (discrete problems) and fitness landscape analysis (continous problems). You'll be able to crash your problem way much faster. If you'd like to know more about it, get some papers, please create separate thread on https://datascience.stackexchange.com/ or https://stats.stackexchange.com/. We can talk there.