Credit Risk Analysis¶

This page demonstrates the 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, while allowing for flexible design of each block—using built-in functions or user-defined functions. The gate-level details are assigned to the synthesis engine.

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 \(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 focuses 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 computing with a quadratic speedup, by employing quantum amplitude estimation [1], [2].

This example highlights 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 on the Classiq platform.

Practical schemes for VaR involve a (classical) search over the CDF, such as a 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 the five main blocks listed below. The uncertainty model, encoding the joint probability of the indicators \(X_{k}\) (Eq.1), is Gaussian Conditional Independence (GCI). The Appendix contains a detailed explanation of the model.

The blocks:

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

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

Credit_Risk_example

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

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

Step 1 - Creating the Circuit¶

The example begins by designing the model, using \(14\) qubits and a maximal depth constraint of \(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 2 - Preparing the Gaussian State¶

Use Classiq state preparation to prepare a standard normal distribution state, which creates a superposition of \(2^{4}\) equally spaced states in \(\left [-5,5\right)\). See the state preparation documentation.

Specify the state preparation output and feed it 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 3 - Gaussian Conditional Independence Model¶

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

Use LinearPauliRotations to employ this step. Calculate the slope and the offset using 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 being the Gaussian state from the previous step (see Appendix).

For demonstration purposes, the Classiq platform includes the LinearGCI function, 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, taking both parameters 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. Specify the output, the indicators \(X_{k}\), to feed 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 4 - Calculating the Loss¶

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

For input, connect the indicators from the previous block.

For the outputs, the loss realization is fed into the CDF encoder block, and the indicators are also specified in the outputs since they are required for the loss uncomputation 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 5 - Encoding the Cumulative Distribution Function¶

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

This step employs 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 to 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 fed to another block and therefore is not required for the circuit synthesis process. (It is 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 6 - Determining the Loss by Uncomputation¶

The last step is to uncompute the loss register so it can be used again in the next oracle call. Call the WeightedAdder function with exactly the same parameters, but using the inverse flag and setting its value 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¶

Connect the building blocks 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 you explicitly specifying them.

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()

Textual Model

{
  "function_library": {
    "functions": [
      {
        "name": "main",
        "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"
            }
          }
        ]
      }
    ]
  },
  "constraints": {
    "max_width": 14,
    "max_depth": 1000
  }
}

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 the previous section.
2. Perform amplitude estimation using the Classiq Executor. Here, choose \(\alpha=0 05\) and \(\epsilon=0.03\).
3. Collect the amplitude estimation result.
Plot the CDF.

The amplitude estimation requires specifying the objective qubit. 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 a calculation to account for the relations between different blocks, such as shared qubits. The Classiq platform gives you the objective qubit 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

These are the corresponding output values:

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 this 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, write this:

\[ 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 you 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 like this:

\[ \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}} \]

Define this:

\[ \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}} \]

The overall wave function after the previous step is this:

\[ \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]\). Obtain this probability using amplitude estimation. A classical search algorithm iteratively modifies \(x\) to eventually obtain the VaR (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}} \]

References¶

[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.