{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Credit Risk Analysis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Introduction\n",
"This tutorial shows how quantum algorithms can be used for credit risk analysis.\n",
"More precisely, how Quantum Amplitude Estimation (QAE) can be used to estimate risk measures with a quadratic speed-up over classical Monte Carlo simulation.\n",
"The tutorial is based on the following papers:\n",
"\n",
"- [Quantum Risk Analysis. Stefan Woerner, Daniel J. Egger.](https://www.nature.com/articles/s41534-019-0130-6) [Woerner2019]\n",
"- [Credit Risk Analysis using Quantum Computers. Egger et al. (2019)](https://arxiv.org/abs/1907.03044) [Egger2019]\n",
"\n",
"A general introduction to QAE can be found in the following paper:\n",
"\n",
"- [Quantum Amplitude Amplification and Estimation. Gilles Brassard et al.](http://arxiv.org/abs/quant-ph/0005055)\n",
"\n",
"The structure of the tutorial is as follows:\n",
"\n",
"1. [Problem Definition](#Problem-Definition)\n",
"2. [Uncertainty Model](#Uncertainty-Model)\n",
"3. [Expected Loss](#Expected-Loss)\n",
"4. [Cumulative Distribution Function](#Cumulative-Distribution-Function)\n",
"5. [Value at Risk](#Value-at-Risk)\n",
"6. [Conditional Value at Risk](#Conditional-Value-at-Risk)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from qiskit import QuantumRegister, QuantumCircuit\n",
"from qiskit.circuit.library import IntegerComparator\n",
"from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem\n",
"from qiskit_aer.primitives import Sampler"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem Definition\n",
"\n",
"In this tutorial we want to analyze the credit risk of a portfolio of $K$ assets.\n",
"The default probability of every asset $k$ follows a *Gaussian Conditional Independence* model, i.e., given a value $z$ sampled from a latent random variable $Z$ following a standard normal distribution, the default probability of asset $k$ is given by\n",
"\n",
"$$p_k(z) = F\\left( \\frac{F^{-1}(p_k^0) - \\sqrt{\\rho_k}z}{\\sqrt{1 - \\rho_k}} \\right) $$\n",
"\n",
"where $F$ denotes the cumulative distribution function of $Z$, $p_k^0$ is the default probability of asset $k$ for $z=0$ and $\\rho_k$ is the sensitivity of the default probability of asset $k$ with respect to $Z$. Thus, given a concrete realization of $Z$ the individual default events are assumed to be independent from each other.\n",
"\n",
"We are interested in analyzing risk measures of the total loss\n",
"\n",
"$$ L = \\sum_{k=1}^K \\lambda_k X_k(Z) $$\n",
"\n",
"where $\\lambda_k$ denotes the _loss given default_ of asset $k$, and given $Z$, $X_k(Z)$ denotes a Bernoulli variable representing the default event of asset $k$. More precisely, we are interested in the expected value $\\mathbb{E}[L]$, the Value at Risk (VaR) of $L$ and the Conditional Value at Risk of $L$ (also called Expected Shortfall). Where VaR and CVaR are defined as\n",
"\n",
"$$ \\text{VaR}_{\\alpha}(L) = \\inf \\{ x \\mid \\mathbb{P}[L <= x] \\geq 1 - \\alpha \\}$$\n",
"\n",
"with confidence level $\\alpha \\in [0, 1]$, and\n",
"\n",
"$$ \\text{CVaR}_{\\alpha}(L) = \\mathbb{E}[ L \\mid L \\geq \\text{VaR}_{\\alpha}(L) ].$$\n",
"\n",
"For more details on the considered model, see, e.g.,
\n",
"[Regulatory Capital Modeling for Credit Risk. Marek Rutkowski, Silvio Tarca](https://arxiv.org/abs/1412.1183)\n",
"\n",
"\n",
"\n",
"The problem is defined by the following parameters:\n",
"\n",
"- number of qubits used to represent $Z$, denoted by $n_z$\n",
"- truncation value for $Z$, denoted by $z_{\\text{max}}$, i.e., Z is assumed to take $2^{n_z}$ equidistant values in $\\{-z_{max}, ..., +z_{max}\\}$ \n",
"- the base default probabilities for each asset $p_0^k \\in (0, 1)$, $k=1, ..., K$\n",
"- sensitivities of the default probabilities with respect to $Z$, denoted by $\\rho_k \\in [0, 1)$\n",
"- loss given default for asset $k$, denoted by $\\lambda_k$\n",
"- confidence level for VaR / CVaR $\\alpha \\in [0, 1]$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# set problem parameters\n",
"n_z = 2\n",
"z_max = 2\n",
"z_values = np.linspace(-z_max, z_max, 2**n_z)\n",
"p_zeros = [0.15, 0.25]\n",
"rhos = [0.1, 0.05]\n",
"lgd = [1, 2]\n",
"K = len(p_zeros)\n",
"alpha = 0.05"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Uncertainty Model\n",
"\n",
"We now construct a circuit that loads the uncertainty model. This can be achieved by creating a quantum state in a register of $n_z$ qubits that represents $Z$ following a standard normal distribution. This state is then used to control single qubit Y-rotations on a second qubit register of $K$ qubits, where a $|1\\rangle$ state of qubit $k$ represents the default event of asset $k$. The resulting quantum state can be written as\n",
"\n",
"$$ |\\Psi\\rangle = \\sum_{i=0}^{2^{n_z}-1} \\sqrt{p_z^i} |z_i \\rangle \\bigotimes_{k=1}^K \n",
"\\left( \\sqrt{1 - p_k(z_i)}|0\\rangle + \\sqrt{p_k(z_i)}|1\\rangle\\right),\n",
"$$\n",
"\n",
"where we denote by $z_i$ the $i$-th value of the discretized and truncated $Z$ [Egger2019]."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from qiskit_finance.circuit.library import GaussianConditionalIndependenceModel as GCI\n",
"\n",
"u = GCI(n_z, z_max, p_zeros, rhos)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
┌───────┐\n", "q_0: ┤0 ├\n", " │ │\n", "q_1: ┤1 ├\n", " │ P(X) │\n", "q_2: ┤2 ├\n", " │ │\n", "q_3: ┤3 ├\n", " └───────┘" ], "text/plain": [ " ┌───────┐\n", "q_0: ┤0 ├\n", " │ │\n", "q_1: ┤1 ├\n", " │ P(X) │\n", "q_2: ┤2 ├\n", " │ │\n", "q_3: ┤3 ├\n", " └───────┘" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u.draw()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now use the simulator to validate the circuit that constructs $|\\Psi\\rangle$ and compute the corresponding exact values for\n", "\n", "- expected loss $\\mathbb{E}[L]$\n", "- PDF and CDF of $L$ \n", "- value at risk $VaR(L)$ and corresponding probability\n", "- conditional value at risk $CVaR(L)$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "u_measure = u.measure_all(inplace=False)\n", "sampler = Sampler()\n", "job = sampler.run(u_measure)\n", "binary_probabilities = job.result().quasi_dists[0].binary_probabilities()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# analyze uncertainty circuit and determine exact solutions\n", "p_z = np.zeros(2**n_z)\n", "p_default = np.zeros(K)\n", "values = []\n", "probabilities = []\n", "num_qubits = u.num_qubits\n", "\n", "for i, prob in binary_probabilities.items():\n", " # extract value of Z and corresponding probability\n", " i_normal = int(i[-n_z:], 2)\n", " p_z[i_normal] += prob\n", "\n", " # determine overall default probability for k\n", " loss = 0\n", " for k in range(K):\n", " if i[K - k - 1] == \"1\":\n", " p_default[k] += prob\n", " loss += lgd[k]\n", "\n", " values += [loss]\n", " probabilities += [prob]\n", "\n", "\n", "values = np.array(values)\n", "probabilities = np.array(probabilities)\n", "\n", "expected_loss = np.dot(values, probabilities)\n", "losses = np.sort(np.unique(values))\n", "pdf = np.zeros(len(losses))\n", "for i, v in enumerate(losses):\n", " pdf[i] += sum(probabilities[values == v])\n", "cdf = np.cumsum(pdf)\n", "\n", "i_var = np.argmax(cdf >= 1 - alpha)\n", "exact_var = losses[i_var]\n", "exact_cvar = np.dot(pdf[(i_var + 1) :], losses[(i_var + 1) :]) / sum(pdf[(i_var + 1) :])" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Expected Loss E[L]: 0.6641\n", "Value at Risk VaR[L]: 2.0000\n", "P[L <= VaR[L]]: 0.9521\n", "Conditional Value at Risk CVaR[L]: 3.0000\n" ] } ], "source": [ "print(\"Expected Loss E[L]: %.4f\" % expected_loss)\n", "print(\"Value at Risk VaR[L]: %.4f\" % exact_var)\n", "print(\"P[L <= VaR[L]]: %.4f\" % cdf[exact_var])\n", "print(\"Conditional Value at Risk CVaR[L]: %.4f\" % exact_cvar)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [ "nbsphinx-thumbnail" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAElCAYAAAAhjw8JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+rUlEQVR4nO3dd5wV1fn48c9DZ+lFNCqKqBRF5atIBAsoYEQhRERgoybYiBo12MsPFUuMBZAkVogJiGVBhRilKSSLBTWCkgQpVqQIUgQpuyjl+f1x5uLl7i1z6+zd+7xfr3ndvWfOmfvszOw9O3PmnCOqijHGGJOOakEHYIwxJv9ZZWKMMSZtVpkYY4xJm1Umxhhj0maViTHGmLRZZWKMMSZtVpkYkwYRGS8iKiKtAvr8Ed7nd49IVxEpDSKmsBgC3Tcmt6wyMSnzvijyvqOSiCwP/S7eslNENorI/0RkooicLyK1svTZpfm6D2NVZKYw1Qg6AGMqkT8Cm3H/ZDUE2gLnAhcCn4rIhar674gytwEPAKtzGGe4R4ESYEVAnx9P0PvG5JBVJsb8aIyqLg9PEJFGwL3ANcDrInKSqi4NrVfVNcCanEYZRlU3ABuC+vx4gt43JrfsNpfJCRGpLSK3ereOykRki4i8JSIDY+T/uYjMEZE1IvK9iHwtInNF5KqIfK1FZKyIfCYi5SLyrfcZT4pIs3TjVtXvVPVa4BmgEe4/7fDPj9oukCh+EWnl3d7q5r0Pv81WGrad5d7SUERGez/vFJER3vq4t5pE5EDvVt06b/8sEJFfRsk3xNvOkBjbqRAXcJf39l/h8SfaN966gSLypoh858X1PxG5TURqR8kb2gf1RORhEVnh7dPPROQWEZFoMZvcsisTk3Vee8Ms3BfnUuAxoAgYAEwSkY6qentY/qHAU8Ba4FXcf94tgGOBi4HHvXw/AT7A3ZKaDrwM1AEOAy7C3QLamKFf4x7gV0AfEWmoqlvi/L5+4t8M3A0MAQ71fg5ZHrHJWsA/gabA68AW4EsfMTcB5nmf9TegMTAQeE5EDlLVh31sI5YxwC9wx3RClJhjEpH7cbfANgDPA9uA3sD9wM9E5ExV/SGiWE3cOXQgMAPY5X3+A7hjfjcmWKpqiy0pLYC6Uyhhvtu8vNOBGmHpLXBfQgp0DUtfAHwPtIiyreZhP1/jlf1dlHz1gLo+f49QDK0S5Fvp5Ts9LG18ZFm/8XvvS+Ptw7DYZgP1oqwf4a3vHu3YAJOBamHphwHfAj8ArcPSh3j5h8Q51qV+PjvBvunipa0ADghLr4GreBW4PcY+mB5+TL3zZ7O31Az676HQF7vNZXLhEtyXwfWquiuUqKrrcO0RAJdFlNkF7IzckLo2gkjlUfJtV9UK6WkKNSTv5yNvMvH7cYOqbk+yzG7gFlXdE/b5XwJ/wv2nf1GKsaTjEu/1PlVdGxbXLuAGYA8Vz4WQa8OPqXf+vIK7/dg2O+Eav6wyMVklIg2AI4CvNazhOsw/vdf/C0t7DncbbLGIPCIivxCRaF/g/8DdInlMRF4WkaEicnQW76GHtpvoUV6/8fu1A/hvCuVWeJVHpFLv9f+irMu2473Xf0auUNVPgFXAYd6DD+G+U9XPomxvpffaJHMhmlRYZWKyLfSlEOupnlB641CCqo4Gfg18BVwLTAW+EZF/iUinsHxfAZ2BKUBPXDvFIuArEbk2g79DyIHe6/p4mfzGn4R16t3XSdI3MdJDVwSRX9i5kPT54NkcI3/oSrd66iGZTLDKxGTbd97rATHW/yQiHwCq+oyqngQ0A84BngZOA2aF/5evqktUdZCXrxNwK+68/qOIXJqpX0JEjgAOxn15LUiU32/8PqXaqXH/GOmhYxG+z0O3wio8lCMijVP8/GhSOh9M5WeVickqVd0KfA4cJCJHRslyuvf6YYzym1V1uqpejmvQbYr7Uo7Mt0tVF6jqg0Cxl/yLNMMPd6f3+qr3O/niI/7dACKSjf+sD4n2WC7Q3Xv9KCxtk/faMkr+WFdTu73XZGIPfWb3yBVhFfaXqro5iW2aSsAqE5MLf8W1Nzwc/qUpIs2BO8LyhNJPj9Hu0cJ7LfPynRDl3jr8+B95WbqBe/07/oRrrN6Mu/JJVMZX/J7Qo8uHpBNnDNWBB0Vk79+5iByGu/W2C3g2LO983NXJL0WkKCx/U+ChGNtPJfbQcR4efoXmnRcjcd9JTyexPVNJWD8TkzYRGR9n9VW4L4neQD/gPyIyHddAfT7uC/YhVX07rMxUYJuIvId7LFSAU4ETcbeYZnv5LgJ+IyJv465+NgGHA31xj+aOSfJXGSYim73PCw2nchruMeNPgAu9RuJE/MYPMAe3H6Z4+6Uc+EpVJyYZezT/BX4KLBCR1/mxn0lj4GZV/TyUUVXXiMhzuH26UESm4fbB2cCbRG+s/xeuAvqDiHTAu7pR1ftiBaSq80TkIeBmYJGIvARsx50fHYC3gXT6v5igBP1ssi35u/BjX4Z4S2Mvbx3gdlwDeTmwFffFURxlu1fgvpC/wP0X/y3u9sjNQIOwfD8FngD+4+UpBz7DddDrkMTvsTwi5p3e9v4HTMR1rqwVo+x4Kval8BW/l7c6rrPeF97n7tOfw4tteZzYRxC7n0kp7qGBZ4F1uKfCPgR+GWNbtXFf5Ktw/VA+w/URqhEZV1iZC4GF3r7fp99RtH0Ttm6wd/y3enF9DPw/oE6M4xN1H8T6/W3J/SLeATHGGGNSZm0mxhhj0maViTHGmLRZZWKMMSZtVpkYY4xJW0E8Gty8eXNt1apV0GEYz7KNywBo28zG5qv0trhjRUM7VilZ5u2/tvm5/xYsWLBBVX2N2FAQlUmrVq2YP39+0GEYT/fx3QEoHVIaaBzGh9nd3WvP0iCjyF/du7vX0tIgo0iZiHzlN6/d5jLGGJO2grgyMZXL8NOGBx2C8auDHau0DC+c/WeVicm5nq17Bh2C8esAO1Zp6Vk4+89uc5mcW7h2IQvXLgw6DOPHpoVuMalZuNAtBcCuTEzODZs5DLAG+LywYJh7tQb41Awb5l7ztAE+GXZlYowxJm12ZWKMSdmWLVtYt24dO3fuDDqUyumuu9zrkiXBxhGhZs2atGjRgoYNG2Zsm1aZGGNSsmXLFr755hsOOugg6tatS/T5wApcNe/mTyXqtKiqlJeXs3r1aoCMVSh2m8sYk5J169Zx0EEHUVRUZBVJHhERioqKOOigg1i3bl3GtmtXJgm0unVa0CEEavkD52R8m/f3uD/j2zRZclzsY7Vz507q1q2bw2Dy0EEHBR1BTHXr1s3o7UmrTEzOdW3ZNegQjF/7xT9WdkWSQP36QUcQU6aPnd3mMjk3b+U85q2cF3QYxo/189xiUrNtm1sKgF2ZmJy7fc7tgPUzyQv/ccfK+pmkyGvkrkwN8NliVybGmII2YsQIRCTq8uyzz1JaWrr3fePGjfeWW758OSLCa6+9FnPb3S+6CGnXDhHh0UcfzcFvExy7MjHGFLxGjRoxc+bMCulHHHEEixYtAuC5556jTZs2SW338bvuYsu2bXQZPDgjcVZmVpkYYwpejRo1OOmkk+LmOfbYY+nQoUNS2z3qiCPSCSuv2G0uY4wxacv5lYmIHAX8GegCbAb+Atytqrt9lO0P3AZ0AMqAD4DzVHV71gI2GTfmrDFBh2D8OmFM0BHkzK5duyqk1aiR5ldky5bplc8jOa1MRKQJMBtYDPQDDgdG4a6Q4s4iIyKXAY8CDwE3AU2AM7BbdXmn4wEdgw7B+NWkY9JFQtMyhxt49ECuOvEqynaWcfZzZ1dYP6TjEIZ0HMKGsg0MmDygwvorO13JoA6DWPndSi6aelGF9Td0uYG+bfsmHWvIxo0bqVmzZoX0L7/8MuVtAlBUlF75PJLrL+IrgLpAf1XdArwhIg2BESLykJdWgYg0Bx4BrlHVcWGrpmY9YpNxs7+YDdgkWXlhrTtWVX2SrEaNGjF79uwK6QceeCDLly9PfcNbon6lVUm5rkx6A7MiKo0S4EGgG/BqjHIDvdcJWYzN5Mh9b94HWGWSFxa5Y5VMZRKv/1BRzaK465sXNY+7vmWjllnpn1SjRg06deqU8e2yZk3mt1lJ5boBvh2wNDxBVVfg2j/axSn3U2AZcKmIrBKRnSLyvojYuBzGGFMJ5LoyaYJrdI+0yVsXywFAW1y7yi1AX2A7MFNE9o9WQESGish8EZm/fv36tII2xhgTX740XgtQHzhfVWcCiMg84CvgauCOyAKqOhYYC9CpUyfNXajGmHyza9cu3nvvvQrpLX08jfXOO++wY8eOfdJatWqVndtmlViuK5NNQKMo6U28dfHKKVAaSlDVLSKyADgqkwEaYwrPd999R5cuXSqk33vvvZxyyilxyz7wwAMV0n79618zfvz4TIWXF3JdmSwlom1ERFoCRUS0pURYgrs6iRwzWYA9mQzQZN9TfZ4KOgTjV+eqf6xGjBjBiBEjYq4vLS0FYPfu3ezevZvq1asD7upDNf5Nj90HH5wwT1WR6zaTGcDPRKRBWNogoByYG6dcaCS100MJItIIOAH4T6aDNNnVtnlb2jav+qOoVgkN27rF0LFjR5o1a5ZUmR7nnEPNBg0SZ6wCcn1l8iRwLTBFRB4EWgMjgNHhjwuLyGfAXFW9FEBV54vIK8DTInIrsAG4GdgJPJbbX8Gk69Vl7gnwdDqZmRxZ5T2tf3DhHqsTTjiBDz74AEi+R/xTDz/M1m3boEEDDj300GyEV2nktDJR1U0i0gPXk/1V3JNdj+AqlMi4qkekXQg8DIzG3RZ7BzhDVeO1tZhKaNS7owCrTPLCUnesCrkyadCgQcqN6W0bNoSGDQtiPpOcP82lqotxw6DEy9MqSto24EpvMcYYU4nYqMHGGGPSZpWJMcaYtFllYowxJm350gPeVCETz50YdAjGry52rNJy2GFBR5AzVpmYnGvZqHAmDMp79exYpaVWraAjyBm7zWVybtKiSUxaNCnoMIwfX01yi0nNt9+6pQBYZWJy7on5T/DE/CeCDsP48ekTbqmC+vbtyzHHHBNz/dVXX03jxo35/vvv425n+fLliMjepX79+hx33HH85S9/gfXr3RJm9erVNGjQgM8//3xvmojw6KOPxvyMPn36cO+99/r8zYJhlYkxpiAVFxezaNEiFi9eXGHd7t27eemll+jfvz+1a9f2tb2RI0fy7rvvMnXqVI477jguv/xynv3HPyrku+++++jTpw+HH36471hvueUWRo8ezebNm32XyTWrTIwxBalfv34UFRXxwgsvVFj3r3/9i2+++Ybi4mLf22vbti0nnXQSvXr1YsKECbRv355n/v73ffJs2bKFCRMmcMkllyQV66mnnkqzZs2YOLHyPhBhlYkxpiDVq1ePvn37MmlSxTahkpISWrRowYEHHsjgwYNp2bIlRUVFHH300YwZM4Y9e+IPVi4iHHPMMaxcu3af9MmTJ1O3bl3OOCPuICBRnXfeeTzzzDNJl8sVe5rLGJNZs7tXTDtkILS5CnaVQenZFde3HuKWHRvg7QEV1x95JRw6CLavhHcvqri+3Q0pjR9WXFzMpEmTWLBgASeccAIAO3fuZMqUKVxwwQWsXbuWtm3bcsEFF9CgQQMWLlzIXXfdRXl5Obfddlvcba9YsYLDDj54n7Q5c+bQuXPnvcPYJ6Nr1648/PDDbNq0iSZN4k1MGwyrTEzOvTTwpaBDMH6dUrWPVe/evWncuDElJSV7K5NZs2axadMmiouL6dq1Kz169ABAVTnllFMoKytj3LhxFSqTPXv2sGvXLrZu3cqECRP48MMPeWPGDGjdem+eBQsW0K9fv5RiPe6441BV5s+fT69evVL8jbPHKhOTc82LmgcdgvGrTgrHqmdp7HU1iuKvr9M8/vp6LeOvT1KtWrXo378/kydP5qGHHkJEmDRpEoceeihdunRhx44d/OEPf+C5555jxYoV7Ny5c2/ZXbt27TMkfWQl8cc//pHTIm5nrV27lubNUzv/Q+XWRtw6qyyszcTk3PiF4xm/cHzQYRg/vhjvliqsuLiYFStW8O6777Jjxw5eeeUVBg8ejIhwyy23MHLkSIYOHcr06dP54IMPGD58OECFed8feeQRPvjgA6ZNm0bXrl258cYb+U9pKWzYsDfPjh07fD8dFilULvJzKwu7MjE5F6pIhnQcEmgcxodQRdJ6SJBRZNXpp5/O/vvvT0lJCWvWrGHr1q17n+J68cUXueaaa7j55pv35p82bVrU7RxxxBF75z3p0qULRx55JLfeeSczxo0D76qiadOmKT/eGyrXtGnTlMpnm12ZGGMKWvXq1Rk4cCAvvvgizz//PO3bt+e4444DoLy8fJ8rid27d1NSUpJwm02aNOGWW25h5ltv8d9ly/amt23bli+//DKlOJcvXw5AmzZtUiqfbVaZGGMKXnFxMWvXrmXq1Kn79C3p1asXjz32GBMnTmTatGn07ds3YY/4kCuvvJKmjRrx8NNP7007+eSTWbBgQdT8Cxcu5KWXXtpnmTt37t718+fPp1GjRhx99NEp/pbZZZWJMabgdenShVatWqGq+1Qmf/7znzn11FP57W9/yyWXXEKHDh0SPhIcUr9+fX73q19RMn06K1euBKB///4sXryYFStWVMj/9NNPc/755++z3HXXXXvXz5w5k3PPPZdq1Srn17a1mRhjDES9/bT//vszderUCumXX3753p9DlVA0d/72t9z5299CSzf6cqdOnTjmmGOYNGkSN9100958scqHfPfdd8yaNYvZs2f7+l2CkPMqTkSOEpE5IlImIl+LyD0iErcHj4i0EhGNsiS+eWkqnekXTGf6BdODDsP40X26W0xqjjjCLWGGDx/OY489xq5du3xv5oknnuCkk07ilFNOyXSEGZPTKxMRaQLMBhYD/YDDgVG4Sm24j03cCLwT9n5DrIym8iqqWRR0CMavGnas0hKlp/uAAQP44osvWL16NYceeqivzTRq1Ig//elPmY4uo3J9m+sKoC7QX1W3AG+ISENghIg85KXFs0xV38t6lCarHv/gcQCuOvGqgCMxCX3ijhVt7FilZN0699qixd6kUP+VZFx55ZWZjCorcn2bqzcwK6LSKMFVMN1yHIsJyOSPJzP548lBh2H8WDHZLSY1mza5pQDkujJpBywNT1DVFUCZty6Rv4nIbhFZIyKjRaRuNoI0xhiTnKRuc4lILWA/oA7wraomW+U2ATZHSd/krYvle+Ax4HVgC9AduAXX5hJ11DQRGQoMBTjkkEOSDNMYY0wyElYmInI08CugJ3AMUD1s3UZgHvAS8LKqlmcjSFVdA1wdllQqIt8Aj4vIcar6nyhlxgJjATp16hT/uTtjjDFpiXmbS0ROFpF/Af/FtWfMBS4Hfg78DBgEPATsAB4BvhaRO0SkfpzP2wQ0ipLexFuXjNDY2CckWc4YY0yGxbsymQL8CbhIVVfF24jXT6QnMMxLujdG1qVEtI2ISEugiIi2FB804tXkidIhpUGHYPzK4HDvBalt26AjyJl4DfCHqurvE1UkAKq6W1VnqWpv4OE4WWcAPxORBmFpg4By3JVPMkLTsUUf6MYYY5Lw8ssvc8YZZ9C4cWNq165NmzZtuP7663nnnXcQEV5++eWo5b755htq1KjBgw8+6OtzhgwZgoggIlSrVo2DDz6Y4uLivQM5Rrr22mu5+OKL974fMWJE3DlR5s+fT9OmTfnuu+98xZMpMSsTVU1p0PwE5Z7ENaZPEZGeXiP5CGB0+OPCIvKZiDwd9n6EiIwSkf5euXtwt9amqOp/U4nTBGfkvJGMnDcy6DCMH0tGuqWKu+GGGxg4cCCtW7dm4sSJvP7661x33XXMmTOHkSNHcuSRR8YcLfjFF19kz549DB48uOLKtWvdEqFdu3a8++67vP3229xzzz2UlpZy9tln88MPP+yTb+XKlYwbNy6pfimdOnXi//7v/3jkkUd8l8mEpDstikgN3FNSpwMC/At4SlUTjg2gqptEpAfwKPAq7smuR3AVSmRc4V1Hl+J6v1+G65OyAncF9Ptk4zfBe+2T1wC4seuNAUdiElrtjhXtq+6xevXVVxk9ejRPP/00l1xyyd70bt26MXToUF5//XXee+89Hn74YbZt20b9+vs2C5eUlNClS5fovdlDVwcHHLBPcr169TjppJMAN7d7UVERxcXFzJ8/n65du+7N9+STT3L88cfTrp2fnhM/uvjii7nxxhsZPnz4PrNBZlMq/Uz+iHu6qxT4N3Az7rFdX1R1saqeoap1VfUnqnqHqu6OyNNKVYeEvS9R1U6q2khVa6nqEap6p6r6GwvaGGNieOSRRzj++OP3qUhCqlevTu/evSkuLqa8vJxXXnlln/UrV65k3rx5e0caHjVqFCeeeCKNGjVi//33p+8VV/DZV18ljCE0f0podOGQZ555hgEDBkQrEtfPf/5zvv32W2bNmpV02VTFe5rr+BirBgBnqupjqvoQcBVwfjaCM8aYbNq5cyfz5s3jrLPOipuvXbt2dOzYscKtrkmTJlGtWjXOP999Ba5atYqrr76aV155hXHjxrF7zx66FhcnbL8IDUl/2GGH7U1btmwZq1at2udKxa+GDRty9NFH53SU4XjXPzNE5B/A7aq6Pix9LdALeFlEBHe7a00WYzTG5JPu3SumDRwIV10FZWVw9tkV1w8Z4pYNGyDaf+JXXgmDBsHKlXDRRRXX33AD9O2bdKgbN27k+++/99Wxubi4mDvuuINNmzbRpInrY11SUsIZZ5zB/vvvD7BPO8Xu3bvpdeihtOjalVdeeYVf/epX+2xv165dqCpLlizh1ltv5ayzzqJz585714cm0erQoUPSvxe4q51///vfKZVNRbzbXG1xw5wsEZGbRKSml3418GcRWYfrG/JLoPKPQmYqjbo161K3po2Ekxeq13VLFef+L45v8ODB7Ny5c+/8Jp9//jkLFizYZzKt9957j169etGsWTNq1KhBUceObCsr45NPPtlnWwsWLKBmzZrUqlWL4447ji1btvDCCy/sk2ft2rXUqVOHevXqpfQ7NW/enLVRGv+zJeaViapuBn4nIk/ihokfKiI3qOo/RKQVP/YXWWZtFyYZMy6YEXQIxq/TUzhWpaWx1xUVxV/fvHn89S1bxl+fpGbNmlG7du2oMx9GOuSQQ+jatSslJSVccskllJSUULt2bfr37w+4W1VnnnkmnTt35qmnnuLAAw+kVq1anHPOOezYse9Dru3bt+eZZ55h586dvP3229x+++385je/YdKkSXvz7NixY5/555NVu3btCp+bTQmb+VV1CXC2iJwDjBSRa4Df2SO5xph8V7NmTU4++WRmzZrFfffdlzB/cXExv/vd71i/fj0lJSX07t2bRo3coB4zZ86krKyMV155Ze/VxK5du/j2228rbKeoqIhOnToBbsrgHTt2cOedd3L99dfz05/+FICmTZuyZcsW9uzZk9JUvZs3b6Zp06ZJl0tVwghFpK6INFLVabixuWYAc0XkURHJXaSmyrh37r3cOzfWIAmmUvnfvW6pwoYNG8b8+fOZMGFChXV79uxh5syZe9+HGtrvvvtuFi1atM8trvLycqpVq7bPo7iTn3rK14yKN9xwA82bN9+n42Pbtm1RVb7y8TRYNMuXL6dNmzYplU1FzCsTETkCmAB0AVREvgCuUNXRIjIRN2TKUhG5F3g88vFeY2KZ8+UcAO7odkfAkZiEvnHHimOq7rHq27cv119/PZdeeinvvPMO/fr1o379+ixdupQnn3ySVq1a7X3aq0WLFvTo0YPHH3+c+vXr0zes0f+MM85g9+7dXHzxxVx66aV8/PHHjHzgARo3bJgwhqKiIq677jruuOMOPv30U4488kg6d+5MjRo1WLBgwT5PeQH88MMPvPTSSxW2061bN/bbbz/A9YRPdhKudMS7MnkGWA4cADTGVSxTRKSOqq5X1StwT3WdC/wvy3EaY0zWjBo1ikmTJvHpp5/yy1/+kl69ejFq1Ch69OjBE088sU/e4uJiVJV+/fpRt+6PDyccc8wxjB8/nvfff58+ffrw/PPP8+KYMTSqH2/s2x9dffXVNGzYkJEj3YgD9erV42c/+xkzZlRst9q6dSvnn39+heXjjz8G4KOPPmL9+vV723NyQVSjj5MoIpuB81X1De99U9yc621U9bOIvOeq6tQsx5qyTp066fz581Mq2+rWaRmOJr8sf+CcjG+z+/jugA34mBdmd3evUQZ8XLJkCe3bt89pOHln2TL3muKAj1OnTuWyyy7j66+/Tqox/rbbbuODDz5I2M8k0TEUkQWq2snPZ8a7MnkdeFBEzhORs4G/Ap97yz4qc0VijDH5ql+/fhx44IFMnDjRd5nt27czbtw4hg8fnsXIKor3NNclwHDgVqAWbnTenhrrUsYYn5oVNQs6BONXbTtWaalePXGeOKpVq8a4ceNYFrrC8WHFihXceeeddI/WeTSL4vUz2YarSIzJqJcHRh/K21RCp9qxSssRR6S9iZNOOmnvoJB+tG/fPpDbj/HG5krcJTSD5YwxxuSveG0mn4jIZSLiqy+/iJwgIs9gVzMmgdtm38Zts28LOgzjx8Lb3BKD3fVOYNUqt1RCmT528dpMbgHuBv4oIq8D84BFuCe6vsc9LnwYbg72s4CWwF9wDfXGxPTuqneDDsH4tSH2sapZsybl5eUUFRXlMKA8s3170BHEVF5eTs2aNRNn9Clem8kUEZmKm9v9V7j53X/Cj3OuC/ADrmH+KWCiqq7LWGTGmEqtRYsWrF69moMOOoi6dev6GizRBE9VKS8vZ/Xq1XtHO86EuGNzeU9uveEtiMhPcJ0Y6wDfAsttkEdjClNDr2f3119/zc6dOwOOppIKjdq7Z0+wcUSoWbMm+++//95jmAlJzeeoqmuwuUuMMZ6GDRtm9AupyrnSm50jgyMdV1a5mRzYmDAHNzw46BCMX0V2rNJycOHsv5xXJiJyFPBn3ACSm3GN9nf7HShSRKrh5p4/Aeirqq9lKVSTJc/2fzboEIxfXe1YpeXZwtl/Oa1MRKQJMBtYDPQDDsdNvFUN19vej8uAwqnujTEmDyQ/40p6rgDqAv1V9Q1VfRL3+PH1IpLwxqtXGf0e+H/ZDdNk07CZwxg2c1jQYRg/Fgxzi0nNsGFuKQC+rkxE5BhVzcQw872BWaq6JSytBHgQ6Aa8mqD8vcA7wJwMxGICsnDtwqBDMH5tWhh0BPlt4cKgI8gZv1cm/xGRD0TkShFpnMbntQOWhieo6gqgjB/nlI9KRI7FDT55Yxqfb4wxJgv8ViZn4No5HgK+FpEXRKRXCuNwNcE1ukfa5K2L58/Ao5FzqcQiIkNFZL6IzF+/fn1yURpjjEmKr8pEVUtV9de4DotXAwcBs4CvROReETk8izEiIoOBtsB9fsuo6lhV7aSqnULTWBpjjMmOpBrgVXW7qv5VVU/DfbkvB27HDQo5V0TOTbCJTUCjKOlNvHUViEhN4GFcu0o17zZbqLG+nog0SOZ3MMFr06wNbZq1CToM40eDNm4xqWnTxi0FIOlHg0WkFTAEN15XS2A68HfgZ8AkEXlMVa+LUXwpEW0jItISKCKiLSVMPdyjwKO9JVwJbubH9CcNMDkztu/YoEMwfv3UjlVaxhbO/vP7NFcRMAC4GDgV+BIYB4z3hlgBeFpELgb+CMSqTGYAN4lIA1Xd6qUNAsqBuTHKbANOj0g7AHgBd1X0Tz+/gzHGmOzxe2XyDe6W2BTc1L2lMfJ9AGyMs50ngWuBKSLyINAaGAGMDn9cWEQ+A+aq6qWqugvY5/O8qyOA/6nq+z5/B1NJDH11KGBXKHnhfXes7AolRUO9/VcAVyh+K5ObgedV9bt4mVR1EW6Ok1jrN4lID+BRXJ+SzcAjuAolMq70Jk82ldYnGz8JOgTj11Y7Vmn5pHD2n9/KZD9c20WFysQblv5yVb3Hz4ZUdTHuUeN4eVolWL8cN5+KMcaYSsDv01x3EXs8rAO99cYYYwqU38pE+HGGxUgHE+OxXmOMMYUh5m0uEfk18GvvrQJPiMiWiGx1gGOA17MTnqmKOh7QMegQjF9NOgYdQX7r2DHoCHImXptJGT8+mSW49pJvI/L8gHvc9/HMh2aqqjFnjQk6BOPXCWOCjiC/jRkTdAQ5E7MyUdUXgRcBRORvwL2q+kWuAjPGGJM/fD3NpaoXZzsQUzgunHIhYDMu5oV57ljZjIsputDbfwUw42K8NpOHgD+p6irv53hUVW/JbGimqlq1ZVXQIRi/yuxYpWVV4ey/eFcm5wPPAau8n+NRwCoTY4wpUPHaTA6L9rMxxhgTKddzwBtjjKmC4rWZnJ3MhlR1evrhmELQ5eAuQYdg/GpuxyotXQpn/8VrM3kN1xbiZwwsxQZmND79oecfgg7B+NXRjlVa/lA4+y9eZWLtJMYYY3yJ1wD/VS4DMYXjvMnnAfDywJcDjsQk9JY7Vpxqxyol53n77+Wqv//itZkUqWpZ6OdEGwrlNSaRjWXx5k8zlcr3dqzSsrFw9l+821xbRaSLqv4bN3VurFGDQ6zNxBhjClS8yuQS4POwnxNVJsYYYwpUvDaTCWE/j89JNMYYY/KS32l7ARCRxkAH4CfAGmCRqm7OfFimKutxWI+gQzB+7W/HKi09Cmf/+apMRKQG8Hvgt0B4Y3yZiDwO/D9V3ZmF+EwVdEe3O4IOwfh1jB2rtNxROPvP73Aqo4HfAfcDRwHNvdc/ANcAo/x+oIgcJSJzRKRMRL4WkXtEJG7jvYgcLSIzvfzfi8gKEfmLiPzE7+caY4zJHr+3uS4CblfV0WFp3wK/F5EdwHDg2kQbEZEmwGxgMdAPOBxXEVXzthFLI+BL4Bnga1yHyruAE0TkRFXd5fP3MJVA7+d6AzDjghkBR2IS+pc7Vpxuxyolvb39N6Pq7z+/lcke4OMY6xbh/0mvK4C6QH9V3QK8ISINgREi8pCXVoGqzgPmhSWVisgq3NzzxwIf+vx8UwmU7ywPOgTj1247VmkpL5z95/c210TgshjrLgf8TiPWG5gVUWmU4CqYbj63ERLqDVQryXLGGGMyLF4P+KvC3i4HBojIx8A/gHVAC9ytqgbASJ+f1w74Z3iCqq4QkTJv3avxCotINS/mw4AHgA+Af/v8bGOMMVkS7zbXo1HSDgTaR0kfDfzRx+c1ATZHSd/krUtkOvAz7+cFwNmquidaRhEZCgwFOOSQQ3xs2hhjTKridVqsjBNnXQM0BY7ENdjPEJGTVXVHZEZVHQuMBejUqZP13q9E+rTpE3QIxq+D7FilpU/h7L+kOi1mwCbck1mRmnjr4lLVT70f3xeRt3BPeP0S+GvGIjRZd2PXG4MOwfjV3o5VWm4snP2XbA/4g4E2QJ3IdT5nWlyKaxsJ32ZLXEfIpcnEoqpfici3QOtkyhljjMk8vz3gGwCTgTNDSd5r+O0jP6MGzwBuEpEGqrrVSxsElANz/cQSFlNboBnu6sTkke7juwNQOqQ00DiMD7O7u9eepUFGkb+6d3evpaVBRpETfq9M/gAcApwKvA2ci7stdSFwBlDscztP4jo3ThGRB3FXFSOA0eGPC4vIZ8BcVb3Uez8S2AW8j2vAbw/cjBvVuMTnZxtjjMkSv43sZ+PG5nrfe/+1qr6pqkOBV4Cb/GxEVTcBPXBXMa8CdwOP4Hqzh6vBvlc683EV2dPANFyF9DJwkqpu9/k7GGOMyRK/Vyb7AytVdbeIbMc9URUyHffF7ouqLsZdzcTL0yrifQl2BWKMMZWW3yuTlbjBHQE+BcKfd/spUOHRXGOMMYXD75XJG0BPYCruttQEETkB+B44jSRGDTZm4NEDgw7B+HWIHau0DCyc/ee3MrkFbx4TVZ0oItuAAbgxta4GnspOeKYquurEqxJnMpVDGztWabmqcPafr8pEVcuAsrD3U3FXKcYkrWynO5WKahYlyGkCt8v7s69hxyolZd7+K6r6+y/ZTottgRNx0/Z+DcxX1WXZCMxUXWc/dzZg/UzyQqk7VtbPJEVne/vP+pk43pwj44DzcI3224D6wB4RmQJcFmsuEmOMMVWf36e5Hsf1fv8VUE9VGwL1gF8Dvbz1xhhjCpTf21z9gOtU9flQgqqWA8+JSBFuCHpjjDEFyu+VyTZgTYx1XwPWC90YYwqY3yuTx4AbReSf3hUJAN5VyY3YbS6ThCEdhwQdgvGr9ZCgI8hvQ4YEHUHOxJu296GIpCOBlSLyBj9O29sLN+Lv/KxFaKocq0zyiFUm6bHKBIDzI97v9JaTwtJCw8ifh8/BHo3ZULYBgOZFzRPkNIHb4Y4VdexYpWSDt/+aV/39F2/a3sNyGYgpHAMmDwCsn0leeNsdK+tnkqIB3v4rgH4mlXGed2OMMXnGd2UiIq1F5AkR+Z+IrPZeHxcRmzbXGGMKnN8e8CcA/8INNf8a8A1ujpPzgAtE5HRV/TBrURpjjKnU/D4aPBL4COjtDfoI7H00eLq3Pu6EV8YYY6ouv5VJZ2BgeEUCbjRhb372SRmPzFRZV3a6MugQjF9H2rFKy5WFs//8ViblQLMY65piMy2aJAzqMCjoEIxfh9qxSsugwtl/fhvgpwEPiMgp4Yne+z8Ar/r9QBE5SkTmiEiZiHwtIveISPUEZU4Ukb+JyGdeuWUicpeI1PH7uabyWPndSlZ+tzLoMIwf21e6xaRm5Uq3FAC/VybXA68Ac0VkHT/2gG8BvAvc4GcjItIEmA0sxg0eeThuyt9qwPA4RQd5eR/EzUF/LHCv93qez9/BVBIXTb0IsH4meeFdd6ysn0mKLvL2XwH0M/E70+JG4BQROYsfJ8daA7yvqq8n8XlX4Kb67e/Nf/KGN1fKCBF5KM6cKA+o6oaw96UisgN4SkQOVdWvkojBmLzR6tZpgX5+SeuNAAwOMI7lD5wT2Gcb/xJWJt6tpP8C16rqTGBmGp/XG5gVUWmU4K44uhHjdllERRLykfd6IGCViTHGBChhm4mq7gAaA3sy8HntgKUR21+Bm1++XZLb6uLF9HkG4jLGGJMGvw3wzwEXZ+DzmgCbo6Rv8tb5IiIH4NpYJqrquhh5horIfBGZv379+lRiNcYY45PfBvgVwEAR+QCYgesBr2HrVVWfyHRw0YhILWAybsKu62LlU9WxwFiATp06aax8Jvdu6OLreQ1TCYxbf27QIeS3GwrnXPdbmYzyXn8CnBBlvQJ+KpNNQKMo6U28dXGJiADPAEcDJ6tqwjKm8unbtm/QIRif5mz9adAh5Le+hXOu+32aK1OjCy8lom1ERFoCRUS0pcQwBvdIcS9V9ZPfVELLNiwDoG3ztgFHYhJpXXsVAF98f3DAkeSpZe5cp23VP9f9XplkygzgJhFpoKqhibUG4XrYz41XUERuA67GDevydnbDNNn0m9d+A1g/k3xw/0GPAjD4iwcCjiRP/cad69bPJIzXVjEEN07X3n4mwARV/cHnZp4ErgWmiMiDQGtgBDA6/HFhEfkMmKuql3rvfwncD4wHVotI+GyPn6uqtbAbY0yAfN2+EpH2uJ7njwEdgN3e62PAZyJylJ/teG0cPYDquD4ldwOPAHdFZK3h5Qk503sdgutxH75YjyZjjAmY3yuTscB3wKlevxAAROQQ3PwmTwKn+dmQqi4mwXD1qtoq4v0QXEVijDGmEvLbsN4JuDO8IoG9HQ7vwg2xYowxpkD5vTJZDsQaobcOrh+KMb4MPy3emJ6mMvnzusFBh5DfhhfOue63MrkVGCUiX6rq+6FEryH8XuDGbARnqqaerXsGHYLx6Z1tHYMOIb/1LJxz3W9lMhxoCMyLMgT9RuB2Ebk9lFlVO2c6UFN1LFy7EICOB3QMNA6T2FF1vgBg8Y7WAUeSpxYudK8dOwYZRU74rUwWeYsxaRs2cxhg/UzywZ0HjgWsn0nKhg1zr9bPxFHVTAzyaIwxporK1DApxhhjCphVJsYYY9JmlYkxxpi05XqgR2O4v8f9QYdgfHpo7a+DDiG/3V8457pVJibnurbsGnQIxqcPy9oHHUJ+61o457rd5jI5N2/lPOatnBd0GMaH44uWcHzRkqDDyF/z5rmlANiVicm52+e4/q3Wz6Tyu/mACYD1M0nZ7V5f7gLoZ2JXJsYYY9JmlYkxxpi0WWVijDEmbVaZGGOMSZs1wJucG3PWmKBDMD7d8/XQoEPIb2PGBB1BzlhlYnLOhp7PHzb0fJoKYOj5kJzf5hKRo0RkjoiUicjXInKPiFRPUKaWiDwsIm+JSLmIaK7iNZk3+4vZzP5idtBhGB9Orr+Qk+svDDqM/DV7tlsKQE6vTESkCTAbWAz0Aw4HRuEqtXjzWxYBlwH/BuYBZ2Q3UpNN9715H2AzLuaDa1qUADbjYsruc+d6Icy4mOvbXFcAdYH+qroFeENEGgIjROQhL60CVd0sIk1VVUXkaqwyMcaYSiXXt7l6A7MiKo0SXAXTLV5BVbVbW8YYU0nlujJpBywNT1DVFUCZt84YY0weyvVtribA5ijpm7x1GSMiQ4GhAIccckgmN22S0OrWaRXS1tbaGHNdVbP8gXOCDsGYnKiyjwar6lhgLECnTp3sFlkl0mzn1UGHYHy6fbUdq7Q89VTQEeRMriuTTUCjKOlNvHWmANTUg4MOwfj0xfd2rNLStm3QEeRMrttMlhLRNiIiLXGP/i6NWsJUOWXV3qes2vtBh2F86NHgfXo0sGOVsldfdUsByPWVyQzgJhFpoKpbvbRBQDkwN8exmIBsqTEVgKIffhpwJCaRy/dzx2rOVjtWKRk1yr327RtsHDmQ6yuTJ4HvgSki0tNrJB8BjA5/XFhEPhORp8MLikhvERkAdPTeD/CWQ3MWvTHGmKhyemWiqptEpAfwKPAq7smuR3AVSmRckUOsPAGEVxwveq8XA+MzHKoxxpgk5PxpLlVdTIIe7Krayk+aMcaYysHmMzHGGJO2KtvPxFRezXfeEHQIxqfrVtqxSsvEiUFHkDNWmZicq6H7BR2C8WnNTjtWaWnZMugIcsZuc5mc2179TbZXfzPoMIwPfRq9SZ9GdqxSNmmSWwqAXZmYnNtafToA9XafFnAkJpELm7lj9dp3dqxS8sQT7nXQoGDjyAG7MjHGGJM2q0yMMcakzSoTY4wxabPKxBhjTNqsAd7k3H4/3BZ0CManK7/K/2MV5CRsTY4dCsCmAGPI1QRtVpmYnKsedUobUxlt2m3HKh2bigpn/9ltLpNz26rPZlv12UGHYXwY0GQ2A5rYsUrVgP/NZsD/CmP/WWVics4qk/xhlUl6rDIxxhhjkmCViTHGmLRZZWKMMSZtVpkYY4xJmz0abHKuxQ8jgg7B+DTkyxFBh5DXhpw/IugQcsYqE5Nz1agTdAjGpx1qxyodO2oWzv6z21wm57ZWn8bW6sH1CDb+XdhsGhc2s2OVqgs/nMaFHxbG/st5ZSIiR4nIHBEpE5GvReQeEanuo1wjEfmbiGwSke9E5DkRaZaLmE1mba/+FturvxV0GMaHPo3eok8jO1ap6rP0LfosLYz9l9PbXCLSBJgNLAb6AYcDo3CV2vAExScDbYDLgD3Ag8DfgVOzFK4xxhifct1mcgVQF+ivqluAN0SkITBCRB7y0ioQkS7AmUA3VX3TS1sNvC8iPVW1MLqYGmNMJZXr21y9gVkRlUYJroLplqDcN6GKBEBV/w186a0zxhgToFxXJu2ApeEJqroCKPPW+S7nWZKgnDHGmBwQVc3dh4nsBG5S1TER6auAZ1T19hjl3gC2q+ovItKfBVqratcoZYYCQ723bYFlMcJqDmxI4tfINYsvPRZf+ip7jBZfeuLFd6iq7udnI1W2n4mqjgXGJsonIvNVtVMOQkqJxZceiy99lT1Giy89mYov17e5NkHUmZGaeOsyXc4YY0wO5LoyWUpEG4eItASKiN4mErOcJ1ZbijHGmBzKdWUyA/iZiDQISxsElANzE5Q7QEROCSWISCegtbcuHQlvhQXM4kuPxZe+yh6jxZeejMSX6wb4JrgOi4twnQ5bA6OBMao6PCzfZ8BcVb00LG0WcCRwIz92WlynqtZp0RhjApbTKxNV3QT0AKoDrwJ3A48Ad0VkreHlCTcId/XyV+AZYAFwbjbjNcYY409Or0yMMcZUTQU5arCIXC4in4rIDhFZICI9fJQZISIaZTkrxRgq/YCXqcQoIq1i7KeSDMd2hIg8JSL/FZHdIlLqs1xO9l8q8eVq33mfdb6I/ENEVovINu/voNhHudoiMkpE1onIdhGZJiKtKlF80fbfe1mIb4CIzBORjd73yDIRGS4itRKUy9X5l3R86Z5/VbafSSzeCfkkMAJ4G7gYeE1ETlTVRQmKfwdEVh5LUoih0g94mWaM4Nq23gl7n+lOW0cDZwPvATWTKJerAUNTjQ+yv+8ArscNR3Sdt/2zgedFpLmq/jlOuT8BA7xy63F/R2+IyDGquqMSxAfuPH0p7P3WDMYV0gz4J/AwsBnojNsXBwBXxymXq/Mv1fgg1fNPVQtqwfWE/2vY+2rA/4BnE5QbAWzIUAy34frHNAxLuxk3rEzDOOW6AAqcFpbW2UvrmeH9lGqMrbx4+mT5OFYL+/kloNRHmVzuv1Tiy8m+8z6reZS054Ev45Q5GNgF/Cos7SDgB+CyoOPz8ihwdbb3X4zP/j3ui1tirM/Z+ZdifGmdfwV1m0tEWuP+K5gcSlPVPcCL5HbAyHwY8DLVGHPCO27Jytn+SzG+nFHVaP9tfgQcGKfYmd7rlLDtrMZd4Wd6/6USX9A2AvFucwU9YG2i+NJSUJUJP3Z8jOzouARoKiKJxqBpLCIbRGSniHwkIv3TiKOyD3iZaowhf/PaCtaIyGgRqZvh+FKRLwOGBrXvugCfxFnfDlilqtsi0nO1/xLFFzJCRHZ5f6t/FZGm2QpIRKqLSJG4PnDXAk+o929+FDk//5KMLySl86/Q2kyaeK+bI9I3ha1fH6PsZ7jbPB8BDYDfAC+LyHmqOiVGmXhxRMYQiqNJlHQ/5VonGUMiqcb4PfAY8DqwBegO3IJrc+mX0QiTl8v9l4rA9p24h1B+AVwSJ1uq50TafMYHMAHX7WA90Am4AzhORDqr6u4shLYdqO39/AxwU5y8QZx/ycSX1vmX95WJiDQCfpIon6qmNeyKqj4b8bmvAvOAOwm77C90qrqGfRv4SkXkG+BxETlOVf8TUGiVXlD7znsa63ngFVUdn43PSEcy8anqkLC3b4rIEmA60BfX0J1pXXHDQXXGfRc8ClyVhc9Jle/40j3/qsJtrvNxl4mJFvjxCiRy0MgmEesT8i4VpwDHio9HeiPkw4CXmfys0JM1J6QVUfryccDQrO477xbQDOAr4IIE2XO+/5KML5qZwDbg+EzGFaKqH6rq26o6Gncb6UoROTxG9pzvvyTji8b3+Zf3lYmq/kVVJdHiZQ9dnUTen2wHfKuqsW5xxfx4b0lWPgx4mWqM0WjEa1DyccDQrO07ESkCXsM1yvZR1bIERZYCLUWkXkR6VvZfCvFVENY+kItz70Pv9bAY64M+/xLFF43v/Zf3lUkyVPULXAPe+aE0EanmvU9qwEgREeA84D8p3IutjANeZirGaAZ4rwsyEVgacrn/MiUr+05EauCeYjwSOEtV1/ko9rr3uncYIxE5ENdHIqP7L8X4om3nLKA+uTn3TvZev4yxPujzL1F80fg//7L9bHNlW4BiYDeu493pwHjcF2SHsDzdcM/TdwtLm4u7TDwT98c0Hdfp6OcpxNAEWAO8AfTEzQi5DbgvIt9nwNMRabOAL4D+uAbJZcBbWdhPKcWI648zyouvJ3CPt39fznB8Rd6JPgB4F/g47H1RJdh/SceXq33nfdZY3H+b1wInRSy1vTxzgDkR5Z7CdWK7CNeB9z3gU6BO0PF55+hYYCBwBq7z3WbgfaB6huOb6W2/N+474W7v76Mk1t9Gjs+/pONL9/zL6C+QLwtwubcjv8dd+vWIWN/dO5G7h6U97Z0E5bgnJN4CeqcRw1G4HqrluC/teyNPeGA5MD4irTHwN++PZAuuYbJCB68M7aekYwQGA/NxowX84O3ne0JfABmMrRU/3maMXFoFvf9SiS9X+y7ssxPFV0pEZ0vck0GjcU9Lbcf9U3VYZYgPN4jsO7j+FDuBlbge+42yEN+9uNHPt3nn0ofANUDNWH8bOT7/ko4v3fPPBno0xhiTtoJqMzHGGJMdVpkYY4xJm1Umxhhj0maViTHGmLRZZWKMMSZtVpkYY4xJm1UmpuCIm4I5G7MXpkVExovI/BTKiYgsFJFfx1hfGiN9gDeda7JjyxlTgVUmxuS/gUBTXAe4ZEwBBNeb3Zi0WGViTP67FpioqjtDCSJygIhMEpGNQDcR2Soi80SkQyiPutkgn8H1jDYmLVaZGBOFiJwhIu+LyA4R+UZEHheR+mHra4rISBFZISLfi8jXIjJVRGp56xuLyF+89B1evnFZiPMI3JwVL0WsGocbY+5q3CB9g3HT60ZOYvUycLyIHJ3p2ExhyfvJsYzJNO+LdSZukMvzgJbAA7jRXc/yst2Gm1/jVtworAcAZwOh9ofRuC/564C13jZOy0K4PXBjZEVOXNQNeFBVXxCR36jqNGBaZGFVXSIim3AD+32chfhMgbDKxJiK7sBNxvRz9aYXEJFvgUki0kVV38XNXPe8qk4IKzc57OfOwGOqOiksbZ/ZOjPkBGCJd8sq3Bqgi4jUjlIm0n9x8RqTMrvNZUxFnYGpuu88NS/jpiUIzUWxEBgiIjeLyLHe/DbhFgI3ichVItImi7EegBsSPtK1uPkrVgMdROTWODPsbfC2Y0zKrDIxpqKfAN+EJ3gVy0bcU1MA9wGP4ebT/g+wUkR+F1bkatyc43cCy0TkUxEZnIVY6+CmUtiHqs7C3Za71lt/CfCxiJwfmddbXycLsZkCYpWJMRWtAVqEJ3h9MZoB3wKo6g5VvVNVWwFtgEnAGG9mP1R1s6peq6oHAMfhJmh6TkSOynCs3+LmyKhAVTep6vO4yava4dpMHoiStbG3HWNSZpWJMRW9D5wb0ZmvP66N8e3IzKr6KW5Wu+9xE4pFrv8vcBPu7y3aHODpWEaUOb0jb7t5bSof4irESK1w01kbkzJrgDeFqpaIDIiSPhd3C+sj4O8i8gRwMPAgMMtrfEdEpuIeuf0INxPlANzf05ve+reBqbjZ7hQ3u+d24N8J4moSI67pqloWJf0d4E4R2U9V14elfyQij3jx1RaRvsBvgdnhhUWkHq6CuyNBXMbEZZWJKVQNgBejpJ+uqqUi0hu4H9dLfAvwAnBzWL55wCB+vOJYDJynqqHhUN4FhuD+69+N+1LvraqrEsTVOkZch+GmWY1UirtFdRYwMSy9BNee0wZ3G+tp3KPOwyLKnwmU4eYmNyZlNm2vMXlORP4IHKGq58RYX6qq3WOsewHYrqqXZTFEUwDsysSY/Pcw8ImItFFV320fItIS6Accm7XITMGwBnhj8px36+wS3CPN0YyPkX4wcIWqfpaNuExhsdtcxhhj0mZXJsYYY9JmlYkxxpi0WWVijDEmbVaZGGOMSZtVJsYYY9L2/wG7RT2NjOgL5AAAAABJRU5ErkJggg==\n", "text/plain": [ "
┌───────┐┌────────┐ ┌───────────┐\n", " state_0: ┤0 ├┤0 ├──────┤0 ├\n", " │ ││ │ │ │\n", " state_1: ┤1 ├┤1 ├──────┤1 ├\n", " │ P(X) ││ │ │ │\n", " state_2: ┤2 ├┤2 ├──────┤2 ├\n", " │ ││ │ │ │\n", " state_3: ┤3 ├┤3 ├──────┤3 ├\n", " └───────┘│ adder │┌────┐│ adder_dg │\n", "objective: ─────────┤ ├┤2 ├┤ ├\n", " │ ││ ││ │\n", " sum_0: ─────────┤4 ├┤0 F ├┤4 ├\n", " │ ││ ││ │\n", " sum_1: ─────────┤5 ├┤1 ├┤5 ├\n", " │ │└────┘│ │\n", " carry: ─────────┤6 ├──────┤6 ├\n", " └────────┘ └───────────┘" ], "text/plain": [ " ┌───────┐┌────────┐ ┌───────────┐\n", " state_0: ┤0 ├┤0 ├──────┤0 ├\n", " │ ││ │ │ │\n", " state_1: ┤1 ├┤1 ├──────┤1 ├\n", " │ P(X) ││ │ │ │\n", " state_2: ┤2 ├┤2 ├──────┤2 ├\n", " │ ││ │ │ │\n", " state_3: ┤3 ├┤3 ├──────┤3 ├\n", " └───────┘│ adder │┌────┐│ adder_dg │\n", "objective: ─────────┤ ├┤2 ├┤ ├\n", " │ ││ ││ │\n", " sum_0: ─────────┤4 ├┤0 F ├┤4 ├\n", " │ ││ ││ │\n", " sum_1: ─────────┤5 ├┤1 ├┤5 ├\n", " │ │└────┘│ │\n", " carry: ─────────┤6 ├──────┤6 ├\n", " └────────┘ └───────────┘" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# define the registers for convenience and readability\n", "qr_state = QuantumRegister(u.num_qubits, \"state\")\n", "qr_sum = QuantumRegister(agg.num_sum_qubits, \"sum\")\n", "qr_carry = QuantumRegister(agg.num_carry_qubits, \"carry\")\n", "qr_obj = QuantumRegister(1, \"objective\")\n", "\n", "# define the circuit\n", "state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name=\"A\")\n", "\n", "# load the random variable\n", "state_preparation.append(u.to_gate(), qr_state)\n", "\n", "# aggregate\n", "state_preparation.append(agg.to_gate(), qr_state[:] + qr_sum[:] + qr_carry[:])\n", "\n", "# linear objective function\n", "state_preparation.append(objective.to_gate(), qr_sum[:] + qr_obj[:])\n", "\n", "# uncompute aggregation\n", "state_preparation.append(agg.to_gate().inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])\n", "\n", "# draw the circuit\n", "state_preparation.draw()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we use QAE to estimate the expected loss, we validate the quantum circuit representing the objective function by just simulating it directly and analyzing the probability of the objective qubit being in the $|1\\rangle$ state, i.e., the value QAE will eventually approximate." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "state_preparation_measure = state_preparation.measure_all(inplace=False)\n", "sampler = Sampler()\n", "job = sampler.run(state_preparation_measure)\n", "binary_probabilities = job.result().quasi_dists[0].binary_probabilities()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact Expected Loss: 0.6641\n", "Exact Operator Value: 0.3789\n", "Mapped Operator value: 0.5749\n" ] } ], "source": [ "# evaluate the result\n", "value = 0\n", "for i, prob in binary_probabilities.items():\n", " if prob > 1e-6 and i[-(len(qr_state) + 1) :][0] == \"1\":\n", " value += prob\n", "\n", "print(\"Exact Expected Loss: %.4f\" % expected_loss)\n", "print(\"Exact Operator Value: %.4f\" % value)\n", "print(\"Mapped Operator value: %.4f\" % objective.post_processing(value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we run QAE to estimate the expected loss with a quadratic speed-up over classical Monte Carlo simulation." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact value: \t0.6641\n", "Estimated value:\t0.6863\n", "Confidence interval: \t[0.6175, 0.7552]\n" ] } ], "source": [ "# set target precision and confidence level\n", "epsilon = 0.01\n", "alpha = 0.05\n", "\n", "problem = EstimationProblem(\n", " state_preparation=state_preparation,\n", " objective_qubits=[len(qr_state)],\n", " post_processing=objective.post_processing,\n", ")\n", "# construct amplitude estimation\n", "ae = IterativeAmplitudeEstimation(\n", " epsilon_target=epsilon, alpha=alpha, sampler=Sampler(run_options={\"shots\": 100, \"seed\": 75})\n", ")\n", "result = ae.estimate(problem)\n", "\n", "# print results\n", "conf_int = np.array(result.confidence_interval_processed)\n", "print(\"Exact value: \\t%.4f\" % expected_loss)\n", "print(\"Estimated value:\\t%.4f\" % result.estimation_processed)\n", "print(\"Confidence interval: \\t[%.4f, %.4f]\" % tuple(conf_int))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cumulative Distribution Function\n", "\n", "Instead of the expected loss (which could also be estimated efficiently using classical techniques) we now estimate the cumulative distribution function (CDF) of the loss.\n", "Classically, this either involves evaluating all the possible combinations of defaulting assets, or many classical samples in a Monte Carlo simulation. Algorithms based on QAE have the potential to significantly speed up this analysis in the future.\n", "\n", "To estimate the CDF, i.e., the probability $\\mathbb{P}[L \\leq x]$, we again apply $\\mathcal{S}$ to compute the total loss, and then apply a comparator that for a given value $x$ acts as\n", "\n", "$$ \\mathcal{C}: |L\\rangle_n|0> \\mapsto \n", "\\begin{cases} \n", "|L\\rangle_n|1> & \\text{if}\\quad L \\leq x \\\\\n", "|L\\rangle_n|0> & \\text{if}\\quad L > x.\n", "\\end{cases} $$\n", "\n", "The resulting quantum state can be written as\n", "\n", "$$ \\sum_{L = 0}^{x} \\sqrt{p_{L}}|L\\rangle_{n_s}|1\\rangle + \n", "\\sum_{L = x+1}^{2^{n_s}-1} \\sqrt{p_{L}}|L\\rangle_{n_s}|0\\rangle, $$\n", "\n", "where we directly assume the summed up loss values and corresponding probabilities instead of presenting the details of the uncertainty model.\n", "\n", "The CDF($x$) equals the probability of measuring $|1\\rangle$ in the objective qubit and QAE can be directly used to estimate it." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
┌──────┐\n", "state_0: ┤0 ├\n", " │ │\n", "state_1: ┤1 ├\n", " │ cmp │\n", "compare: ┤2 ├\n", " │ │\n", " a15: ┤3 ├\n", " └──────┘" ], "text/plain": [ " ┌──────┐\n", "state_0: ┤0 ├\n", " │ │\n", "state_1: ┤1 ├\n", " │ cmp │\n", "compare: ┤2 ├\n", " │ │\n", " a15: ┤3 ├\n", " └──────┘" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# set x value to estimate the CDF\n", "x_eval = 2\n", "\n", "comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)\n", "comparator.draw()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "def get_cdf_circuit(x_eval):\n", " # define the registers for convenience and readability\n", " qr_state = QuantumRegister(u.num_qubits, \"state\")\n", " qr_sum = QuantumRegister(agg.num_sum_qubits, \"sum\")\n", " qr_carry = QuantumRegister(agg.num_carry_qubits, \"carry\")\n", " qr_obj = QuantumRegister(1, \"objective\")\n", " qr_compare = QuantumRegister(1, \"compare\")\n", "\n", " # define the circuit\n", " state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name=\"A\")\n", "\n", " # load the random variable\n", " state_preparation.append(u, qr_state)\n", "\n", " # aggregate\n", " state_preparation.append(agg, qr_state[:] + qr_sum[:] + qr_carry[:])\n", "\n", " # comparator objective function\n", " comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)\n", " state_preparation.append(comparator, qr_sum[:] + qr_obj[:] + qr_carry[:])\n", "\n", " # uncompute aggregation\n", " state_preparation.append(agg.inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])\n", "\n", " return state_preparation\n", "\n", "\n", "state_preparation = get_cdf_circuit(x_eval)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we first use quantum simulation to validate the quantum circuit." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
┌───────┐┌────────┐ ┌───────────┐\n", " state_0: ┤0 ├┤0 ├────────┤0 ├\n", " │ ││ │ │ │\n", " state_1: ┤1 ├┤1 ├────────┤1 ├\n", " │ P(X) ││ │ │ │\n", " state_2: ┤2 ├┤2 ├────────┤2 ├\n", " │ ││ │ │ │\n", " state_3: ┤3 ├┤3 ├────────┤3 ├\n", " └───────┘│ adder │┌──────┐│ adder_dg │\n", "objective: ─────────┤ ├┤2 ├┤ ├\n", " │ ││ ││ │\n", " sum_0: ─────────┤4 ├┤0 ├┤4 ├\n", " │ ││ cmp ││ │\n", " sum_1: ─────────┤5 ├┤1 ├┤5 ├\n", " │ ││ ││ │\n", " carry: ─────────┤6 ├┤3 ├┤6 ├\n", " └────────┘└──────┘└───────────┘" ], "text/plain": [ " ┌───────┐┌────────┐ ┌───────────┐\n", " state_0: ┤0 ├┤0 ├────────┤0 ├\n", " │ ││ │ │ │\n", " state_1: ┤1 ├┤1 ├────────┤1 ├\n", " │ P(X) ││ │ │ │\n", " state_2: ┤2 ├┤2 ├────────┤2 ├\n", " │ ││ │ │ │\n", " state_3: ┤3 ├┤3 ├────────┤3 ├\n", " └───────┘│ adder │┌──────┐│ adder_dg │\n", "objective: ─────────┤ ├┤2 ├┤ ├\n", " │ ││ ││ │\n", " sum_0: ─────────┤4 ├┤0 ├┤4 ├\n", " │ ││ cmp ││ │\n", " sum_1: ─────────┤5 ├┤1 ├┤5 ├\n", " │ ││ ││ │\n", " carry: ─────────┤6 ├┤3 ├┤6 ├\n", " └────────┘└──────┘└───────────┘" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "state_preparation.draw()" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "state_preparation_measure = state_preparation.measure_all(inplace=False)\n", "sampler = Sampler()\n", "job = sampler.run(state_preparation_measure)\n", "binary_probabilities = job.result().quasi_dists[0].binary_probabilities()" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Operator CDF(2) = 0.9551\n", "Exact CDF(2) = 0.9521\n" ] } ], "source": [ "# evaluate the result\n", "var_prob = 0\n", "for i, prob in binary_probabilities.items():\n", " if prob > 1e-6 and i[-(len(qr_state) + 1) :][0] == \"1\":\n", " var_prob += prob\n", "\n", "print(\"Operator CDF(%s)\" % x_eval + \" = %.4f\" % var_prob)\n", "print(\"Exact CDF(%s)\" % x_eval + \" = %.4f\" % cdf[x_eval])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we run QAE to estimate the CDF for a given $x$." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact value: \t0.9521\n", "Estimated value:\t0.9590\n", "Confidence interval: \t[0.9577, 0.9603]\n" ] } ], "source": [ "# set target precision and confidence level\n", "epsilon = 0.01\n", "alpha = 0.05\n", "\n", "problem = EstimationProblem(state_preparation=state_preparation, objective_qubits=[len(qr_state)])\n", "# construct amplitude estimation\n", "ae_cdf = IterativeAmplitudeEstimation(\n", " epsilon_target=epsilon, alpha=alpha, sampler=Sampler(run_options={\"shots\": 100, \"seed\": 75})\n", ")\n", "result_cdf = ae_cdf.estimate(problem)\n", "\n", "# print results\n", "conf_int = np.array(result_cdf.confidence_interval)\n", "print(\"Exact value: \\t%.4f\" % cdf[x_eval])\n", "print(\"Estimated value:\\t%.4f\" % result_cdf.estimation)\n", "print(\"Confidence interval: \\t[%.4f, %.4f]\" % tuple(conf_int))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Value at Risk\n", "\n", "In the following we use a bisection search and QAE to efficiently evaluate the CDF to estimate the value at risk." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "def run_ae_for_cdf(x_eval, epsilon=0.01, alpha=0.05):\n", "\n", " # construct amplitude estimation\n", " state_preparation = get_cdf_circuit(x_eval)\n", " problem = EstimationProblem(\n", " state_preparation=state_preparation, objective_qubits=[len(qr_state)]\n", " )\n", " ae_var = IterativeAmplitudeEstimation(\n", " epsilon_target=epsilon, alpha=alpha, sampler=Sampler(run_options={\"shots\": 100, \"seed\": 75})\n", " )\n", " result_var = ae_var.estimate(problem)\n", "\n", " return result_var.estimation" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "def bisection_search(\n", " objective, target_value, low_level, high_level, low_value=None, high_value=None\n", "):\n", " \"\"\"\n", " Determines the smallest level such that the objective value is still larger than the target\n", " :param objective: objective function\n", " :param target: target value\n", " :param low_level: lowest level to be considered\n", " :param high_level: highest level to be considered\n", " :param low_value: value of lowest level (will be evaluated if set to None)\n", " :param high_value: value of highest level (will be evaluated if set to None)\n", " :return: dictionary with level, value, num_eval\n", " \"\"\"\n", "\n", " # check whether low and high values are given and evaluated them otherwise\n", " print(\"--------------------------------------------------------------------\")\n", " print(\"start bisection search for target value %.3f\" % target_value)\n", " print(\"--------------------------------------------------------------------\")\n", " num_eval = 0\n", " if low_value is None:\n", " low_value = objective(low_level)\n", " num_eval += 1\n", " if high_value is None:\n", " high_value = objective(high_level)\n", " num_eval += 1\n", "\n", " # check if low_value already satisfies the condition\n", " if low_value > target_value:\n", " return {\n", " \"level\": low_level,\n", " \"value\": low_value,\n", " \"num_eval\": num_eval,\n", " \"comment\": \"returned low value\",\n", " }\n", " elif low_value == target_value:\n", " return {\"level\": low_level, \"value\": low_value, \"num_eval\": num_eval, \"comment\": \"success\"}\n", "\n", " # check if high_value is above target\n", " if high_value < target_value:\n", " return {\n", " \"level\": high_level,\n", " \"value\": high_value,\n", " \"num_eval\": num_eval,\n", " \"comment\": \"returned low value\",\n", " }\n", " elif high_value == target_value:\n", " return {\n", " \"level\": high_level,\n", " \"value\": high_value,\n", " \"num_eval\": num_eval,\n", " \"comment\": \"success\",\n", " }\n", "\n", " # perform bisection search until\n", " print(\"low_level low_value level value high_level high_value\")\n", " print(\"--------------------------------------------------------------------\")\n", " while high_level - low_level > 1:\n", "\n", " level = int(np.round((high_level + low_level) / 2.0))\n", " num_eval += 1\n", " value = objective(level)\n", "\n", " print(\n", " \"%2d %.3f %2d %.3f %2d %.3f\"\n", " % (low_level, low_value, level, value, high_level, high_value)\n", " )\n", "\n", " if value >= target_value:\n", " high_level = level\n", " high_value = value\n", " else:\n", " low_level = level\n", " low_value = value\n", "\n", " # return high value after bisection search\n", " print(\"--------------------------------------------------------------------\")\n", " print(\"finished bisection search\")\n", " print(\"--------------------------------------------------------------------\")\n", " return {\"level\": high_level, \"value\": high_value, \"num_eval\": num_eval, \"comment\": \"success\"}" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--------------------------------------------------------------------\n", "start bisection search for target value 0.950\n", "--------------------------------------------------------------------\n", "low_level low_value level value high_level high_value\n", "--------------------------------------------------------------------\n", "-1 0.000 1 0.752 3 1.000\n", " 1 0.752 2 0.959 3 1.000\n", "--------------------------------------------------------------------\n", "finished bisection search\n", "--------------------------------------------------------------------\n" ] } ], "source": [ "# run bisection search to determine VaR\n", "objective = lambda x: run_ae_for_cdf(x)\n", "bisection_result = bisection_search(\n", " objective, 1 - alpha, min(losses) - 1, max(losses), low_value=0, high_value=1\n", ")\n", "var = bisection_result[\"level\"]" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated Value at Risk: 2\n", "Exact Value at Risk: 2\n", "Estimated Probability: 0.959\n", "Exact Probability: 0.952\n" ] } ], "source": [ "print(\"Estimated Value at Risk: %2d\" % var)\n", "print(\"Exact Value at Risk: %2d\" % exact_var)\n", "print(\"Estimated Probability: %.3f\" % bisection_result[\"value\"])\n", "print(\"Exact Probability: %.3f\" % cdf[exact_var])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Conditional Value at Risk\n", "\n", "Last, we compute the CVaR, i.e. the expected value of the loss conditional to it being larger than or equal to the VaR.\n", "To do so, we evaluate a piecewise linear objective function $f(L)$, dependent on the total loss $L$, that is given by\n", "\n", "$$\n", "f(L) = \\begin{cases} \n", "0 & \\text{if}\\quad L \\leq VaR \\\\\n", "L & \\text{if}\\quad L > VaR.\n", "\\end{cases}\n", "$$\n", "\n", "To normalize, we have to divide the resulting expected value by the VaR-probability, i.e. $\\mathbb{P}[L \\leq VaR]$." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
┌────┐\n", "q158_0: ┤0 ├\n", " │ │\n", "q158_1: ┤1 ├\n", " │ │\n", " q159: ┤2 F ├\n", " │ │\n", " a83_0: ┤3 ├\n", " │ │\n", " a83_1: ┤4 ├\n", " └────┘" ], "text/plain": [ " ┌────┐\n", "q158_0: ┤0 ├\n", " │ │\n", "q158_1: ┤1 ├\n", " │ │\n", " q159: ┤2 F ├\n", " │ │\n", " a83_0: ┤3 ├\n", " │ │\n", " a83_1: ┤4 ├\n", " └────┘" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# define linear objective\n", "breakpoints = [0, var]\n", "slopes = [0, 1]\n", "offsets = [0, 0] # subtract VaR and add it later to the estimate\n", "f_min = 0\n", "f_max = 3 - var\n", "c_approx = 0.25\n", "\n", "cvar_objective = LinearAmplitudeFunction(\n", " agg.num_sum_qubits,\n", " slopes,\n", " offsets,\n", " domain=(0, 2**agg.num_sum_qubits - 1),\n", " image=(f_min, f_max),\n", " rescaling_factor=c_approx,\n", " breakpoints=breakpoints,\n", ")\n", "\n", "cvar_objective.draw()" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "
Software | Version |
---|---|
qiskit | None |
qiskit-terra | 0.45.0.dev0+c626be7 |
qiskit_ibm_provider | 0.6.1 |
qiskit_aer | 0.12.0 |
qiskit_algorithms | 0.2.0 |
qiskit_finance | 0.4.0 |
System information | |
Python version | 3.9.7 |
Python compiler | GCC 7.5.0 |
Python build | default, Sep 16 2021 13:09:58 |
OS | Linux |
CPUs | 2 |
Memory (Gb) | 5.778430938720703 |
Fri Aug 18 16:24:38 2023 EDT |
© Copyright IBM 2017, 2023.
This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.