7-3. Portfolio optimization using HHL algorithm
In this section, let us calculate the optimal portfolio (asset allocation) based on the data of past stock price fluctuations, referring to the paper [1]. Portfolio optimization is one of the problems that is expected to be solved faster than before by using the HHL algorithm learned in section 7-1. Specifically, this time, we consider the problem of what kind of asset allocation will give the highest return with the lowest risk when investing in stocks of the four GAFA (Google, Apple, Facebook, and Amazon) companies.
Stock Price Data Acquisition
First, we obtain the stock price data of each company.
Use daily data of GAFA 4 companies
To obtain stock price data, use pandas_datareader to obtain from Yahoo! Finance database.
Adjusted closing price in dollars (Adj. Close) is used for the stock price.
[1]:
# Install pandas, pandas_datareader required for data acquisition
# !pip install pandas pandas_datareader
[2]:
import numpy as np
import pandas as pd
import pandas_datareader.data as web
import datetime
import matplotlib.pyplot as plt
[ ]:
# Issue Selection
codes = ['GOOG', 'AAPL', 'FB', 'AMZN'] # GAFA
# Using one year of data from 2017
start = datetime.datetime(2017, 1, 1)
end = datetime.datetime(2017, 12, 31)
# Obtain daily stock price data from Yahoo! Finance
data = web.DataReader(codes, 'yahoo', start, end)
df = data['Adj Close']
## Display of most recent data
display(df.tail())
[4]:
## Plot the stock price
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))
df.loc[:,['AAPL', 'FB']].plot(ax=axes[0])
df.loc[:,['GOOG', 'AMZN']].plot(ax=axes[1])
[4]:
<matplotlib.axes._subplots.AxesSubplot at 0x114957850>
*The reason why the four stocks are divided into two groups here is because the stock price values are close to each other and easy to see when plotting, and there is no deeper meaning.
Data Preprocessing
Next, the acquired stock prices are converted to daily returns to obtain some statistics.
Conversion to Daily Returns
The daily return (rate of change) \(y_t\) (where \(t\) is the date) of an individual stock is defined as follows
This is obtained with the pct_change() method of the pandas DataFrame.
[5]:
daily_return = df.pct_change()
display(daily_return.tail())
| Symbols | GOOG | AAPL | FB | AMZN |
|---|---|---|---|---|
| Date | ||||
| 2017-12-22 | -0.003300 | 0.000000 | -0.001409 | -0.005448 |
| 2017-12-26 | -0.003188 | -0.025370 | -0.006828 | 0.007190 |
| 2017-12-27 | -0.006974 | 0.000176 | 0.009262 | 0.004674 |
| 2017-12-28 | -0.001172 | 0.002814 | 0.001689 | 0.003248 |
| 2017-12-29 | -0.001660 | -0.010814 | -0.008206 | -0.014021 |
Expected Returns
Find the expected return \(\vec R\) for each stock. Here we use the arithmetic mean of the historical returns:
[6]:
expected_return = daily_return.dropna(how='all').mean() * 252 # Multiply by 252 annual operating days for annualization
print(expected_return)
Symbols
GOOG 0.300215
AAPL 0.411192
FB 0.430156
AMZN 0.464567
dtype: float64
Variance/Covariance Matrix
The sample unbiased variance/covariance matrix \(\Sigma\) of the returns is defined by
[7]:
cov = daily_return.dropna(how='all').cov() * 252 # For annualized basis
display(cov)
| Symbols | GOOG | AAPL | FB | AMZN |
|---|---|---|---|---|
| Symbols | ||||
| GOOG | 0.023690 | 0.013303 | 0.018382 | 0.021614 |
| AAPL | 0.013303 | 0.031119 | 0.016291 | 0.018877 |
| FB | 0.018382 | 0.016291 | 0.028855 | 0.023337 |
| AMZN | 0.021614 | 0.018877 | 0.023337 | 0.044120 |
Portfolio Optimization
Now that we are ready, let’s tackle portfolio optimization.
First, the portfolio (i.e., asset allocation) is represented by a 4-component vector \(\vec{w} = (w_0,w_1,w_2,w_3)^T\). For example, if \(\vec{w}=(1,0,0,0)\), it means a portfolio with 100% of its assets invested in Google stock.
Let us consider a portfolio that satisfies the following equation.
This formula means
“minimizing the variance of the portfolio’s return”
under the following condition
“The expected return (average of returns) of the portfolio is \(\mu\)”
“the sum of the weights invested in the portfolio is 1” (\(\vec 1 = (1,1,1,1)^T\))
In other words, the best portfolio is the one that minimizes the variance (risk) of the desired future return of only \(\mu\) as much as possible. Such a problem setup is called Markowitz’s mean-variance approach and is one of the foundational ideas of modern financial engineering.
Using Lagrange’s undetermined multiplier method, \(\vec{w}\) satisfying the above conditions can be obtained by solving the linear equation
where \(\eta, \theta\) are the parameters of Lagrange’s undetermined multiplier method. Therefore, to find the optimal portfolio \(\vec w\), we can solve the simultaneous equation (1) for \(\vec w\). We have now attributed the portfolio optimization problem to a linear equation for which the HHL algorithm can be used.
Creating the matrix W
[8]:
R = expected_return.values
Pi = np.ones(4)
S = cov.values
row1 = np.append(np.zeros(2), R).reshape(1,-1)
row2 = np.append(np.zeros(2), Pi).reshape(1,-1)
row3 = np.concatenate([R.reshape(-1,1), Pi.reshape(-1,1), S], axis=1)
W = np.concatenate([row1, row2, row3])
np.set_printoptions(linewidth=200)
print(W)
[[0. 0. 0.30021458 0.41119151 0.43015563 0.46456748]
[0. 0. 1. 1. 1. 1. ]
[0.30021458 1. 0.02369003 0.01330333 0.01838175 0.0216144 ]
[0.41119151 1. 0.01330333 0.03111917 0.01629131 0.01887668]
[0.43015563 1. 0.01838175 0.01629131 0.02885482 0.02333747]
[0.46456748 1. 0.0216144 0.01887668 0.02333747 0.04412049]]
[9]:
## Check eigenvalue of W -> fits in [-pi, pi]
print(np.linalg.eigh(W)[0])
[-2.11207187 -0.10947986 0.01121933 0.01864265 0.11919724 2.20027702]
Creating the right-hand side vector
Given the expected return \(\mu\) of the portfolio below, we can calculate the least risky portfolio that will yield such a return. The \(\mu\) can be set freely. In general, the higher the expected return, the greater the risk, but we will use 10% as an example (this is a very bearish one, since it is a time when GAFA stocks are going gangbusters).
[10]:
mu = 0.1 # Portfolio returns (parameters to be put in hand)
xi = 1.0
mu_xi_0 = np.append(np.array([mu, xi]), np.zeros_like(R)) ## Vector on the right side of equation (1)
print(mu_xi_0)
[0.1 1. 0. 0. 0. 0. ]
Extend the matrix so that it can be handled by a quantum system
Since \(W\) is 6-dimensional, it is computable in a quantum system with 3 qubits (\(2^3 = 8\)). Therefore, we also make a matrix and a vector with the extended two dimensions filled with zeros.
[11]:
nbit = 3 ## Number of bits used for state
N = 2**nbit
W_enl = np.zeros((N, N)) ## enl stands for enlarged
W_enl[:W.shape[0], :W.shape[1]] = W.copy()
mu_xi_0_enl = np.zeros(N)
mu_xi_0_enl[:len(mu_xi_0)] = mu_xi_0.copy()
We are now ready to solve the simultaneous equation (1).
Calculating the minimum variance portfolio using the HHL algorithm
Now, let us solve the simultaneous linear equations (1) using the HHL algorithm. First, as a preliminary step, we prepare the following functions.
A function
input_state_gatethat returns a quantum circuit that converts the quantum state to \(|0\cdots0\rangle \to \sum_i x_i |i \rangle\) according to the classical data \(\mathbf{x}\) This function should be built using the qRAM concept, but since we are using a simulator, we will implement it as a non-unitary gate this time. Also, standardization is ignored.)Function
CPhaseGatethat returns a control phase gate.function
QFT_gatethat returns a gate to perform quantum Fourier transform.
[12]:
# Install Qulacs
# !pip install qulacs
## Run only if you are in a Google Colaboratory / (Linux or Mac) jupyter notebook environment.
## Qulacs errors will be output normally.
!pip3 install wurlitzer
%load_ext wurlitzer
[13]:
import numpy as np
from qulacs import QuantumCircuit, QuantumState, gate
from qulacs.gate import merge, Identity, H, SWAP
def input_state_gate(start_bit, end_bit, vec):
"""
Making a quantum gate which transform |0> to \sum_i x[i]|i>m where x[i] is input vector.
!!! this uses 2**n times 2**n matrix, so it is quite memory-cosuming.
!!! this gate is not unitary (we assume that the input state is |0>)
Args:
int start_bit: first index of qubit which the gate applies
int end_bit: last index of qubit which the gate applies
np.ndarray vec: input vector.
Returns:
qulacs.QuantumGate
"""
nbit = end_bit - start_bit + 1
assert vec.size == 2**nbit
mat_0tox = np.eye(vec.size, dtype=complex)
mat_0tox[:,0] = vec
return gate.DenseMatrix(np.arange(start_bit, end_bit+1), mat_0tox)
def CPhaseGate(target, control, angle):
"""
Create controlled phase gate diag(1,e^{i*angle}) with controll. (Qulacs.gate is requried)
Args:
int target: index of target qubit.
int control: index of control qubit.
float64 angle: angle of phase gate.
Returns:
QuantumGateBase.DenseMatrix: diag(1, exp(i*angle)).
"""
CPhaseGate = gate.DenseMatrix(target, np.array( [[1,0], [0,np.cos(angle)+1.j*np.sin(angle)]]) )
CPhaseGate.add_control_qubit(control, 1)
return CPhaseGate
def QFT_gate(start_bit, end_bit, Inverse = False):
"""
Making a gate which performs quantum Fourier transfromation between start_bit to end_bit.
(Definition below is the case when start_bit = 0 and end_bit=n-1)
We associate an integer j = j_{n-1}...j_0 to quantum state |j_{n-1}...j_0>.
We define QFT as
|k> = |k_{n-1}...k_0> = 1/sqrt(2^n) sum_{j=0}^{2^n-1} exp(2pi*i*(k/2^n)*j) |j>.
then, |k_m > = 1/sqrt(2)*(|0> + exp(i*2pi*0.j_{n-1-m}...j_0)|1> )
When Inverse=True, the gate represents Inverse QFT,
|k> = |k_{n-1}...k_0> = 1/sqrt(2^n) sum_{j=0}^{2^n-1} exp(-2pi*i*(k/2^n)*j) |j>.
Args:
int start_bit: first index of qubits where we apply QFT.
int end_bit: last index of qubits where we apply QFT.
bool Inverse: When True, the gate perform inverse-QFT ( = QFT^{\dagger}).
Returns:
qulacs.QuantumGate: QFT gate which acts on a region between start_bit and end_bit.
"""
gate = Identity(start_bit) ## make empty gate
n = end_bit - start_bit + 1 ## size of QFT
## loop from j_{n-1}
for target in range(end_bit, start_bit-1, -1):
gate = merge(gate, H(target)) ## 1/sqrt(2)(|0> + exp(i*2pi*0.j_{target})|1>)
for control in range(start_bit, target):
gate = merge( gate, CPhaseGate(target, control, (-1)**Inverse * 2.*np.pi/2**(target-control+1)) )
## perform SWAP between (start_bit + s)-th bit and (end_bit - s)-th bit
for s in range(n//2): ## s runs 0 to n//2-1
gate = merge(gate, SWAP(start_bit + s, end_bit - s))
## return final circuit
return gate
First, set the necessary parameters for the HHL algorithm. Set the clock register quantum bit number reg_nbit to 7 and the coefficient scale_fac used for scaling the matrix \(W\) to 1 (i.e., do not scale it). Also, the coefficient \(c\) used for the control rotation gate is kept at half of the smallest number of nonzeros that can be represented by the reg_nbit bits.
[14]:
# number of registers used for phase estimation
reg_nbit = 7
## Factor to scale W_enl
scale_fac = 1.
W_enl_scaled = scale_fac * W_enl
## Minimum value to be assumed as an eigenvalue of W_enl_scaled
## In this case, since the projection succeeds 100%, we set the value as a constant multiple of the minimum value that can be represented by the register.
C = 0.5*(2 * np.pi * (1. / 2**(reg_nbit) ))
The core of the HHL algorithm will be written. In this article, we use the simulator qulacs, so various simplifications have been made. This is an implementation to get a sense of how the HHL algorithm works.
The part of preparing the input state \(|\mathbf{b}\rangle\) is simplified.
The part of \(e^{iA}\) used in the quantum phase estimation algorithm is the diagonalization of \(A\) by a classical computer.
Control rotation gate that takes inverses is also implemented by a classical computer..
Measure the projection to the auxiliary bit \(|0 \rangle{}_{S}\) and treat only the states where the measurement result
0is obtained. (For convenience of implementation, the definition of the action of the control rotation gate is reversed from section 7-1)
[15]:
from functools import reduce
## diagonalization. AP = PD <-> A = P*D*P^dag
D, P = np.linalg.eigh(W_enl_scaled)
#####################################
### Create an HHL quantum circuit. Starting from the 0th bit,
### we have the bits in the space where A acts (0th ~ nbit-1th),
### the register bits (nbit th ~ nbit+reg_nbit-1 th), and
### the bits for conditional rotation (nbit+reg_nbit th).
#####################################
total_qubits = nbit + reg_nbit + 1
total_circuit = QuantumCircuit(total_qubits)
## ------ Prepare vector b input to the 0th~(nbit-1)th bit ------
## Normally we should use qRAM algorithm, but here we use our own input gate.
## In qulacs, you can also implement it with state.load(b_enl).
state = QuantumState(total_qubits)
state.set_zero_state()
b_gate = input_state_gate(0, nbit-1, mu_xi_0_enl)
total_circuit.add_gate(b_gate)
## ------- Hadamard gate on register bit -------
for register in range(nbit, nbit+reg_nbit): ## from nbit to nbit+reg_nbit-1
total_circuit.add_H_gate(register)
## ------- Implement phase estimation -------
## U := e^{i*A*t) and its eigenvalues are diag( {e^{i*2pi*phi_k}}_{k=0, ... N-1)).
## Implement \sum_j |j><j| exp(i*A*t*j) to register bits
for register in range(nbit, nbit+reg_nbit):
## Implement U^{2^{register-nbit}}.
## Use diagonalized results.
U_mat = reduce(np.dot, [P, np.diag(np.exp( 1.j * D * (2**(register-nbit)) )), P.T.conj()] )
U_gate = gate.DenseMatrix(np.arange(nbit), U_mat)
U_gate.add_control_qubit(register, 1) ## add control bit
total_circuit.add_gate(U_gate)
## ------- Perform inverse QFT to register bits -------
total_circuit.add_gate(QFT_gate(nbit, nbit+reg_nbit-1, Inverse=True))
## ------- multiply conditional rotation -------
## The eigenvalue of A*t corresponding to the register |phi> is l = 2pi * 0. phi = 2pi * (phi / 2**reg_nbit).
## The definition of conditional rotation is (opposite of the text)
## |phi>|0> -> C/(lambda)|phi>|0> + sqrt(1 - C^2/(lambda)^2)|phi>|1>.
## Since this is a classical simulation, the gate is made explicitly.
condrot_mat = np.zeros( (2**(reg_nbit+1), (2**(reg_nbit+1))), dtype=complex)
for index in range(2**reg_nbit):
lam = 2 * np.pi * (float(index) / 2**(reg_nbit) )
index_0 = index ## integer which represents |index>|0>
index_1 = index + 2**reg_nbit ## integer which represents |index>|1>
if lam >= C:
if lam >= np.pi: ## Since we have scaled the eigenvalues in [-pi, pi] beforehand, [pi, 2pi] corresponds to a negative eigenvalue
lam = lam - 2*np.pi
condrot_mat[index_0, index_0] = C / lam
condrot_mat[index_1, index_0] = np.sqrt( 1 - C**2/lam**2 )
condrot_mat[index_0, index_1] = - np.sqrt( 1 - C**2/lam**2 )
condrot_mat[index_1, index_1] = C / lam
else:
condrot_mat[index_0, index_0] = 1.
condrot_mat[index_1, index_1] = 1.
## Convert to DenseGate and implement
condrot_gate = gate.DenseMatrix(np.arange(nbit, nbit+reg_nbit+1), condrot_mat)
total_circuit.add_gate(condrot_gate)
## ------- Perform QFT to register bits -------
total_circuit.add_gate(QFT_gate(nbit, nbit+reg_nbit-1, Inverse=False))
## ------- Implement the inverse of phase estimation (U^\dagger = e^{-iAt}) -------
for register in range(nbit, nbit+reg_nbit): ## from nbit to nbit+reg_nbit-1
## implement {U^{\dagger}}^{2^{register-nbit}}
## use diagonalized results.
U_mat = reduce(np.dot, [P, np.diag(np.exp( -1.j* D * (2**(register-nbit)) )), P.T.conj()] )
U_gate = gate.DenseMatrix(np.arange(nbit), U_mat)
U_gate.add_control_qubit(register, 1) ## add a control bit
total_circuit.add_gate(U_gate)
## ------- act Hadamard gate on register bit -------
for register in range(nbit, nbit+reg_nbit):
total_circuit.add_H_gate(register)
## ------- Project auxiliary bits to 0. Implemented in qulacs as a non-unitary gate -------
total_circuit.add_P0_gate(nbit+reg_nbit)
#####################################
### Run the HHL quantum circuit and retrieve the result
#####################################
total_circuit.update_quantum_state(state)
## The 0th to (nbit-1)th bit corresponds to the calculation result |x>.
result = state.get_vector()[:2**nbit].real
x_HHL = result/C * scale_fac
Comparing the solution x_HHL by the HHL algorithm with the solution x_exact by diagonalization of the usual classical computation, we can see that they are roughly in agreement. (There are several parameters that determine the accuracy of the HHL algorithm (e.g. reg_nbit), so please try changing them and experiment with them).
[16]:
## Exact solution
x_exact = np.linalg.lstsq(W_enl, mu_xi_0_enl, rcond=0)[0]
print("HHL: ", x_HHL)
print("exact:", x_exact)
rel_error = np.linalg.norm(x_HHL- x_exact) / np.linalg.norm(x_exact)
print("rel_error", rel_error)
HHL: [ 0.09580738 -0.04980738 2.36660125 0.09900883 -0.47774813 -0.98438791 0. 0. ]
exact: [ 0.15426894 -0.07338059 2.29996915 0.17711988 -0.66526695 -0.81182208 0. 0. ]
rel_error 0.11097291393510306
If we take out only the actual weighting part
[17]:
w_opt_HHL = x_HHL[2:6]
w_opt_exact = x_exact[2:6]
w_opt = pd.DataFrame(np.vstack([w_opt_exact, w_opt_HHL]).T, index=df.columns, columns=['exact', 'HHL'])
w_opt
[17]:
| exact | HHL | |
|---|---|---|
| Symbols | ||
| GOOG | 2.299969 | 2.366601 |
| AAPL | 0.177120 | 0.099009 |
| FB | -0.665267 | -0.477748 |
| AMZN | -0.811822 | -0.984388 |
[18]:
w_opt.plot.bar()
[18]:
<matplotlib.axes._subplots.AxesSubplot at 0x1147add90>
*Stocks with negative weights represent “short selling” (borrowing shares and selling them. It is a technique whereby profits can be made during periods of falling stock prices). In this case, the target return was set at 10%, which is quite small for GAFA stocks (expected return of 30-40% alone), so short selling is likely to have lowered the overall expected return.
Appendix: Back Testing
The verification of an investment rule derived from historical data using subsequent data is called “back-testing,” and is important for measuring the effectiveness of the investment rule. In this section, we will observe how much asset values change in the following year, 2018, when invested in a portfolio constructed from 2017 data as described above.
[19]:
# Using one year data of 2018
start = datetime.datetime(2017, 12, 30)
end = datetime.datetime(2018, 12, 31)
# Obtain daily stock price data from Yahoo! Finance
data = web.DataReader(codes, 'yahoo', start, end)
df2018 = data['Adj Close']
display(df2018.tail())
| Symbols | GOOG | AAPL | FB | AMZN |
|---|---|---|---|---|
| Date | ||||
| 2018-12-24 | 976.219971 | 144.656540 | 124.059998 | 1343.959961 |
| 2018-12-26 | 1039.459961 | 154.843475 | 134.179993 | 1470.900024 |
| 2018-12-27 | 1043.880005 | 153.838562 | 134.520004 | 1461.640015 |
| 2018-12-28 | 1037.079956 | 153.917389 | 133.199997 | 1478.020020 |
| 2018-12-31 | 1035.609985 | 155.405045 | 131.089996 | 1501.969971 |
[20]:
## Plot the stock price
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))
df2018.loc[:,['AAPL', 'FB']].plot(ax=axes[0])
df2018.loc[:,['GOOG', 'AMZN']].plot(ax=axes[1])
[20]:
<matplotlib.axes._subplots.AxesSubplot at 0x1173c5cd0>
[21]:
# Changes in portfolio asset values
pf_value = df2018.dot(w_opt)
pf_value.head()
[21]:
| exact | HHL | |
|---|---|---|
| Date | ||
| 2017-12-29 | 1368.986737 | 1257.159152 |
| 2018-01-02 | 1393.124262 | 1279.864457 |
| 2018-01-03 | 1418.828873 | 1304.724242 |
| 2018-01-04 | 1423.832841 | 1308.934869 |
| 2018-01-05 | 1443.032292 | 1326.138623 |
[22]:
# Since the initial amount may differ between exact and HHL, we look at returns normalized by the value at the beginning of the period.
pf_value.exact = pf_value.exact / pf_value.exact[0]
pf_value.HHL = pf_value.HHL / pf_value.HHL[0]
print(pf_value.tail())
exact HHL
Date
2018-12-24 0.801548 0.749625
2018-12-26 0.828918 0.766234
2018-12-27 0.841539 0.781597
2018-12-28 0.821053 0.756478
2018-12-31 0.805599 0.735876
[23]:
pf_value.plot(figsize=(9, 6))
[23]:
<matplotlib.axes._subplots.AxesSubplot at 0x11475b8d0>
The stocks of all GAFA companies except Amazon were soft in 2018, resulting in a loss of roughly -20%. The EXACT solution seems to be somewhat better. Incidentally, since what we originally did was risk minimization, we also calculated the risk for the past year, and the EXACT solution resulted in smaller risk.
[24]:
pf_value.pct_change().std() * np.sqrt(252) ## annualized
[24]:
exact 0.402005
HHL 0.501925
dtype: float64
Reference
[1] P. Rebentrost and S. Lloyd, “Quantum computational finance: quantum algorithm for portfolio optimization“, https://arxiv.org/abs/1811.03975