# ! pip install qiskit qiskit-finance==0.2.1

import warnings
warnings.filterwarnings("ignore")
import qiskit.tools.jupyter
%qiskit_version_table

from qiskit import IBMQ
IBMQ.save_account('MY_API_TOKEN', overwrite=True)

# ! pip install eikon

import eikon as ek
print("Eikon version: ", ek.__version__)

eikon_key = open("eikon.txt","r")
ek.set_app_key(str(eikon_key.read()))
eikon_key.close()

import numpy as np
import pandas as pd
import datetime
import matplotlib.pyplot as plt
import seaborn as sns 
import qiskit
from qiskit import Aer, QuantumCircuit
from qiskit_finance.data_providers._base_data_provider import BaseDataProvider
from qiskit.finance.applications.ising import portfolio
from qiskit.circuit.library import TwoLocal
from qiskit.aqua import QuantumInstance
from qiskit.optimization.applications.ising.common import sample_most_likely
from qiskit.aqua.algorithms import VQE, QAOA, NumPyMinimumEigensolver
from qiskit.aqua.components.optimizers import COBYLA

class EikonDataProvider(BaseDataProvider):
    """
    The abstract base class for Eikon data_provider.
    """
    def __init__(self, stocks_list, start_date, end_date):
        '''
        stocks -> List of interested assets
        start -> start date to fetch historical data
        end -> 
        '''
        super().__init__()
        self._stocks = stocks_list
        self._start = start_date
        self._end = end_date
        self._data = []
        self.stock_data = pd.DataFrame()
        
    def run(self):
        self._data = []
        stocks_notfound = []
        stock_data = ek.get_timeseries(self._stocks,
                  start_date=self._start, 
                  end_date=self._end, 
                  interval='daily',
                  corax='adjusted')
        for ticker in self._stocks:
            stock_value = stock_data['CLOSE']
            self.stock_data[ticker] = stock_data['CLOSE']
            if stock_value.dropna().empty:
                stocks_notfound.append(ticker)
            self._data.append(stock_value)

# List of stocks 
stock_list = ['FB.O']

# Start Date
start_date = datetime.datetime(2020,12,1)

# End Date
end_date = datetime.datetime(2021,12,1)

# Set number of equities to the number of stocks
num_assets = len(stock_list)

# Set the risk factor
risk_factor = 0.7

# Set budget
budget = 2

# Scaling of budget penalty term will be dependant on the number of assets
penalty = num_assets

data = EikonDataProvider(stocks_list = stock_list, start_date=start_date, end_date=end_date)
data.run()

# Top 5 rows of data 
df = data.stock_data
df.head()

df.describe()

# Closing Price History
fig, ax = plt.subplots(figsize=(15, 8))
ax.plot(df)
plt.title('Close Price History')
plt.xlabel('Date',fontsize =20)
plt.ylabel('Price in USD',fontsize = 20)
ax.legend(df.columns.values)
plt.show()

df.tail()

strike_price = 315   # agreed upon strike price
T = 40 / 253         # 40 days to maturity

S = 310.6              # initial spot price
# vol = 0.4            # volatility of 40%
vol = df['FB.O'].pct_change().std()
r = 0.05             # annual interest rate of 4%

# resulting parameters for log-normal distribution
mu = ((r - 0.5 * vol**2) * T + np.log(S))
sigma = vol * np.sqrt(T)
mean = np.exp(mu + sigma**2/2)
variance = (np.exp(sigma**2) - 1) * np.exp(2*mu + sigma**2)
stddev = np.sqrt(variance)

# lowest and highest value considered for the spot price; in between, an equidistant discretization is considered.
low  = np.maximum(0, mean - 3*stddev)
high = mean + 3*stddev

low

high

import numpy as np
from qiskit.algorithms import AmplitudeEstimation
from qiskit_finance.circuit.library import EuropeanCallPricingObjective  # F
from qiskit.circuit.library import LogNormalDistribution, NormalDistribution
from qiskit_finance.applications import EuropeanCallPricing
from qiskit_finance.applications.estimation import EuropeanCallDelta
from qiskit import Aer
from qiskit.algorithms import IterativeAmplitudeEstimation
sim = QuantumInstance(Aer.get_backend("aer_simulator"), shots=1000)

# number of qubits to represent the uncertainty
num_uncertainty_qubits = 5

distribution = LogNormalDistribution(num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high))

european_call_pricing = EuropeanCallPricing(num_state_qubits=num_uncertainty_qubits,
                                            strike_price=strike_price,
                                            rescaling_factor=0.25,  # approximation constant for payoff function
                                            bounds=(low, high),
                                            uncertainty_model=distribution)

problem = european_call_pricing.to_estimation_problem()

problem.state_preparation.draw('mpl', style='iqx')
plt.show()

epsilon = 0.01  # determines final accuracy
alpha = 0.05  # determines how certain we are of the result

