I start new series about Quantum computing which is aimed to sharing knowledge regarding to this new and excited technology to the public. Hope it helpful for people who interest in this field, So let’s get started.

This is 1st article about Quantum computing (QC) applied in financial domain for portfolio optimization using Modern Portfolio Theory (MPT) technique in order to help investors minimize portfolio’s risk and maximize return in long-term. It will be separated in to 5 parts as follow

• What is Modern Portfolio Theory (MPT)?
• Classical approach
• Python implementation of stochastic and deterministic method
• Quantum algorithm approach
• Python implementation of quantum algorithm using Qiskit library

What is Modern Portfolio Theory (MPT)?

Basically the modern portfolio theory (MPT) is a practical method for selecting investments in order to maximize overall returns within an acceptable risk (Asset’s volatility).

This mathematical framework is used to build a portfolio of investments that maximize the amount of expected return for the given risk level.

Harry Markowitz an American economist who develop this theory in his paper “Portfolio Selection,” which was published in the Journal of Finance in 1952. He was later awarded a Nobel Prize for his work on modern portfolio theory. Original paper can be found here https://www.jstor.org/stable/2975974

However every mathematical model has it’s own assumptions and limitations. Here are some limitations of this model.

• The limitations of Markowitz model include over reliance on historical data, irrelevant assumptions, and the use of mean-variance instead of potential risks.
• Markowitz‘s assumptions become irrelevant; this is especially the case with volatile markets.

Classical approach

There are 2 ways in order to build portfolio in classical approach which is have it’s own pros and cons.

• Stochastic method using Monte Carlo simulation technique which will find approximate global maximum in term of sharpe ratio also approximately global minimum in term of volatility as well.
• Deterministic method using SLSQP algorithm for example.

Deterministic method can find global maximum or minimum solution for small size dimension (number of assets) but have downside in term of computational cost for larger portfolio (e.g 1000++ dimensions or assets). On the other hand Stochastic method like Monte Carlo technique can find approximated global solution in acceptable computation cost. This is about combinatorial related problem which is the main challenging of portfolio optimization task which is where Quantum computing algorithm comes in.

Next section we will implement the portfolio optimization for a specified group of crypto currencies asset with python using both classical (Stochastic and Deterministic) method as well as quantum algorithm method.

Python implementation of stochastic and deterministic method

First we will download crypto currencies historical data using pandas datareader.

from pandas_datareader import data as web
from datetime import datetime
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
import pandas as pd
import scipy.optimize as opt
import seaborn as sns

import yfinance as yf
yf.pdr_override()

symbols = ["ETH-USD", "BNB-USD", "ADA-USD",'LINK-USD','LTC-USD']

start_date = '2020-12-30'
end_date = '2023-12-30'

data = web.DataReader(symbols,
                    start_date,
                    end_date)

adj_close = data['Adj Close']

Then declare return calculation functions as follow.

def calc_returns(price_data, resample=None, ret_type="arithmatic"):
    if ret_type=="arithmatic":
        ret = price_data.pct_change().dropna()
    elif ret_type=="log":
        ret = np.log(price_data/price_data.shift()).dropna()
    else:
        raise ValueError("ret_type: return calculation type is not valid. use \"arithmatic\" or \"log\"")

    if resample != None:
        if ret_type=="arithmatic":
            ret = ret.resample(resample).apply(lambda df: (df+1).cumprod(axis=0).iloc[-1]) - 1
        elif ret_type=="log":
            ret = ret.resample(resample).apply(lambda df: df.sum(axis=0))
    return(ret)

def calc_returns_stats(returns):
    mean_returns = returns.mean()
    cov_matrix = returns.cov()
    return(mean_returns, cov_matrix)

