# Portfolio optimization using genetic algorithm

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?

• @Emma I don't know if I can write LaTeX here, but the problem is: maximize Sum(i=1, n) w_i * r_i subjected to w^T * Cov * w <= max_risk, Sum(i=1, n) w_i = 1, w_i >= 0, where w is the weight vector, r_i is the return for the i-th stock, Cov is the covariance matrix between stocks and max_risk is the maximum risk the client is willing to tolerate. – Heathcliff Nov 27 '19 at 1:23
• I barely know any Python, so it won't risk posting an answer, but I have to say this is probably some of the most readable Python I've ever seen. About the algorithm itself though... your population size (num_solutions, 5000) seems way too high. I've found personally that in general, sizes from 20 to 100 work best, or about 10 times your number of dimensions, but I don't know how many dimensions you have here. – cliesens Nov 27 '19 at 2:03
• Some other things you could try : use other selection methods (again, personally, I've found roulette selection to not be the most efficient in most cases) such as tournament selection, or some custom ones (keep the best n% for crossover,...). They shouldn't be too hard to implement, and they will often provide the most noticeable changes in performance. Finally, you could also try making your mutation rate vary over time, although I would suggest trying this only after you've noticed actual improvements in your algorithm. – cliesens Nov 27 '19 at 2:10
• @cliesens You are not strictly required to know the programming language as long as you can give feedback with regard to the algorithm. – AlexV Nov 27 '19 at 8:02
• The current question title, which states your concerns about the code, is too general to be useful here. Please edit to the site standard, which is for the title to simply state the task accomplished by the code. Please see How to get the best value out of Code Review: Asking Questions for guidance on writing good question titles. – Toby Speight Nov 27 '19 at 11:28

• 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