Black-Litterman Portfolio Optimization with Noisy Intermediate-Scale Quantum Computers

BL model can be generally understood as a method to construct a new combined return vector $\mathbf{\mu_{BL}}$, which is a combination of market implied return and investors’ view, to replace the $\mu$ in equation 1, and do portfolio optimization. The market implied return $\mathbf{\Pi}$ is the return we should have if the market capitalization weight $\mathbf{w_{mkt}}$ is the optimal solution for a portfolio optimization, and is calculated via reverse optimization:

where $\gamma$ is the risk aversion coefficient as in equation (1) that depends on investor’s preference. For a general case, we choose $\gamma$ based on the ratio of excess return and covariance of the market index over the pass 10 years. Since $\mathbf{\Pi}$ is fully determined by the market, we can incorporate quantum technology into BL model via investor’s views. Investor’s views can be quantized by the following three terms :
$\mathbf{P}$ : a matrix that identifies the assets involved in the views.
$\mathbf{Q}$ : the view vector.
$\mathbf{\Omega}$ : a diagonal matrix representing the uncertainty of views.
For example, if we have a view about two assets A and B, say A will have a return of 0.05 with uncertainty 0.001, and B will have a return of 0.1 with uncertainty 0.0005. The corresponding matrices are written as

Note that there are also relative views such as return of A asset will out perform B by some amount, however, in this work we consider only a direct view of every asset (thus $\mathbf{P}$ is an identity matrix). In general, these views can be generate from any kind of resource, including financial analysts, news, or personal experience…etc. However, to take advantage of the quantum classification methods, we should treat the views quantitatively. According to [10], the uncertainty of a view about $k^{th}$ asset can be taken as proportional to the variance of the asset itself. We can thus write

where $\omega_{k}=\Omega_{kk}$, $\mathbf{p_{k}}$ is the $k^{th}$ row of $\mathbf{P}$, and $\tau$ is a constant that will not affect the final combined return vector if we use equation 4 to construct $\mathbf{\Omega}$[9]. For the view vector $\mathbf{Q}$, [10]also suggested that we can set the $k^{th}$ view in $\mathbf{Q}$ to be

where $\eta_{k}\in\{-2,-1,1,2\}$ denotes the view as ”very bearish”, ”bearish”, ”bullish”, and ”very bullish” respectively. Here, we further defined $\mathbf{\eta}$ for two binary classification models:

where $Y_{1}\in\{+1,-1\}$, $Y_{2}\in\{1,2\}$, and $s_{i}$ is the accurate rate of the testing data of the classifier for $Y_{i}$. Now, $\eta_{k}\in[-2,2]$ since we have multiplied it with a scale $s_{1}s_{2}$. This means that we will have a stronger view for the classifier with higher accuracy, and vice versa. After we constructed these matrices, the combined return vector is calculated with the result formula in the original paper [8],

By replacing $\mathbf{\mu}$ with $\mathbf{\mu_{BL}}$ in equation (1), the objective function now becomes

and the corresponding QUBO form is

Many of the previous quantum finance studies either applied simple return, i.e. $R_{i}=\frac{p_{i}-p_{i-1}}{p_{i-1}}$, by default or without giving specific form for return. However, in finance, log return is more commonly used due to it’s statistical properties (follows normal distribution) and the ability of capturing compound return (time additive). Log return is defined as

and when we add them along a time sequence $t=0$ to $t=T$, it becomes $r_{0T}=log(\frac{p_{T}}{p_{0}})$. The rate of return (RoR) during this period is $e^{r_{0T}}-1$. On the other hand, the sum of simple return will not capture the compound return. Thus, it is more preferable to adopt log return when we are optimizing the data for a buy and hold strategy.

The risk aversion coefficient $\gamma$ is another thing that is also not discussed in quantum finance studies. In the discrete optimization case, the budget constraint, $1^{T}x=B$, is needed otherwise the solutions are not comparable. However, when we construct a portfolio in real cases, $x$ should always be re-normalized to 1 after we obtain the solution, i.e. $x_{r}=\frac{1}{B}x$. Since the variance in the objective function (equation (8)) is quadratic and the return is linear, the re-normalization will result in equivalently solving for a risk aversion coefficient $\frac{\gamma}{B}$. For cases that $\gamma$s are chosen randomly, this may not matter, nevertheless, if $\gamma$ is estimated with historical data and has financial meaning, which is the case here, it should be carefully considered. We thus propose an effective $\gamma$:

where $\gamma$ can be estimated with, for example, the ratio of the risk premium to the variance of index over the pass 10 years, i.e.

