Quantum Monte Carlo using Qiskit

Demo 1 — Quantum Random Number Generation

Why This Matters

Classical random number generators are pseudo-random — they use deterministic algorithms that appear random but follow a predictable pattern given the seed. Quantum mechanics provides true randomness from the fundamental nature of quantum measurement.

In finance, true randomness is valuable for:

• Monte Carlo simulation integrity
• Cryptographic key generation
• Fair lottery / selection mechanisms

The Quantum Circuit

The simplest QRNG uses a Hadamard gate to create a perfect 50/50 superposition, then measures:

How it works:

1. Start with qubit in state $|0\rangle$
2. Apply Hadamard (H) gate → creates equal superposition: $(|0\rangle + |1\rangle)/\sqrt{2}$
3. Measure → exactly 50% chance of 0, 50% chance of 1
4. This randomness is guaranteed by quantum mechanics

Python Code: Quantum Random Number Generator

Python — Qiskitfrom qiskit import QuantumCircuit
from qiskit_aer import AerSimulator

# ─── Create Quantum Circuit ───
def quantum_random_bits(n_bits=8, shots=1):
    """Generate truly random bits using quantum mechanics."""
    qc = QuantumCircuit(n_bits, n_bits)

    # Apply Hadamard to each qubit → creates superposition
    for i in range(n_bits):
        qc.h(i)

    # Measure all qubits
    qc.measure(range(n_bits), range(n_bits))
    return qc

# Create and visualize the circuit
qc = quantum_random_bits(8)
print("Quantum Random Number Generator Circuit:")
print(qc.draw(output='text'))

# ─── Run on Simulator ───
simulator = AerSimulator()
n_random_numbers = 1000
qc = quantum_random_bits(8, shots=n_random_numbers)

result = simulator.run(qc, shots=n_random_numbers).result()
counts = result.get_counts()

# Convert to random numbers
random_numbers = []
for bitstring, count in counts.items():
    value = int(bitstring, 2)  # Convert binary to decimal
    random_numbers.extend([value] * count)

import numpy as np
import matplotlib.pyplot as plt

print(f"\nGenerated {len(random_numbers)} quantum random numbers")
print(f"Range: [{min(random_numbers)}, {max(random_numbers)}]")
print(f"Mean: {np.mean(random_numbers):.2f} (expected: 127.5)")
print(f"Std: {np.std(random_numbers):.2f}")

# ─── Visualization ───
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].hist(random_numbers, bins=32, density=True, alpha=0.7,
             color='steelblue', edgecolor='white')
