Credit Risk Analysis¶

Introduction¶

This guide demonstrates a functional level construction of a quantum algorithm for Value at Risk, described in [1], [2], using the Classiq platform. The functional level description puts the design focus on assembling the building blocks of the algorithm and generating the logical connections between these blocks, all while allowing for flexible design of each block - either by using built-in functions (see built-in functions) or user defined functions. The gate-level details are left to the synthesis engine.

Credit Risk Analysis Overview¶

Given a portfolio of \(K\) assets, the overall loss is given by summing the individual losses, \(L=\sum_{k=1}^{K}\lambda_{k}X_{k}\), where \(\lambda_{k}\) is the loss value of each asset \(k\) and \(X_{k}\) is an indicator for a loss (default), getting values of either \(0\) or \(1\).

For any non-trivial joint probability distribution of the indicators \(X_{k}\), classical algorithms rely on computationally intensive Monte Carlo simulations to obtain quantities of interest related to statistics of the overall loss.

This example will focus on Value at Risk (VaR). The VaR is defined as the maximal loss up to a given confidence,

\[ \rm{VaR}_{\alpha}\left(L\right) = \inf \left\{ x \mid \mathbb{P} \left[L<= x\right] \geq \alpha \right\}, \tag{1} \]

with \(\alpha \in \left[0, 1\right]\) being the confidence level. This computation can potentially benefit from quantum computer by a quadratic speedup, by employing quantum amplitude estimation [1], [2].

This example will highlight the benefits of designing the amplitude estimation operator \(\mathcal{A}\), which encodes the cumulative distribution function (CDF), at a functional level, and then executing the amplitude estimation process, in the Classiq platform.

Practical schemes for Value at Risk involve a (classical) search over the CDF, such as bisection search. Here, for demonstration purposes and to keep the example simple, the CDF is calculated for all possible values.

Building the Amplitude Estimation Operator¶

The construction is divided into 5 main blocks, as listed below. The uncertainty model, encoding the joint probability of the indicators \(X_{k}\) (Eq.1), is Gaussian Conditional Independence (GCI), and the Appendix contains the detailed explanation of the model.

The blocks are:

Gaussian state preparation (state preparation)
Gaussian Conditional Independence probability calculation (linear Pauli rotations)
Loss calculation (weighted adder)
Cumulative distribution function encoding (comparator)
Loss un-computation (inverse of the weighted adder)

The parentheses denote the function used for each block. The blocks are connected according to the following schematics:

Credit_Risk_example

For demonstration purposes, each block is a Classiq's built-in function. However, user-created functions can be used in each block instead.

A step-by-step Python example will be shown below. The textual model implementation is shown afterwards.

Step 0 - Circuit Creation¶

The example begins by designing the model. Here, \(14\) qubits are used and the maximal depth constraint is \(1000\).

from classiq import Model
from classiq.model import Constraints

cdf_encoder_max_width = 14
cdf_encoder_max_depth = 1100
constraints = Constraints(
    max_width=cdf_encoder_max_width, max_depth=cdf_encoder_max_depth
)


def create_circuit():
    model = Model(constraints=constraints)
    return model

Step 1 - Gaussian State Preparation¶

We prepare a standard normal distribution state by employing Classiq's state preparation, which creates a superposition of \(2^{4}\) equally spaced states in \(\left [-5,5\right)\). Refer to the state preparation documentation for more details.

We specify the state preparation output, and it will be fed into the next block.

from classiq import QReg
from classiq.builtin_functions import StatePreparation

num_sp_qubits = 4
probabilities = {
    "gaussian_moment_list": [{"mu": 0, "sigma": 1}],
    "num_qubits": num_sp_qubits,
}
error_metric = {"L2": {"upper_bound": 0.0005}}

sp_output_gauss_dist = QReg(size=num_sp_qubits)


def add_state_preparation(model):
    state_preparation_params = StatePreparation(
        probabilities=probabilities, error_metric=error_metric
    )
    model.StatePreparation(
        state_preparation_params, out_wires={"OUT": sp_output_gauss_dist}
    )

Step 2 - The Gaussian Conditional Independence Model¶

This step encodes the joint probability distribution of the indicators \(X_{k}\) (Eq.1).

To employ this step, we use LinearPauliRotations. The slope and the offset can be calculated by a Taylor approximation as a function of the GCI model parameters \(\rho_{k}\) and \(p^{0}_{X_{k}}\), with \(k\) denoting each asset, and the control state is the Gaussian state from the previous step (see Appendix for more details).