def portfolio(weights, mean_returns, cov_matrix):

    portfolio_return = np.dot(weights.reshape(1,-1), mean_returns.values.reshape(-1,1))
    portfolio_var = np.dot(np.dot(weights.reshape(1,-1), cov_matrix.values), weights.reshape(-1,1))
    portfolio_std = np.sqrt(portfolio_var)

    return(np.squeeze(portfolio_return),np.squeeze(portfolio_var),np.squeeze(portfolio_std))

1. Stochastic approach using Monte-Carlo simulation method

epoch = 100000

porfolio_var_list = []
porfolio_ret_list = []
w_list =[]

max_sharpe = 0
max_sharpe_var = None
max_sharpe_ret = None
max_sharpe_w = None

daily_ret = calc_returns(adj_close, resample=None, ret_type="log")
mean_returns, cov_matrix = calc_returns_stats(daily_ret)

for i in range(1,epoch+1):
    rand_weights = np.random.random(len(symbols))
    rand_weights = rand_weights/np.sum(rand_weights)

    porfolio_ret, porfolio_var, portfolio_std = portfolio(rand_weights, mean_returns, cov_matrix)

    porfolio_ret = porfolio_ret * 252
    porfolio_var = porfolio_var * 252
    portfolio_std = portfolio_std * (252**0.5)

    sharpe = (porfolio_ret/(porfolio_var**0.5)).item()
    if sharpe > max_sharpe:
        max_sharpe = sharpe
        max_sharpe_var = porfolio_var.item()
        max_sharpe_ret = porfolio_ret.item()
        max_sharpe_w = rand_weights

    porfolio_var_list.append(porfolio_var)
    porfolio_ret_list.append(porfolio_ret)
    w_list.append(rand_weights)

Calculate returns and volatility.

stat = np.hstack([np.array(porfolio_ret_list).reshape(-1,1),
                 np.array(porfolio_var_list).reshape(-1,1),
                 (np.array(porfolio_ret_list)/np.sqrt(np.array(porfolio_var_list))).reshape(-1,1)]
                 )
stat=pd.DataFrame(stat, columns=['Return','Variance','Sharpe ratio'])
stat.head(3)

Then calculate minimum risk portfolio for crypto assets.

min_vol = stat[stat['Variance'] == stat['Variance'].min()]
min_vol_index = min_vol.index.item()
min_vol_weights = w_list[min_vol_index]
print("Portfolio with minimum volatility:\n")
print(f"Annual Sharpe Ratio: {round(stat['Sharpe ratio'][min_vol_index],2)} | Annual Return: % {round(stat['Return'][min_vol_index]*100,2)} | Annual Volatility: % {round(stat['Variance'][min_vol_index]*100,2)}\n")
for index,symbol in enumerate(symbols):
    print(f'{symbol}:\t% {round(min_vol_weights[index]*100,2)}')

Now let’s plot efficient frontier graph of portfolio.

daily_ret_var = daily_ret.var()
plt.figure(figsize=(10,8))
plt.scatter(porfolio_var_list,porfolio_ret_list,c=stat['Sharpe ratio'], alpha=0.8, cmap='magma')
for sym in symbols:
    plt.scatter(daily_ret_var.loc[sym]*252, mean_returns.loc[sym]*252, marker='x', s=100, label=sym)

plt.scatter([max_sharpe_var], [max_sharpe_ret], marker='*', s=250, label='max sharpe porfolio', c='blue')
plt.scatter(min_vol['Variance'].values, min_vol['Return'].values, marker='*', s=250, label='minimum volatility porfolio', c='green')

plt.xlabel('Variance')
plt.ylabel('Return')
plt.title('Portfolio Optimization using Monte Carlo Simulation method')
plt.legend(loc='upper right')

cmap = mpl.cm.magma
norm = mpl.colors.Normalize(vmin=stat['Sharpe ratio'].min(), vmax=stat['Sharpe ratio'].max())
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap),label='Sharpe Ratio',);

2. Deterministic approach using SLSQP method in Scipy library

First we declare functions for optimizer.