axes[0].set_title('Quantum Random Number Distribution', fontweight='bold')
axes[0].set_xlabel('Value (0-255)')
axes[0].set_ylabel('Density')
axes[0].axhline(y=1/256, color='red', linestyle='--', label='Uniform (expected)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

top_10 = dict(list(sorted(counts.items(), key=lambda x: -x[1]))[:10])
axes[1].bar(range(len(top_10)), list(top_10.values()),
            color='darkorange', alpha=0.7)
axes[1].set_xticks(range(len(top_10)))
axes[1].set_xticklabels(list(top_10.keys()), rotation=45, fontsize=8)
axes[1].set_title('Top 10 Most Frequent Outcomes', fontweight='bold')
axes[1].set_xlabel('Bitstring')
axes[1].set_ylabel('Count')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Demo 2 — Quantum Amplitude Estimation: Convergence Comparison

Estimating Expected Payoff of a European Call Option

Pythonimport numpy as np
import matplotlib.pyplot as plt

# Problem: Estimate expected payoff of a European Call Option
S0 = 100       # Current stock price
K = 105        # Strike price
T = 1.0        # Maturity
r = 0.05       # Risk-free rate
sigma = 0.20   # Volatility

# ─── Classical Monte Carlo ───
def classical_monte_carlo_call(n_sims):
    """Classical Monte Carlo pricing of European Call."""
    Z = np.random.standard_normal(n_sims)
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
    payoff = np.maximum(ST - K, 0)
    price = np.exp(-r * T) * np.mean(payoff)
    std_error = np.exp(-r * T) * np.std(payoff) / np.sqrt(n_sims)
    return price, std_error

# ─── Convergence Analysis ───
sim_counts = [10, 50, 100, 500, 1000, 5000, 10000, 50000, 100000]
classical_prices = []
classical_errors = []

np.random.seed(42)
for n in sim_counts:
    price, error = classical_monte_carlo_call(n)
    classical_prices.append(price)
    classical_errors.append(error)

# Theoretical quantum errors (quadratic speedup)
quantum_errors_theoretical = [
    classical_errors[-1] * (sim_counts[-1] / n)
    for n in sim_counts
]

# ─── Visualization ───
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Price convergence
axes[0].errorbar(sim_counts, classical_prices, yerr=classical_errors,
                 fmt='o-', color='steelblue', capsize=3,
                 label='Classical MC Price ± SE')
axes[0].axhline(y=classical_prices[-1], color='red', linestyle='--',
                alpha=0.5, label=f'Best Estimate: ₹{classical_prices[-1]:.2f}')
axes[0].set_xscale('log')
axes[0].set_title('Classical MC: Price Convergence', fontweight='bold')
axes[0].set_xlabel('Number of Simulations (log scale)')
axes[0].set_ylabel('Estimated Call Price (₹)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Error comparison
axes[1].loglog(sim_counts, classical_errors, 'o-', color='steelblue',
               label='Classical MC: O(1/√N)', linewidth=2)
axes[1].loglog(sim_counts, quantum_errors_theoretical, 's--',
               color='darkorange',
               label='Quantum AE: O(1/N) (theoretical)', linewidth=2)
axes[1].set_title('Error: Classical vs Quantum', fontweight='bold')
axes[1].set_xlabel('Samples / Oracle Calls (log scale)')
axes[1].set_ylabel('Estimation Error (log scale)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].annotate('Quadratic\nSpeedup!',
    xy=(5000, quantum_errors_theoretical[5]),
    xytext=(500, quantum_errors_theoretical[3]),
    arrowprops=dict(arrowstyle='->', color='red'),
    fontsize=12, color='red', fontweight='bold')

plt.tight_layout()
plt.show()

Summary Table

Python# ─── Summary Table ───
print("\n" + "=" * 70)
print(f"{'CONVERGENCE COMPARISON: Classical MC vs Quantum AE':^70}")
print("=" * 70)
print(f"{'Samples/Calls':<15} {'Classical Error':<20} {'Quantum Error (theory)':<25}")
print("-" * 70)
for n, ce, qe in zip(sim_counts, classical_errors, quantum_errors_theoretical):
    print(f"{n:<15,} {ce:<20.6f} {qe:<25.6f}")
print("=" * 70)

Demo 3 — Qiskit Probability Estimation

Estimate the probability of a biased coin using quantum amplitude estimation concepts:

Python — Qiskit# ─── Simple Quantum Circuit for Probability Estimation ───
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
import numpy as np

def create_bernoulli_oracle(p):
    """
    Create a quantum oracle that prepares a state with probability p
    of measuring |1⟩. This simulates a Bernoulli random variable.
    In finance: P(loss > VaR threshold)
    """
    theta = 2 * np.arcsin(np.sqrt(p))
    qc = QuantumCircuit(1)
    qc.ry(theta, 0)
    return qc

# Create an oracle for p = 0.3 (30% probability of "success")
p_true = 0.3
oracle = create_bernoulli_oracle(p_true)

print("Oracle Circuit (encodes probability):")
print(oracle.draw(output='text'))

# ─── Classical Estimation: Sample Many Times ───
n_shots_list = [10, 100, 1000, 10000]

print(f"\nTrue probability: {p_true}")
print(f"\n{'Shots':<10} {'Estimated p':<15} {'Error':<10}")
print("-" * 35)

simulator = AerSimulator()
for n_shots in n_shots_list:
    qc = QuantumCircuit(1, 1)
    qc.ry(2 * np.arcsin(np.sqrt(p_true)), 0)
    qc.measure(0, 0)

    result = simulator.run(qc, shots=n_shots).result()
    counts = result.get_counts()

    p_estimated = counts.get('1', 0) / n_shots
    error = abs(p_estimated - p_true)
    print(f"{n_shots:<10} {p_estimated:<15.4f} {error:<10.4f}")

print("\nNote: Classical sampling error decreases as O(1/√N)")
print("Quantum amplitude estimation would achieve O(1/N)")

Presenter Guidance