# データ取得に必要なpandas, pandas_datareaderのインストール
!pip install pandas pandas_datareader

import numpy as np
import pandas as pd
import pandas_datareader.data as web
import datetime
import matplotlib.pyplot as plt

from qiskit import Aer
from qiskit.circuit.library import QFT
from qiskit.aqua import QuantumInstance, aqua_globals
from qiskit.quantum_info import state_fidelity
from qiskit.aqua.algorithms import HHL, NumPyLSsolver
from qiskit.aqua.components.eigs import EigsQPE
from qiskit.aqua.components.reciprocals import LookupRotation
from qiskit.aqua.operators import MatrixOperator
from qiskit.aqua.components.initial_states import Custom

# 銘柄選択
codes = ["GOOG", 'AAPL', 'FB', 'AMZN'] # GAFA

# 2017年の1年間のデータを使用
start = datetime.datetime(2017, 1, 1)
end = datetime.datetime(2017, 12, 31)

# Yahoo! Financeから日次の株価データを取得
data = web.DataReader(codes, 'yahoo', start, end)

df = data['Adj Close'] 

## 直近のデータの表示
display(df.tail())

## 株価をプロットしてみる
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])

daily_return = df.pct_change()
display(daily_return.tail())

expected_return = daily_return.dropna(how='all').mean() * 252 # 年率換算のため年間の営業日数252を掛ける
print(expected_return)

cov = daily_return.dropna(how='all').cov() * 252 # 年率換算のため
display(cov)

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)

mu = 0.1 # ポートフォリオのリターン
xi = 1.0 
mu_xi_0 = np.append(np.array([mu, xi]), np.zeros_like(R))
print(mu_xi_0)

matrix, vector, truncate_powerdim, truncate_hermitian = HHL.matrix_resize(W, mu_xi_0)

# 2**3=8次元の行列になっているか確認
print(matrix)

def create_eigs(matrix, num_ancillae, num_time_slices, negative_evals):
    ne_qfts = [None, None]
    if negative_evals:
        num_ancillae += 1
        # The QFT and IQFT circuits for handling negative eigenvalues
        ne_qfts = [QFT(num_ancillae - 1), QFT(num_ancillae - 1).inverse()]

    """
    Specifically, this class is based on PhaseEstimationCircuit with no measurements 
    and has additional handling of negative eigenvalues
    """
    return EigsQPE(MatrixOperator(matrix=matrix), # The Hamiltonian Operator object
                   QFT(num_ancillae).inverse(),
                   num_time_slices=num_time_slices,
                   num_ancillae=num_ancillae,
                   expansion_mode='trotter', 
                   expansion_order=2, # The suzuki expansion order, has a minimum value of 1
                   evo_time=None,  # This is t, can set to: np.pi*3/4
                   negative_evals=negative_evals, # Set True to indicate negative eigenvalues need to be handled
                   ne_qfts=ne_qfts)

# 扱う問題の次元
orig_size = len(vector)

# 固有値探索の初期化
eigs = create_eigs(matrix, 7, 50, True)
num_q, num_a = eigs.get_register_sizes()

# 制御回転の初期化
c = 0.5*(2 * np.pi * (1. / 2**(num_a) ))
reci = LookupRotation(scale=c, negative_evals=eigs._negative_evals, evo_time=eigs._evo_time)

# 初期状態の設定
init_state = Custom(num_q, state_vector=vector)

algo = HHL(
    matrix,
    vector, 
    truncate_powerdim, 
    truncate_hermitian, 
    eigs, # The eigenvalue estimation instance
    init_state, # The initial quantum state preparation
    reci, # The eigenvalue reciprocal and controlled rotation instance
    num_q, # Number of qubits required for the matrix Operator instance
    num_a, # Number of ancillary qubits for Eigenvalues instance
    orig_size # The original dimension of the problem (if truncate_powerdim)
    )

## 結果
x_HHL = algo.run(QuantumInstance(Aer.get_backend('statevector_simulator')))
## 厳密解
# x_exact = np.linalg.lstsq(matrix, vector, rcond=0)[0]
x_exact = NumPyLSsolver(matrix, vector).run()

# 解の信頼度を求める
def fidelity(hhl, ref):
    solution_hhl_normed = hhl / np.linalg.norm(hhl)
    solution_ref_normed = ref / np.linalg.norm(ref)
    fidelity = state_fidelity(solution_hhl_normed, solution_ref_normed)
    print("Fidelity:\t\t %f" % fidelity)

fidelity(x_HHL.solution, x_exact.solution)

# the number of qubits required
print("circuit_width:\t", x_HHL['circuit_info']['width'])
# the maximum number of gates
print("circuit_depth:\t", x_HHL['circuit_info']['depth'])
# To verify whether running the circuit on current real hardware is feasible
print("CNOT gates:\t", x_HHL['circuit_info']['operations']['cx'])

w_opt_HHL = x_HHL.solution[2:6] 
w_opt_exact = x_exact.solution[2:6] 
w_opt = pd.DataFrame(np.vstack([w_opt_exact, w_opt_HHL]).T, index=df.columns, columns=['exact', 'HHL'])
w_opt

w_opt.plot.bar()

# 2018年の1年間のデータを使用
start = datetime.datetime(2017, 12, 30)
end = datetime.datetime(2018, 12, 31)

# Yahoo! Financeから日次の株価データを取得
data = web.DataReader(codes, 'yahoo', start, end)

df2018 = data['Adj Close'] 

display(df2018.tail())

## 株価をプロットしてみる
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])

# ポートフォリオの資産価値の推移
pf_value = df2018.dot(w_opt)
pf_value.head()

# exact と HHLで初期金額が異なることがありうるので、期初の値で規格化したリターンをみる。
pf_value.exact = pf_value.exact / pf_value.exact[0] 
pf_value.HHL = pf_value.HHL / pf_value.HHL[0] 
print(pf_value.tail())

pf_value.plot(figsize=(9, 6))

pf_value.pct_change().std() * np.sqrt(252) ## 年率換算