• func_negative_sharpe_ratio: Calculate negative sharpe ratio which will minimized by optimizer.
• func_portfolio_variance: Calculate portfolio’s annualize volatility.

def func_negative_sharpe_ratio(weights, mean_returns, cov_matrix, risk_free_rate=0):
    portfolio_return, portfolio_var, portfolio_std = portfolio(weights, mean_returns, cov_matrix)
    sr = ((portfolio_return - risk_free_rate)/portfolio_std) * (252**0.5)
    return(-sr)

def func_portfolio_variance(weights, mean_returns, cov_matrix):
    portfolio_return, portfolio_var, portfolio_std = portfolio(weights, mean_returns, cov_matrix)
    return(portfolio_var*252)

Calculate covariance matrix and mean returns from daily returns.

daily_ret = calc_returns(adj_close, resample=None, ret_type="log")
mean_returns, cov_matrix = calc_returns_stats(daily_ret)

Create optimization function to find portfolio’s weights in order to maximize sharpe ratio by try to minimize portfolio’s negative sharpe ratio, Then run this function.

def func_optimize_sharpe_ratio(mean_returns, cov_matrix, risk_free_rate=0, w_bounds=(0,1)):

    init_guess = np.array([1/len(mean_returns) for _ in range(len(mean_returns))])
    args = (mean_returns, cov_matrix, risk_free_rate)
    constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1})
    result = opt.minimize(fun=func_negative_sharpe_ratio,
                          x0=init_guess,
                          args=args,
                          method='SLSQP',
                          bounds=tuple(w_bounds for _ in range(len(mean_returns))),
                          constraints=constraints,
                          )
    
    if result['success']:
        print(result['message'])
        opt_sharpe = - result['fun']
        opt_weights = result['x']
        opt_return, opt_variance, opt_std = portfolio(opt_weights, mean_returns, cov_matrix)
        return(opt_sharpe, opt_weights, opt_return.item()*252, opt_variance.item()*252, opt_std.item()*(252**0.5))
    else:
        print("Optimization operation was not succesfull!")
        print(result['message'])
        return(None)

opt_sharpe, opt_weights, opt_return, opt_variance, opt_std = func_optimize_sharpe_ratio(mean_returns, cov_matrix, risk_free_rate=0, w_bounds=(0,1))

Optimization terminated successfully

Check the result.

print("Portfolio with maximum sharpe ratio:\n")
print(f"Annual Sharpe Ratio: {round(opt_sharpe,2)} | Annual Return: % {round(opt_return*100,2)} | Annual Volatility: % {round(opt_variance*100,2)}\n")
for index,symbol in enumerate(symbols):
    print(f'{symbol}:\t% {round(opt_weights[index]*100,2)}')

Next we will minimize portfolio’s volatility.

def func_minimize_portfolio_variance(mean_returns, cov_matrix, w_bounds=(0,1)):
    
    init_guess = np.array([1/len(mean_returns) for _ in range(len(mean_returns))])
    args = (mean_returns, cov_matrix)
    constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1})
    result = opt.minimize(fun=func_portfolio_variance,
                          x0=init_guess,
                          args=args,
                          method='SLSQP',
                          bounds=tuple(w_bounds for _ in range(len(mean_returns))),
                          constraints=constraints,
                          )
    
    if result['success']:
        print(result['message'])
        min_var = result['fun']
        min_var_weights = result['x']
        min_var_return, min_var_variance, min_var_std = portfolio(min_var_weights, mean_returns, cov_matrix)
        min_var_sharpe = (min_var_return/min_var_std)*(252**0.5)
        return(min_var_sharpe, min_var_weights, min_var_return.item()*252, min_var_variance.item()*252, min_var_std.item()*(252**0.5))
    else:
        print("Optimization operation was not succesfull!")
        print(result['message'])
        return(None)

min_var_sharpe, min_var_weights, min_var_return, min_var_variance, min_var_std = func_minimize_portfolio_variance(mean_returns, cov_matrix, w_bounds=(0,1))