where $\mu$ and $\sigma^{2}$ are the mean return and variance of the index, and $r_{f}$ is the mean risk free rate. Here we refer $r_{f}$ to the 13 week treasury bill (Yahoo Finance : $\hat{}IRX$). We use $\gamma_{eff}$ rather than $\gamma$ for solving equation (8), and when we re-normalize the optimal solution such that the sum of portfolio weights is 1, it is an optimal solution for the risk aversion coefficient $\gamma$ that we want.

In the framework of BL model, the market cap weighted index is assumed to be efficient at present. For a problem of asset pool that are large enough and are all included in $S\&P500$, we can refer the index to the $S\&P500$ index. In our case, however, we have to construct the corresponding market weighted index for our assets pool. Here, 12 individual stocks are considered, we determined the market cap by the data at closest time to the end time of the training period. Since the market cap weight is in unit of money, we first assume a total 100 unit of money (i.e. the initially normalized index level), and calculate the shares for each asset at the time of market cap data. The shares then holds as a constant, and price of the index at different time is determined by

where $P_{i}$ is the price of $i$ asset, and $S_{i}$ is the shares of $i$ asset.

The walk-forward backtest is a common backtest method which one continuously move a window, which is consist of a training period and a following testing period, along a long period of time. We performed a walk-forward backtest for total time period from 2008/01/01 to 2021/12/31, with a 260 week (about 5 year) training period and 52 week (about 1 year) testing period, and a moving interval of 52 week to analyze the performance of BL model. Note that our moving interval is identical to the testing period, thus this backtest can also be seen as a dynamic portfolio optimization without considering the transaction fee.

The goal for quantum classifiers is to classify $Y_{1}\in\{-1,1\}$ and $Y_{2}\in\{1,2\}$ in equation (6) after $t$ trading units (e.g. days/weeks) in the future for each asset in our pool, and predict the given feature vector $\mathbf{X}$ at the latest date of training period to construct $\eta s$ via equation (6), and thus the view matrix $\mathbf{Q}$. We can represent the feature of any given date $D$ with financial indicators of that day, and labeled $Y1$ and $Y2$ according to the mean (log) return $\overline{r}$ of the $t$ trading units since $D$. The labels are given by the following rules:

where $\mathbf{\sigma^{\prime}}=\frac{\sigma}{\sqrt{t}}$, and $\sigma$ is the standard deviation of $\{r\}$ of the corresponding assets. Since the testing period is 52 weeks, $t=52$ here.

In this work, we implemented two kinds of common quantum supervised learning models, which are quantum neural network (QNN)[12, 13], and quantum kernel method for support vector machine (QSVM)[11], and compared them with classical NN and SVM. Both quantum classifiers rely on a quantum circuit feature map $\mathcal{U}_{\Phi}(\overrightarrow{x})$ that encodes classical data into Hilbert space, but the training part works in different ways. For QNN, the feature map circuit is directly followed by an ansatz circuit with trainable parameters. Every data point that are encoded with the feature map and go through the unitary transformation of the ansatz, can be measured in computational basis and output a collection samples $\{b\}$. We interpret these measurement bit strings with a parity function that maps the samples to $sum(b)\ mod\ 2$, resulting in a binary probability distribution. We then use classical optimizer to minimize the squared loss function and update the parameters in the ansatz circuit. For QSVM, the quantum kernel can be represented by $\mathcal{U}_{\Phi}(\overrightarrow{x})^{\dagger}\mathcal{U}_{\Phi}(\overrightarrow{x})$, which the same feature maps as shown in Fig 2, followed by its own dagger. We measure the probability of $\ket{0}^{\otimes n}$, which has been shown in [11] to be equivalent to the inner product of the data points in feature space. Thus, we can apply this quantum kernel to the classical support vector machine (SVM).

We search among several kinds of feature map circuits, including the simple embedding 2(b) and the Pauli feature map 2(a), with different Pauli gates (e.g. [Z, YY], [Y, ZZ]…etc.) and entanglement structure (e.g. linear, circular, and full), finding that in our case, a simple embedding feature map without any entanglement structure (n=0) is enough to obtain high testing accurcies. QSVM with simple embedding feature map turns out to be slightly better than classical SVM with Gaussian kernel. The Pauli feature map, on the other hand, includes two-qubits rotational gates. It is more complex but did not work well on this task since we are training with only a small size and simple feature data. For predictions of more complex situations, which the feature dimension may be larger and more complex, one may be able to leverage the complexity of entangled quantum feature map. On the other hand, QNN and NN did not perform as well as SVM based ML model here.