For demonstration purposes, Classiq platform includes a function LinearGCI which receives \(\rho_{k}\) and \(p^{0}_{X_{k}}\), performs the Taylor approximation and builds the corresponding LinearPauliRotations. LinearGCI also requires the number of qubits forming the Gaussian grid and the truncation value, both parameters are from the previous step.

The example uses 3 assets, with rhos and p_zeros corresponding to \(\rho_{k}\) and \(p^{0}_{X_{k}}\), respectively, and \(k=1,2,3\) denoting each asset.

The input is identical to the StatePreparation output, thus forming a wire that connects between them. We also specify the output, the indicators \(X_{k}\), which will be fed into the next block to calculate the loss.

from classiq.builtin_functions import LinearGCI

truncation_value = 5 * probabilities["gaussian_moment_list"][0]["sigma"]
rho_per_asset = [0.1, 0.4, 0.1]
p_zero_per_asset = [0.4, 0.2, 0.3]
num_assets = len(rho_per_asset)

gci_output_indicators = QReg(size=num_assets)


def add_gci(model):
    linear_gci_params = LinearGCI(
        num_state_qubits=num_sp_qubits,
        truncation_value=truncation_value,
        rhos=rho_per_asset,
        p_zeros=p_zero_per_asset,
    )
    model.LinearGCI(
        linear_gci_params,
        in_wires={"state": sp_output_gauss_dist},
        out_wires={"target": gci_output_indicators},
    )

Step 3 - Loss Calculation¶

the previous two steps calculated the probabilities of each possible realization of the \(X_{k}\)'s according to the Gaussian Conditional Independence model. This steps calculates the overall loss for each realization, by using the WeightedAdder function.

As an input, we connect the indicators from the previous block. Concerning the outputs: the loss realization is fed into the CDF encoder block; the indicators are also specified in the outputs, since they are required for the loss un-computation block.

from classiq.builtin_functions import WeightedAdder


num_assets = len(p_zero_per_asset)
loss_per_asset = [2, 1, 3]
num_loss_realization_qubits = sum(loss_per_asset).bit_length()  # max loss

lc_output_indicators = QReg(size=num_assets)
lc_output_loss_realization = QReg(size=num_loss_realization_qubits)


def add_loss_realization_calculation(model):
    weighted_adder_params = WeightedAdder(
        num_state_qubits=num_assets, weights=loss_per_asset
    )
    model.WeightedAdder(
        weighted_adder_params,
        in_wires={"state": gci_output_indicators},
        out_wires={"state": lc_output_indicators, "sum": lc_output_loss_realization},
    )

Step 4 - Cumulative Distribution Function Encoding¶

To calculate the VaR, it is basically required to obtain the cumulative distribution function (CDF) and search (e.g. binary) for the lowest possible loss value for a given confidence.

For this step, we employ a comparator, which performs a less than or equal comparison between the following terms: the loss realization register, which is fed from the previous block, and the loss constant at which we calculate the CDF, which is given as a parameter for the Python function.

The size of the loss realization register is the number of bits required to express the maximal possible loss.

Note that the comparison result bit, which encodes the CDF, is not specified, since it is not being fed to another block and therefore is not required for the circuit synthesis process. It will be referred to during the amplitude estimation step.

from classiq.builtin_functions import LessEqual


cdfe_output_loss_realization = QReg(size=num_loss_realization_qubits)


def add_cdf_encoder(model, loss):
    cdf_comparator_params = LessEqual(
        left_arg={"size": num_loss_realization_qubits},
        right_arg=loss,
        output_name="is_less_equal",
    )

    less_equal_outputs = model.LessEqual(
        cdf_comparator_params,
        in_wires={"left_arg": lc_output_loss_realization},
        out_wires={"left_arg": cdfe_output_loss_realization},
    )
    model.set_outputs({"objective": less_equal_outputs["is_less_equal"]})

Step 5 - Loss Un-computation¶

The last step is to un-compute the loss register, so it could be used again in the next oracle call. The WeightedAdder function is called again with the exact same parameters, however the inverse flag is now used and its value is set to True.

Note that the inputs for this block are the indicators, and the loss realization, the latter is passed after going through the comparator.

def add_loss_uncomputation(model):
    weighted_adder_inverse_params = WeightedAdder(
        num_state_qubits=num_assets, weights=loss_per_asset
    )
    model.WeightedAdder(
        weighted_adder_inverse_params,
        is_inverse=True,
        in_wires={"state": lc_output_indicators, "sum": cdfe_output_loss_realization},
    )