Optimization terminated successfully

Let’s check the result.

print("Portfolio with minimum volatility (variance):\n")
print(f"Annual Sharpe Ratio: {round(min_var_sharpe,3)} | Annual Return: % {round(min_var_return*100,2)} | Annual Volatility: % {round(min_var_variance*100,3)}\n")
for index,symbol in enumerate(symbols):
    print(f'{symbol}:\t% {round(min_var_weights[index]*100,2)}')

We also able to find minimum variance portfolio according to specified return we want.

def func_calc_portfolio_return(weights, mean_returns, cov_matrix):
    portfolio_return, portfolio_var, portfolio_std = portfolio(weights, mean_returns, cov_matrix)
    return(portfolio_return.item()*252)

def func_efficient_portfolio(mean_returns, cov_matrix, target_return, w_bounds=(0,1)):
    """retuens the portfolio weights with minimum variance for a specific level of expected portfolio return"""    

    init_guess = np.array([1/len(mean_returns) for _ in range(len(mean_returns))])
    args = (mean_returns, cov_matrix)
    constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x)-1},
                   {'type': 'eq', 'fun': lambda x: 252*np.squeeze(np.dot(x.reshape(1,-1),mean_returns.values.reshape(-1,1))) - target_return})
    result = opt.minimize(fun=func_portfolio_variance,
                          x0=init_guess,
                          args=args,
                          method='SLSQP',
                          bounds=tuple(w_bounds for _ in range(len(mean_returns))),
                          constraints=constraints,
                          )
    if not result['success']:
        print(result['message'])
    efficient_variance = result['fun']
    efficient_weights = result['x']
    efficient_return, _ , efficient_std = portfolio(efficient_weights, mean_returns, cov_matrix)
    efficient_sahrpe = (efficient_return/efficient_return)*(252**0.5)
    return(efficient_sahrpe, efficient_weights, efficient_return.item()*252, efficient_variance, efficient_std.item()*(252**0.5))

# We want 40% return for example.
expected_return = 0.4
efficient_sharpe, efficient_weights, efficient_return, efficient_variance, efficient_std = func_efficient_portfolio(mean_returns,
                                                                                                                cov_matrix,
                                                                                                                target_return=expected_return,
                                                                                                                w_bounds=(0,1))

Result here.

print("Portfolio for 40% return level:\n")
print(f"Annual Sharpe Ratio: {round(efficient_sharpe ,3)} | Annual Return: % {round(efficient_return*100,2)} | Annual Volatility: % {round(efficient_variance*100,3)}\n")
for index,symbol in enumerate(symbols):
    print(f'{symbol}:\t% {round(efficient_weights[index]*100,2)}')

Plot efficient frontier.

daily_ret_var = daily_ret.var()
plt.figure(figsize=(10,8))

target_rets = np.arange(0.23,0.35,0.002)
efficient_vars = np.array([func_efficient_portfolio(mean_returns,cov_matrix,target_return=x ,w_bounds=(0,1))[3] for x in target_rets])
plt.plot(efficient_vars,target_rets,marker='.',label="efficient frontier")

plt.scatter(porfolio_var_list,porfolio_ret_list,c=stat['Sharpe ratio'], alpha=0.8, cmap='magma')
for sym in symbols:
    plt.scatter(daily_ret_var.loc[sym]*252, mean_returns.loc[sym]*252, marker='x', s=100, label=sym)

plt.scatter([opt_variance], [opt_return], marker='*', s=250, label='max sharpe p', c='blue')
plt.scatter(min_var_variance, min_var_return, marker='*', s=250, label='min vol p', c='green')

plt.xlabel("Variance")
plt.ylabel("Return")
plt.title("Portfolio's Efficient Frontier")
plt.legend(loc="upper right")