Before actually solving the QUBO objective function, i.e. equation (9), we must first determine a proper penalty $\lambda$. $\lambda$ should be large enough to implement the constraint, but not unnecessarily large to dominate the Hamiltonian, making the risk-return part requires extreme precision to solve. Moreover, a large $\lambda$ also directly increases the coefficient of the Hamiltonian, and thus we need more measurements (circuit repetitions) to evaluate the expectation value to same precision. Here we choose $\lambda=1.0$ by empirically analyze the distribution of eigenvalue of the Hamiltonian. We believe that choosing $\lambda$ such that the eigenvalue of the lowest eigenstate that doesn’t satisfy the constraint is at least larger than the $50\%$ but less than $100\%$ of the eigenvalue of eigenstates that satisfy the constraint, as shown in Fig 4, will be preferable for VQEs.

We solved equation (9) with gate-based quantum computer simulator using a heuristic circuit VQE, and QAOA[14]. QAOA is also a type of VQE that has a specific variational form that depends on the problem Hamiltonian, and is related to the quantum adiabatic algorithm. QUBO problems are natively solvable for quantum computers since the objective function can easily be mapped to an Ising Hamiltonian by $x=\frac{1-Z}{2}$, where $x\in\{0,1\}$ and $Z\in\{1,-1\}$ is the Pauli Z operator. The mapped Ising Hamiltonian of equation (9) is

where $h_{ij}$ and $h_{i}$ are coefficients of corresponding terms and $C$ is the offset constant. This Ising Hamiltonian should in general be fully connected, unless the covariance between some of the assets are zero. With this problem Hamiltonian, we construct the heuristic ansatz for VQE shown in Fig 5, and QAOA ansatz is built following the instruction stated in [14]. From there, we then update the parameters in the circuit with Sequential Least Squares Programming (SLSQP) optimizer in order to minimize the Hamiltonian expectation value. After convergence, the optimal parameters is then fixed in the ansatz to perform final samplings, and the solution is chosen as the outcome bit string with lowest eigenvalue among the final measurements. We take only 5 measurement shots as final samplings here. There are two reasons to use such few shots, first is that once we optimized the parameter in the ansatz successfully, the probability distribution should concentrate on low energy eigenstates, thus we should expect needing only few final samplings to get a good solution with high probability. Secondly, this is also to mimic the situation of solving a large size problem, which the final shots we can get is always a tiny portion of the total search space, and should always not exceed the search space (otherwise one should rather perform a exhaustive search). Since the search space grows exponentially, a large amount of final samplings (compare to search space) is neither justified nor possible. We further analyzed the probability of getting a ”good” solution with random sampling in appendix A. We show that the probability of reaching such AR within 5 final samplings is very low, and provide a upper bound for maximum number of final sampling.

The heuristic VQE form we used requires $3N(p+1)$ parameters for an N assets optimization problem, where p is the repetition of the circuit. For QAOA, the entanglement has to be fully connected since the problem itself is. It might be more hardware demanding for most cases, but requires only $2p$ parameters, and the performance is guaranteed to improve (or at least the same) as $p$ increases[14]. We choose $p=4$ for heuristic VQE and $p=8$ for QAOA, and take 10 and 500 initial parameter guesses for each respectively. The initial guess that is eventually optimized to the lowest expectation value will served as the optimal parameter and used to perform final sampling.

We applied quantum algorithms to BL model and implement it for real world stock data, with walk-forward back test framework, i.e. a 262-week (5 years) in-sample period and a 52-week (1 year) out-sample period, starting from $2008-01-01$ and ends at $2021-12-31$, creating a total 9 continuous time segments. The training and testing stock price data are all downloaded from Yahoo finance weekly data, and the market capitalization data is obtained from [15]. The financial indicators for predicting investor’s views (TABLE I) are downloaded from ”Federal Reserve Bank of St. Louis” [16]. The asset pool contains 12 individual stocks randomly chosen from $S\&P500$ companies of different industries.

It requires several preprocessing before the financial indicators can become a feature vector for prediction tasks.
1. Backward moving average with a 3-week rolling window.
2. Standard normalization.
3. PCA reduce to 4 dimension.
4. Re-normalize each feature to $(0,2\pi]$.
We then label the these features of each week according to the rule in equation (14). In order to improve the ”trainability” of our data, we remove both feature vectors if they are too close in feature space but having different labels. The main idea of doing this is that we don’t want those feature vectors that give no information, or more specifically, require extremely high precision to distinguish each other, to distract our models. The feature vectors can then be encoded into the quantum circuit via our quantum feature map to train predictive models. We further extend the procedure to a global portfolio optimization case in Appendix B, which the asset pool contains 16 ETFs that approximately represents the global economy.