Full Circuit for Amplitude Estimation Operator¶

The building blocks are connected together, by specifying the relevant inputs and outputs for each block. The function receives the loss as a parameter and encodes the CDF for that value. The generated circuit image for \(\rm{loss} =3\) is presented above. Note how the auxiliary qubits of the two weighted adders and the comparator are reused, without explicit specification from the user.

SDK

def cdf_amplitude_estimation_encoder(loss):
    model = create_circuit()
    add_state_preparation(model)
    add_gci(model)
    add_loss_realization_calculation(model)
    add_cdf_encoder(model, loss)
    add_loss_uncomputation(model)
    circuit = model.synthesize()
    return circuit


cdf_amplitude_estimation_encoder(loss=3).show_interactive()

Textual Model

{
  "constraints": {
    "max_width": 14,
    "max_depth": 1000
  },
  "logic_flow": [
    {
      "function": "StatePreparation",
      "function_params": {
        "probabilities": {
          "gaussian_moment_list": [
            {
              "mu": 0.0,
              "sigma": 1.0
            }
          ],
          "num_qubits": 4
        },
        "error_metric": {
          "L2": {
            "upper_bound": 0.0005
          }
        }
      },
      "outputs": "out_sp"
    },
    {
      "function": "LinearGCI",
      "function_params": {
        "num_state_qubits": 4,
        "truncation_value": 5.0,
        "p_zeros": [0.4, 0.2, 0.3],
        "rhos": [0.1, 0.4, 0.1]
      },
      "inputs": {
        "state": "out_sp"
      },
      "outputs": {
        "target": "target_gci"
      }
    },
    {
      "function": "WeightedAdder",
      "function_params": {
        "num_state_qubits": 3,
        "weights": [2, 1, 3]
      },
      "inputs": {
        "state": "target_gci"
      },
      "outputs": {
        "sum": "sum_wa",
        "state": "state_wa"
      }
    },
    {
      "function": "LessEqual",
      "function_params": {
        "left_arg": {
          "size": 3
        },
        "right_arg": 3.0
      },
      "inputs": {
        "left_arg": "sum_wa"
      },
      "outputs": {
        "left_arg": "left_arg_c"
      }
    },
    {
      "function": "WeightedAdder",
      "function_params": {
        "num_state_qubits": 3,
        "weights": [2, 1, 3]
      },
      "is_inverse": true,
      "inputs": {
        "state": "state_wa",
        "sum": "left_arg_c"
      }
    }
  ]
}

Executing the Amplitude Estimation - Calculating the Cumulative Distribution Function¶

The final step involves executing the amplitude estimation, to calculate the cumulative distribution function (CDF), on the Classiq platform.

The example works as follows:

Loop over all possible loss values
1. Generate the operator \(\mathcal{A}\), the CDF encoder, for the specific loss, using the Classiq synthesis engine, similar to previous section.
2. Perform amplitude estimation using the Classiq Executor. Here, we chose \(\alpha=0. 05\) and \(\epsilon=0.03\).
3. Collect the amplitude estimation result
Plot the CDF

The amplitude estimation requires the objective qubit to be specified. When building circuits at the gate level, obtaining the objective qubit is a non-trivial task: it requires good knowledge of the specific implementations of each block, and calculation to account for the relations between different blocks, such as shared qubits. The Classiq platform allows to obtain it using the model outputs.

import numpy as np
from classiq import Executor
from classiq.execution import (
    ExecutionPreferences,
    AmplitudeEstimation,
    IBMBackendPreferences,
    QuantumProgram,
)


def run_ae(loss):
    sp_circuit_x = cdf_amplitude_estimation_encoder(loss)
    sp_quantum_program = QuantumProgram(code=sp_circuit_x.qasm)
    objective_qubits = sp_circuit_x.data.qubit_mapping.logical_outputs["objective"]
    preferences = ExecutionPreferences(
        amplitude_estimation=AmplitudeEstimation(
            alpha=0.05, epsilon=0.03, objective_qubits=objective_qubits
        ),
        backend_preferences=IBMBackendPreferences(backend_name="aer_simulator"),
    )

    res = Executor(preferences=preferences).execute(sp_quantum_program)
    cdf_res = res.vendor_format_result["estimation"]
    return cdf_res


max_loss = sum(loss_per_asset)
loss_vals = range(max_loss + 1)
cdf = np.zeros(len(loss_vals))

