Global Data Quantumのポートフォリオ・オプティマイザーによる動的なポートフォリオ最適化

!pip install qiskit-ibm-catalog
!pip install pandas
!pip install matplotlib
!pip install yfinance

from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(
    channel="ibm_quantum_platform",
    instance="INSTANCE_CRN",
    # For `token`, use the 44-character API_KEY you created
    # and saved from the IBM Quantum Platform Home dashboard
    token="YOUR_API_KEY",
)

# Access function
dpo_solver = catalog.load("global-data-quantum/quantum-portfolio-optimizer")

import yfinance as yf
import pandas as pd

# List of IBEX 35 symbols
symbols = [
    "ACS.MC",
    "ITX.MC",
    "FER.MC",
    "ELE.MC",
    "SCYR.MC",
    "AENA.MC",
    "AMS.MC",
]

start_date = "2022-11-01"
end_date = "2023-4-01"

series_list = []
symbol_names = [symbol.replace(".", "_") for symbol in symbols]

# Create a full date index including weekends
full_index = pd.date_range(start=start_date, end=end_date, freq="D")

for symbol, name in zip(symbols, symbol_names):
    print(f"Downloading data for {symbol}...")
    data = yf.download(symbol, start=start_date, end=end_date)["Close"]
    data.name = name

    # Reindex to include weekends
    data = data.reindex(full_index)

    # Fill missing values (for example, weekends or holidays)
    # by forward/backward fill
    data.ffill(inplace=True)
    data.bfill(inplace=True)

    series_list.append(data)

# Combine all series into a single DataFrame
assets = pd.concat(series_list, axis=1)

# Convert index to string for consistency
assets.index = assets.index.astype(str)

if (
    assets.isna().values.any()
):  # If has nan values, load the data from the fallback online database
    print(
        "WARNING: Asset data format is invalid or incomplete. "
        "Loading data from the remote backup database."
    )

    url = "https://raw.githubusercontent.com/Global-Data-Quantum/PortfolioData/1f3ebb95fdef245014a04d4273f688a2951e0061/data/qpo_tutorial_ibex.csv"
    assets = pd.read_csv(url, index_col=0)

qubo_settings = {
    "nt": 4,
    "nq": 2,
    "dt": 30,
    # maximum investment per asset is 2**nq/max_investment = 80%
    "max_investment": 5,
    "risk_aversion": 1000.0,
    "transaction_fee": 0.01,
    "restriction_coeff": 1.0,
}

optimizer_settings = {
    "de_optimizer_settings": {
        "num_generations": 20,
        "population_size": 40,
        "recombination": 0.4,
        "max_parallel_jobs": 5,
        "max_batchsize": 4,
        "mutation_range": [0.0, 0.25],
    },
    "optimizer": "differential_evolution",
    "primitive_settings": {
        "estimator_shots": 25_000,
        "estimator_precision": None,
        "sampler_shots": 100_000,
    },
}

ansatz_settings = {
    "ansatz": "optimized_real_amplitudes",
    "multiple_passmanager": False,
}

dpo_job = dpo_solver.run(
    assets=assets.to_dict(),  # Insert assets as a dictionary
    qubo_settings=qubo_settings,
    optimizer_settings=optimizer_settings,
    ansatz_settings=ansatz_settings,
    backend_name="ibm_torino",
    previous_session_id=[],
    apply_postprocess=True,
)

# Get the results of the job
dpo_result = dpo_job.result()

# Show the solution strategy
dpo_result["result"]

import pandas as pd

# Get results from the job
dpo_result = dpo_job.result()

# Convert metadata to a DataFrame, excluding 'session_id'
df = pd.DataFrame(dpo_result["metadata"]["all_samples_metrics"])

# Find the minimum objective cost
min_cost = df["objective_costs"].min()
print(f"Minimum Objective Cost Found: {min_cost:.2f}")

# Extract the row with the lowest cost
best_row = df[df["objective_costs"] == min_cost].iloc[0]

# Display the results associated with the best solution
print("Best Solution:")
print(f"  - Restriction Deviation: {best_row['rest_breaches']}%")
print(f"  - Sharpe Ratio: {best_row['sharpe_ratios']:.2f}")
print(f"  - Return: {best_row['returns']:.2f}")