We here showcase in TABLE II the average prediction accuracy for $Y_{1}$, and $Y_{2}$ over total 9 time segments of each asset. We found that in average, $QSVM\approx SVM>NN>QNN$ in terms of testing accuracy. QSVM is also much faster to train compare to QNN. Here, we chose QSVM models to predict the investor’s view for future optimizations. The corresponding $\eta s$ predicted by QSVM for each assets is shown as a colored matrix in Fig 6.

Since the problem size here is not large, both classical and quantum kernel performed well. However, we may look forward to potential advantage for quantum kernel in predictions for more complex situations, as in real world use cases.

As stated above, we applied (a) VQE circuit in Fig 5 with $p=4$, (b) QAOA with $p=8$, and (c) exhaustive search to solve all QUBO problems. Before we show the results, we have to first re-define the standard of measuring how good a solution is. In general, people usually refer to the approximation ratio(AR), i.e. $\frac{E}{E^{*}}$ for a minimization tasks, or the relative error, i.e. $\frac{|E-E^{*}|}{E^{*}}$, where $E$ is the eigenvalue obtained from the algorithm, and $E^{*}$ is the exact minimum eigenvalue. However, for a portfolio optimization problem, $E^{*}$ may often be around $0$, and thus results in unreasonably large ARs or relative errors. Moreover, after we include the $1^{T}x=B$ constraint into the penalty term (i.e. equation (9)), the energy level of the Hamiltonian splits into chuncks, due to the large amount of extra energy a solution will gain if it violates the constraint (to different degrees), see Fig 4 as an example. Thus, if we look at the whole energy range of the Hamiltonian with penalty, which may be orders larger than the one without penalty, and construct an AR based on how close it is to the ground state, the solutions that simply meets the constrain will turn out to have very high AR and thus becomes hard to distinct a good solution from a bad one if they both meet the constraint. We thus here define AR as

which is equivalent to shifting the worse solution that meets the constraint $E_{w}$ to $0$. Note that $E^{*}\leq E\leq E_{w}=0$, and the solutions that violates the constraint will have negative AR (but positive AR do not necessarily meets the constraint).

The optimization result is shown in Fig 7. First, we look at the AR of the optimal circuits (7(a)), which is the corresponding AR of $\expectationvalue{H}{\psi(\overrightarrow{\theta^{*}})}$, where $\psi(\overrightarrow{\theta^{*}})$ is the ansatz with optimal parameters $\theta^{*}$. We see that VQE heuristic ansatzs are all solved to at least $0.9$, and the mean of AR is $0.96$. When we perform the final samplings, we found the outcome contains only a single eigenstates, meaning that the optimal ansatz has probability distribution highly concentrated around one or a few eigenstates. On the other hand, the optimal QAOA ansatzs wasn’t solve to AR as good as heuristic ansatzs, and the final samplings has a distribution over several eigenstates. However, we are still able to get a relatively good solution from only 5 samplings in our case, but it is not guaranteed. We further analyzed the optimal ansatzs by sampling 500000 shots and calculate the variance of eigenvalues in 7(b). As expected, the variance of VQE heuristic ansatzs are almost 0 and QAOA ansatz is large. Furthermore, although the heuristic ansatz still suffers from local minimum and is not easy to solve it to the ground state, a few number of initial guess is sufficient to solve it to a good AR. For QAOA, the result of optimization varies tremendously with the initial point, and we may usually get bad solutions (AR of ansatz $<0$) if we don’t provide enough number of initial guesses. We are not able to determine which ansatz form is better by only looking at the AR and variance of the optimal ansatzs, since they are both able to obtain a solution with good AR within 5 final samplings. However, in terms of the amount of initial guesses, circuit depth and connectivity, and the concentration of probability, we came to a conclusion that VQE heuristic ansatz should be preferred.

For backtesting, we adopt the view predicted by QSVM, and the show both the exact solution and the optimal solution obtain by VQE and QAOA. Here, the portfolio performance is measured by certainty-equivalent return (CER):

where, $\mu$ is the sample mean of (log) return, and $\sigma$ is the standard deviation of the (log) return. The backtesting period covers 9 continuous 52-week time segments, and is compared to 3 kinds of common portfolios, MPT model, naive portfolio ($1/N$ portfolio), and the market-cap-weighted index. In our case, BL model out performs MPT model in terms of both pure return and mean CER in a long continuous becktesting period. We also show that for VQE/QAOA solutions with high enough ARs could performed similarly as the optimal solution, and even may have a chance to outperform exact solutions. Note that we are not claiming a slightly worse solution will outperform exact solution, but a solution with high AR is enough to perform as good as exact solution. For example in the extended case in Appendix B, the performance is proportional to ARs. This means that we may not necessarily have to solve this problem to the ground state, and we can save computational cost (in terms of Hamiltonian evaluation and classical optimization) by requiring less precision.