ae = IterativeAmplitudeEstimation(epsilon, alpha=alpha, quantum_instance=sim)

result = ae.estimate(problem)

conf_int_result = np.array(result.confidence_interval_processed)
print("Esimated value: \t%.4f" % european_call_pricing.interpret(result))
print("Confidence interval: \t[%.4f, %.4f]" % tuple(conf_int_result))

x = distribution.values
y = distribution.probabilities
plt.figure(figsize=(15,6))
plt.bar(x, y, width=0.3)
plt.xticks(x, size=15, rotation=90)
plt.yticks(size=10)
plt.grid()
plt.xlabel('Spot Price at Maturity $S_T$ (\$)', size=15)
plt.ylabel('Probability ($\%$)', size=15)
plt.show()

# plot exact payoff function (evaluated on the grid of the uncertainty model)
x = distribution.values
y = np.maximum(0, x - strike_price)
plt.figure(figsize=(15,6))
plt.plot(x, y, "ro-")
plt.grid()
plt.title("Payoff Function for Call Option", size=15)
plt.xlabel("Spot Price", size=15)
plt.ylabel("Payoff", size=15)
plt.xticks(x, size=15, rotation=90)
plt.yticks(size=15)
plt.show()

# evaluate exact expected value (normalized to the [0, 1] interval)
exact_value = np.dot(distribution.probabilities, y)
exact_delta = sum(distribution.probabilities[x >= strike_price])
print("exact expected value:\t%.4f" % exact_value)
print("exact delta value:   \t%.4f" % exact_delta)

european_call_delta = EuropeanCallDelta(
    num_state_qubits=num_uncertainty_qubits,
    strike_price=strike_price,
    bounds=(low, high),
    uncertainty_model=distribution,
)

european_call_delta._objective.decompose().draw('mpl', style='iqx')
plt.show()

european_call_delta_circ = QuantumCircuit(european_call_delta._objective.num_qubits)
european_call_delta_circ.append(distribution, range(num_uncertainty_qubits))
european_call_delta_circ.append(
    european_call_delta._objective, range(european_call_delta._objective.num_qubits)
)

european_call_delta_circ.draw('mpl', style='iqx')
plt.show()

# set target precision and confidence level
epsilon = 0.01
alpha = 0.05

problem = european_call_delta.to_estimation_problem()

# construct amplitude estimation
ae_delta = IterativeAmplitudeEstimation(epsilon, alpha=alpha, quantum_instance=sim)

result_delta = ae_delta.estimate(problem)

conf_int = np.array(result_delta.confidence_interval_processed)
print("Exact delta:    \t%.4f" % exact_delta)
print("Esimated value: \t%.4f" % european_call_delta.interpret(result_delta))
print("Confidence interval: \t[%.4f, %.4f]" % tuple(conf_int))

put_distribution = LogNormalDistribution(num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high))

put_european_call_pricing = EuropeanCallPricing(num_state_qubits=num_uncertainty_qubits,
                                            strike_price=strike_price,
                                            rescaling_factor=0.25,  # approximation constant for payoff function
                                            bounds=(low, high),
                                            uncertainty_model=put_distribution)

put_problem = put_european_call_pricing.to_estimation_problem()

epsilon = 0.01  # determines final accuracy
alpha = 0.05  # determines how certain we are of the result

ae = IterativeAmplitudeEstimation(epsilon, alpha=alpha, quantum_instance=sim)

put_result = ae.estimate(put_problem)

put_conf_int_result = np.array(put_result.confidence_interval_processed)
print("Esimated value: \t%.4f" % european_call_pricing.interpret(put_result))
print("Confidence interval: \t[%.4f, %.4f]" % tuple(put_conf_int_result))

x = put_distribution.values
y = put_distribution.probabilities
plt.figure(figsize=(15,6))
plt.bar(x, y, width=0.3)
plt.xticks(x, size=12, rotation=90)
plt.yticks(size=12)
plt.xlabel('Spot Price at Maturity $S_T$ (\$)', size=15)
plt.ylabel('Probability ($\%$)', size=15)
plt.show()

# plot exact payoff function (evaluated on the grid of the uncertainty model)
x = distribution.values
y = np.maximum(0, strike_price-x)
plt.figure(figsize=(15,6))
plt.plot(x, y, "ro-")
plt.title("Payoff Function for Put Option", size=15)
plt.xlabel("Spot Price", size=15)
plt.ylabel("Payoff", size=15)
plt.xticks(x, size=12, rotation=90)
plt.yticks(size=12)
plt.show()

# evaluate exact expected value (normalized to the [0, 1] interval)
put_exact_value = np.dot(put_distribution.probabilities, y)
conf_int = np.array(put_result.confidence_interval_processed)
print("Exact Expected Value:\t%.4f" % exact_value)
print("Esimated value: \t%.4f" % european_call_delta.interpret(put_result))
print("Confidence interval: \t[%.4f, %.4f]" % tuple(conf_int))