import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
import matplotlib.patheffects as patheffects


def plot_normalized(dpo_x, dpo_y_normalized, random_x, random_y_normalized):
    """
    Plots normalized results for two sampling results.

    Parameters:
        dpo_x (array-like): X-values for the VQE Post-processed curve.
        dpo_y_normalized (array-like): Y-values (normalized) for the
        VQE Post-processed curve.
        random_x (array-like): X-values for the Noise (Random) curve.
        random_y_normalized (array-like): Y-values (normalized) for the
        Noise (Random) curve.
    """
    plt.figure(figsize=(6, 3))
    plt.tick_params(axis="both", which="major", labelsize=12)

    # Define custom colors
    colors = ["#4823E8", "#9AA4AD"]

    # Plot DPO results
    (line1,) = plt.plot(
        dpo_x, dpo_y_normalized, label="VQE Postprocessed", color=colors[0]
    )
    line1.set_path_effects(
        [patheffects.withStroke(linewidth=3, foreground="white")]
    )

    # Plot Random results
    (line2,) = plt.plot(
        random_x, random_y_normalized, label="Noise (Random)", color=colors[1]
    )
    line2.set_path_effects(
        [patheffects.withStroke(linewidth=3, foreground="white")]
    )

    # Set X-axis ticks to increment by 5 units
    plt.gca().xaxis.set_major_locator(MultipleLocator(5))

    # Axis labels and legend
    plt.xlabel("Objective cost", fontsize=14)
    plt.ylabel("Normalized Counts", fontsize=14)

    # Add DOCPLEX reference line
    plt.axvline(
        x=-4.11, color="black", linestyle="--", linewidth=1, label="DOCPlex"
    )  # DOCPlex value
    plt.ylim(bottom=0)

    plt.legend()

    # Adjust layout
    plt.tight_layout()
    plt.show()

import numpy as np
from collections import defaultdict

# ================================
# STEP 1: DPO COST DISTRIBUTION
# ================================

# Extract data from DPO results
counts_list = dpo_result["metadata"]["all_samples_metrics"][
    "counts"
]  # List of how many times each solution occurred
cost_list = dpo_result["metadata"]["all_samples_metrics"][
    "objective_costs"
]  # List of corresponding objective function values (costs)

# Round costs to one decimal and accumulate counts for each unique cost
dpo_counter = defaultdict(int)
for cost, count in zip(cost_list, counts_list):
    rounded_cost = round(cost, 1)
    dpo_counter[rounded_cost] += count

# Prepare data for plotting
dpo_x = sorted(dpo_counter.keys())  # Sorted list of cost values
dpo_y = [dpo_counter[c] for c in dpo_x]  # Corresponding counts

# Normalize the counts to the range [0, 1] for better comparison
dpo_min = min(dpo_y)
dpo_max = max(dpo_y)
dpo_y_normalized = [
    (count - dpo_min) / (dpo_max - dpo_min) for count in dpo_y
]

# ================================
# STEP 2: RANDOM COST DISTRIBUTION
# ================================

# Read the QUBO matrix
qubo = np.array(dpo_result["metadata"]["qubo"])

bitstring_length = qubo.shape[0]
num_random_samples = 100_000  # Number of random samples to generate
random_cost_counter = defaultdict(int)

# Generate random bitstrings and calculate their cost
for _ in range(num_random_samples):
    x = np.random.randint(0, 2, size=bitstring_length)
    cost = float(x @ qubo @ x.T)
    rounded_cost = round(cost, 1)
    random_cost_counter[rounded_cost] += 1

# Prepare random data for plotting
random_x = sorted(random_cost_counter.keys())
random_y = [random_cost_counter[c] for c in random_x]

# Normalize the random cost distribution
random_min = min(random_y)
random_max = max(random_y)
random_y_normalized = [
    (count - random_min) / (random_max - random_min) for count in random_y
]

# ================================
# STEP 3: PLOTTING
# ================================

plot_normalized(dpo_x, dpo_y_normalized, random_x, random_y_normalized)