for ix, loss in enumerate(loss_vals):
    # reset registers
    sp_output_gauss_dist = QReg(size=num_sp_qubits)
    gci_output_indicators = QReg(size=num_assets)
    lc_output_indicators = QReg(size=num_assets)
    lc_output_loss_realization = QReg(size=num_loss_realization_qubits)
    cdfe_output_loss_realization = QReg(size=num_loss_realization_qubits)

    cdf_res = run_ae(loss)
    print(f"CDF[{loss}]={np.round(cdf_res, 5)}")
    cdf[ix] = cdf_res

The corresponding output values are

CDF[0]=0.37708
CDF[1]=0.4549
CDF[2]=0.65208
CDF[3]=0.8236
CDF[4]=0.86427
CDF[5]=0.95492
CDF[6]=0.99895

The resulting plot and the code to generate it are shown below.

Credit_Risk_example

from matplotlib import pyplot as plt

with plt.style.context("seaborn"):
    plt.plot(loss_vals, cdf, "o", markersize=10)
    plt.xlabel("Loss")
    plt.ylabel("CDF")
    plt.title("Loss Cumulative Distribution Function")
    plt.show()

Appendix - Value at Risk Algorithm and the Uncertainty Model¶

Each indicator of the loss \(X_{k}\) (see overview) is a Bernoulli variable with a probability \(p_{X_{k}}\), whose uncertainty is modeled for the analysis. In the simplest version of the model, the indicators are independent variables. More realistically, however, the losses are correlated, and the assumption is relaxed towards conditional independence only: The probabilities \(p_{X_{k}}\) depend on a latent normal random variable \(Z\), such that \(\left\{X_{k} | Z = z\right\}\) is a set of independent variables.

The explicit dependency \(p_{X_{k}}\left (z\right)\) is called Gaussian Conditional Independence and is given by the following expression:

\[ p_{X_{k}}\left(z\right) = F\left( \frac{F^{-1}(p_{X_{k}}^0) - \sqrt{\rho_k}z}{\sqrt{1 - \rho_k}} \right) \tag{1} \]

\(F\) refers to the Gaussian cumulative distribution function. \(\rho_ {k}\) refers to the sensitivity to changes in \(Z\), and \(p_{X_{k}}^{0}\) is the probability given zero sensitivity.

Therefore, we write:

\[ L\left(z\right)=\sum_{k=1}^{K}\lambda_{k}X_{k}\left(z\right) \tag{2} \]

Note that setting \(\rho_{k}=0 \forall k\) brings us back to the simplest model of independent losses.

To reflect the dependency on the latent random variable \(Z\), the first step will be to prepare a state with Gaussian distributed weights. Denote by \(n_{z}\) the number of bits representing the possible values of the random variable \(Z\), the resulting state is given by:

\[ \left|\psi_{1}\right\rangle=\sum_{i=0}^{2^{n_{z}}-1}g\left(z_ {i}\right)\left|z_{i}\right\rangle_{n_{Z}}, \]

where \(g\left(z\right)\) is the standard normal PDF (\(\mu=0\), \(\sigma=1\)).

The next step involves additional \(K\) qubits, representing the probability distribution of the indicators \(X_{k}\) in the Gaussian Conditional Independence (GCI) model:

\[ \left|\psi_{2}\right\rangle=\sum_{i=0}^{2^{n_{z}}-1}g\left(z_ {i}\right)\left|z_{i}\right\rangle_{n_{Z}} \bigotimes_{k=1}^{K}\left (\sqrt{p_ {X_{k}=0}\left(z_{i}\right)}\left|0\right\rangle +\sqrt{p_{X_ {k}=1}\left(z_{i}\right)}\left|1\right\rangle \right) \]

The probability distributions \(p_{X_{k}}\left(z\right)=1\) can be generated by Y rotations. Implicitly, \(\left|\psi_{1}\right\rangle\) includes \(K\) zero valued qubits, each is transformed as follows:

\[ \left|z_{i}\right\rangle _{n_{Z}}\left|0\right\rangle _{K}\rightarrow\left|z_{i}\right\rangle_{n_{Z}} \bigotimes_{k=1}^{K}\left( \cos\left(\alpha_{k}\left(z_{i}\right)\right)\left|0\right\rangle -\sin\left(\alpha_{k}\left(z_{i}\right)\right)\left|1\right\rangle \right) \]

Ideally, \(\alpha_{k}\left(z_{i}\right)=\sin^{-1}\left(\sqrt{p_{X_ {k}=1}\left(z_{i}\right)}\right)\). Here, we perform a linear approximation.

Now, we can interchange between the product and the summation over the \(\left|0\right\rangle\) and \(\left|1\right\rangle\) states by rewriting the expression as