cmap = mpl.cm.magma
norm = mpl.colors.Normalize(vmin=stat['Sharpe ratio'].min(), vmax=stat['Sharpe ratio'].max())
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap),label='Sharpe Ratio',);

Quantum algorithm approach

Since classical optimization algorithms running on classical Von-Neumann machine face computation problem when try to optimize portfolio which have highly dimensions variables (financial assets in our case) Then Quantum algorithms steps in, by using unique characteristics of quantum mechanics such as Superposition, Entanglement and Interference which will be explained in depth in my next articles, which is in theory able to solve large scale optimization problem more efficient.

Python implementation of quantum algorithm using Qiskit library

This section we will solve portfolio optimization problem using quantum algorithms implemented with IBM’s Qiskit library, so let’s get start coding.

This time we will download financial data by using Qiskit Finance package. First import necessary libraries.

%matplotlib inline
from qiskit_finance import QiskitFinanceError
from qiskit_finance.data_providers import *
import datetime
import matplotlib.pyplot as plt
from pandas.plotting import register_matplotlib_converters

register_matplotlib_converters()


from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver, QAOA, SamplingVQE
from qiskit.algorithms.optimizers import COBYLA
from qiskit.circuit.library import TwoLocal
from qiskit.result import QuasiDistribution
from qiskit_aer.primitives import Sampler
from qiskit_finance.applications.optimization import PortfolioOptimization
from qiskit_finance.data_providers import RandomDataProvider
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit.utils import algorithm_globals
import numpy as np

Then download crypto’s historical data as follow.

try:
    data = YahooDataProvider(
        tickers=["ETH-USD", "BNB-USD", "ADA-USD",'LINK-USD','LTC-USD'],
        start=datetime.datetime(2021, 1, 1),
        end=datetime.datetime(2021, 12, 31),
    )
    data.run()
    for (cnt, s) in enumerate(data._tickers):
        plt.plot(data._data[cnt], label=s)
    plt.legend(loc="upper center", bbox_to_anchor=(0.5, 1.1), ncol=3)
    plt.xticks(rotation=90)
    plt.show()
except QiskitFinanceError as ex:
    data = None
    print(ex)

Now we will take a look of means as well as similarity and covariance matrix from crypto assets.

means = data.get_mean_vector()
print("Means:")
print(means)

rho = data.get_similarity_matrix()
print("A time-series similarity measure:")
print(rho)
plt.title("A time-series similarity measure:")
plt.imshow(rho)
plt.colorbar()
plt.show()

cov = data.get_covariance_matrix()
print("A covariance matrix:")
print(cov)
plt.title("A covariance matrix:")
plt.imshow(cov)
plt.colorbar()
plt.show()

Means: [2.77586987e+03 3.77938485e+02 1.50031422e+00 2.66617162e+01 1.85989058e+02] A time-series similarity measure: [[1.00000000e+00 1.14567613e-06 9.90225958e-07 9.99288739e-07 1.06076416e-06] [1.14567613e-06 1.00000000e+00 7.29801853e-06 7.87357416e-06 1.43218006e-05] [9.90225958e-07 7.29801853e-06 1.00000000e+00 1.09185201e-04 1.48911673e-05] [9.99288739e-07 7.87357416e-06 1.09185201e-04 1.00000000e+00 1.72428205e-05] [1.06076416e-06 1.43218006e-05 1.48911673e-05 1.72428205e-05 1.00000000e+00]]

A covariance matrix: [[1.04826889e+06 1.51016435e+05 4.57055885e+02 2.19497519e+03 1.60955598e+04] [1.51016435e+05 2.85843976e+04 6.99916238e+01 5.82798973e+02 4.44786591e+03] [4.57055885e+02 6.99916238e+01 3.80566812e-01 1.39906053e+00 6.86065529e+00] [2.19497519e+03 5.82798973e+02 1.39906053e+00 4.94479249e+01 3.13792028e+02] [1.60955598e+04 4.44786591e+03 6.86065529e+00 3.13792028e+02 2.32601929e+03]]

Calculate mean and covariance matrix.

mu = data.get_period_return_mean_vector()
sigma = data.get_period_return_covariance_matrix()

# plot sigma
plt.imshow(sigma, interpolation="nearest")
plt.colorbar()
plt.show()

Using qiskit’s built-in deterministic algorithm as a classical reference

num_assets = 4
q = 0.5  # set risk factor
budget = num_assets // 2  # set budget
penalty = num_assets  # set parameter to scale the budget penalty term

portfolio = PortfolioOptimization(
    expected_returns=mu, covariances=sigma, risk_factor=q, budget=budget
)
qp = portfolio.to_quadratic_program()
qp

def print_result(result):
    selection = result.x
    value = result.fval
    print("Optimal: selection {}, value {:.4f}".format(selection, value))

    eigenstate = result.min_eigen_solver_result.eigenstate
    probabilities = (
        eigenstate.binary_probabilities()
        if isinstance(eigenstate, QuasiDistribution)
        else {k: np.abs(v) ** 2 for k, v in eigenstate.to_dict().items()}
    )
    print("\n----------------- Full result ---------------------")
    print("selection\tvalue\t\tprobability")
    print("---------------------------------------------------")
    probabilities = sorted(probabilities.items(), key=lambda x: x[1], reverse=True)

    for k, v in probabilities:
        x = np.array([int(i) for i in list(reversed(k))])
        value = portfolio.to_quadratic_program().objective.evaluate(x)
        print("%10s\t%.4f\t\t%.4f" % (x, value, v))

exact_mes = NumPyMinimumEigensolver()
exact_eigensolver = MinimumEigenOptimizer(exact_mes)

result = exact_eigensolver.solve(qp)

print_result(result)

Global optimal selection using deterministic method is [0,1,1,0,0]

Solution using Sampling VQE

We can now use the Sampling Variational Quantum Eigensolver (SamplingVQE) to solve the problem.

We will specify the optimizer and variational form to be used.

algorithm_globals.random_seed = 1234

cobyla = COBYLA()
cobyla.set_options(maxiter=500)
ry = TwoLocal(num_assets, "ry", "cz", reps=3, entanglement="full")
vqe_mes = SamplingVQE(sampler=Sampler(), ansatz=ry, optimizer=cobyla)
vqe = MinimumEigenOptimizer(vqe_mes)
result = vqe.solve(qp)

print_result(result)

VQE algorithm able to find optimal solution as well.

Solution using QAOA

We also show here a result using the Quantum Approximate Optimization Algorithm (QAOA).

This is another variational algorithm and it uses an internal variational form that is created based on the problem.

algorithm_globals.random_seed = 1234

cobyla = COBYLA()
cobyla.set_options(maxiter=250)
qaoa_mes = QAOA(sampler=Sampler(), optimizer=cobyla, reps=3)
qaoa = MinimumEigenOptimizer(qaoa_mes)
result = qaoa.solve(qp)

print_result(result)

As expected QAOA algorithm also able to find optimal solution.

Quantum computing is an exciting domain which is contain a lot more details. Next articles will be deep down to basic foundation of quantum mechanics as well as basic quantum circuits.

I would like to give thanks to valuable articles and resources listed below.

– https://plotly.com/python/v3/ipython-notebooks/markowitz-portfolio-optimization/

– https://www.jstor.org/stable/2975974

– https://medium.datadriveninvestor.com/portfolio-optimization-with-python-using-scipy-optimize-monte-carlo-method-a5b4e89e0548

– https://scipy.github.io/devdocs/tutorial/optimize.html#sequential-least-squares-programming-slsqp-algorithm-method-slsqp

– https://qiskit.org/

– https://www.sciencedirect.com/science/article/pii/S0370157322003118

– https://arxiv.org/pdf/1812.01041.pdf

See you next time.