\[ \left|\psi_{2}\right\rangle =\sum_{i=0}^{2^{n_{z}}-1}g\left(z_ {i}\right)\left|z_{i}\right\rangle \sum_{x_{1,}x_{2},\dots x_ {K}=\left\{ 0,1\right\} ^{k}}\left(\prod_{k=1}^{K}\sqrt{p_{X_{k}=x_ {k}}\left(z_{i}\right)}\right)\left|x_{1,}x_{2},\dots x_ {K}\right\rangle \]

Denote by \(n_{L}\) the number of qubits required to represent the maximal loss (the loss such that all indicators are \(1\)). These additional zero-input qubits form the register representing the loss:

\[ \left|x_{1},x_{2},\dots x_{K}\right\rangle \left|0\right\rangle_{n_{L}}\rightarrow\left|x_{1},x_{2},\dots x_ {K}\right\rangle \otimes\left|\lambda_{1}x_{1}+\dots+\lambda_{k}x_ {K}\right\rangle_{n_{L}} \]

Now, define

\[ \left|\lambda_{1}x_{1}+\dots+\lambda_{K}x_{K}\right\rangle _{n_{L}}\equiv\left|L\left(\left\{ x_{k}\right\} \right)\right\rangle _{n_{L}} \]

Then, the overall wave function after the previous step will be

\[ \left|\psi_{3}\right\rangle =\sum_{i=0}^{2^{n_{z}}-1}\sum_{x_{1}, \dots x_{K}=\left\{ 0,1\right\} ^{K}}\sqrt{g\left(z_{i}\right) \bigotimes_{k=1}^{K}\left(p_{X_{k}=x_{k}}\right)}\left|z_{i} \right\rangle _{n_{z}}\left|x_{1},x_{2},\dots x_{k}\right\rangle \left|L\left(\left\{ x_{k}\right\} \right)\right\rangle_{n_{L}} \]

The next step will be to apply a desired function of the loss \(L\), and encode it in an additional objective qubit to extract it via amplitude estimation. For the Value at Risk, the loss is fed into a comparator against a value \(x\). After the comparison, the resulting state is given by

\[ \left|\psi_{4}\right\rangle =\sqrt{1-\mathbb{P}\left[L\leq x\right]} \left|\psi_{3}\right\rangle _{L>x}\left|0\right\rangle +\sqrt{\mathbb{P}\left[L\leq x\right]}\left|\psi_{3}\right\rangle \_{L\leq x}\left|1\right\rangle + \]

where \(\left|\psi_{3}\right\rangle_{L>x}\), \(\left|\psi_{3} \right\rangle _{L\leq x}\) are the normalized wavefunctions corresponding to the terms in the summation in \(\left|\psi_{3}\right\rangle\) such that \(L>x\), \(L\leq x\), respectively.

The probability that the objective qubit is \(1\) is the CDF at \(x\): \(\mathbb{P}\left[L<= x\right]\). This probability will be obtained via an amplitude estimation. \(x\) will be iteratively modified via a classical search algorithm to eventually obtain the Value at Risk (Eq. 1).

The final step is to un-compute the loss bits. These are auxiliary qubits that have to be cleaned when the circuit is applied consecutively, as part of the amplitude estimation process:

\[ \left|\psi_{5}\right\rangle =\sqrt{1- \mathbb{P}\left[L\leq x\right]}\left|\tilde{\psi}_{3}\right\rangle _{L>x}\left|0\right\rangle +\sqrt{\mathbb{P}\left[L\leq x\right]} \left|\tilde{\psi}_{3}\right\rangle _{L\leq x}\left|1\right\rangle \]

where

\[ \left|\tilde{\psi}_{3}\right\rangle =\sum_{i=0}^{2^{n_{z}}-1} \sum_{x_{1},\dots x_{K}=\left\{ 0,1\right\} ^{K}}\sqrt{g\left(z_{i} \right)\bigotimes_{k=1}^{K}\left(p_{X_{k}=x_{k}}\right)}\left|z_{i} \right\rangle _{n_{z}}\left|x_{1},x_{2},\dots x_{k}\right\rangle \left|0\right\rangle _{n_{L}} \]

[1] Woerner, S., Egger, D.J. Quantum risk analysis. npj Quantum Inf 5, 15 (2019).

[2] D. Egger, R. Garcia Gutierrez, J. Mestre and S. Woerner, "Credit Risk Analysis Using Quantum Computers" in IEEE Transactions on Computers, vol. 70, no. 12, pp. 2136-2145, 2021.