Quantum Computation for Pricing the Collateralized Debt Obligations
Hao Tang, Anurag Pal, Lu-Feng Qiao, Tian-Yu Wang, Jun Gao, Xian-Min Jin
QQuantum Computation for Pricing the Collateral Debt Obligations
Hao Tang,
1, 2, ∗ Anurag Pal,
1, 2
Lu-Feng Qiao,
1, 2
Tian-Yu Wang,
1, 2
Jun Gao,
1, 2 and Xian-Min Jin
1, 2, † Center for Integrated Quantum Information Technologies (IQIT),School of Physics and Astronomy and State Key Laboratory of Advanced Optical Communication Systems and Networks,Shanghai Jiao Tong University, Shanghai 200240, China CAS Center for Excellence and Synergetic Innovation Center in Quantum Information and Quantum Physics,University of Science and Technology of China, Hefei, Anhui 230026, China
Collateral debt obligation (CDO) has been oneof the most commonly used structured financialproducts and is intensively studied in quanti-tative finance. By setting the asset pool intodifferent tranches, it effectively works out andredistributes credit risks and returns to meet therisk preferences for different tranche investors.The copula models of various kinds are normallyused for pricing CDOs, and the Monte Carlosimulations are required to get their numericalsolution. Here we implement two typical CDOmodels, the single-factor Gaussian copula modeland Normal Inverse Gaussian copula model,and by applying the conditional independenceapproach, we manage to load each model ofdistribution in quantum circuits. We then applyquantum amplitude estimation as an alternativeto Monte Carlo simulation for CDO pricing. Wedemonstrate the quantum computation resultsusing IBM Qiskit. Our work addresses a usefultask in finance instrument pricing, significantlybroadening the application scope for quantumcomputing in finance.
Quantum computing for finance applications is anemerging field with quickly growing popularity. Thefinance industry involves various numerical and an-alytical tasks, e.g. , derivative pricing, credit rating,forex algorithm trading, and portfolio optimization, etc. .They all demand heavy quantitative work, and theimproved calculation speed and precision would bringsignificant social value. Quantum computing aims atthese very targets . Early studies focused on improv-ing finance models with basic quantum mechanics .Schrodinger equations and Feynman’s path integralwere suggested to solve stochastic differential equa-tions for pricing interest rate derivatives , and Heisen-berg uncertainty principle was used to interpret theleptokurtic and fat-tailed distribution of stock pricevolatilities . Recent studies tend to utilize quantum ad-vantages as a faster computing machine. Algorithmsthat can be implemented in quantum circuits, such asamplitude estimation , quantum principle componentanalysis (PCA) , quantum generative adversarial net-work (QGAN) , the quantum-classical hybrid variationalquantum eigensolver (VQE) and quantum-approximate-optimization-algorithm (QAOA) , spring up and begin tobe applied to various financial quantitative tasks . Within all sectors of quantitative finance, the Monte-Carlo simulation always plays a significant role , asonly a few stochastic equations for derivative pricing havefound analytical solutions , while most can only besolved numerically by repeating random settings a greatmany times in an uncertainty distribution ( e.g. normal orlog-normal distribution), which therefore consumes muchtime. The quantum amplitude estimation (QAE) algo-rithm was raised in 2002. It is newly suggested as apromising alternative to the Monte Carlo method, as itshows a quadratic speedup comparing to the latter . Sofar, applications of QAE for option pricing and creditrisk analysis have been demonstrated.Considering the wide use of Monte Carlo simulationand the large variety of pricing models, the involvementof quantum techniques in finance is still at its infancy.Credit derivatives are frequently mentioned financial in-struments because of the strong demand for tackling de-fault risks in the finance industry. Collateral debt obliga-tion (CDO) is a multi-name credit derivative backed ona pool of portfolios of defaultable assets (loans, bonds,credits etc.). CDO then packages the portfolio into sev-eral tranches with different returns and priorities to sufferthe default loss . CDO can effectively protect the seniortranche from the loss, but too many default events in thepool would still make the CDO collapsed, which was thecase during the subprime financial crisis in 2008. Manyvoices were then made for improving the CDO pricingmodel and strengthening regulations in various aspects.Nonetheless, the CDO itself is a useful credit instrumentthat can work out and redistribute credit risks in a veryquantitative way, and it is still widely studied in quan-titative finance. So far, however, the implementation ofcomplex credit instruments like CDO in quantum algo-rithms has never been reported.In this work, we present the first quantum circuit im-plementation for CDO pricing using IBM Qiskit. To ad-dress the correlations among a large number of assets inthe CDO pool, we use both the common Gaussian cop-ula model and an improved model, the Normal InverseGaussian copula model that can interpret the skew-ness and kurtosis of the real markets which the Gaussiandistribution cannot portray . We follow a conditionalindependence approach to load the correlated distribu-tions in the quantum circuits, and then use quantumcomparators and QAE algorithm to calculate the lossesin different tranches. We demonstrate the quantum com-putation results for a CDO that matches the classical a r X i v : . [ q -f i n . R M ] A ug FIG. 1:
The CDO tranche structure.
The CDO comprisesEquity Tranche (consisting of unrated or lowly rated securi-ties), Mezzanine Tranche (consisting of intermediately ratedsecurities) and Senior Tranche (consisting of highly rated se-curities), and the tranches have a sequence to bear the loss.
Monte Carlo method, suggesting a promising approachfor pricing various derivatives.
I. THE CDO STRUCTURE AND PRICINGMODELSA. The CDO tranche structures
The CDO pool is normally divided into three tranches:the Equity, Mezzanine and Senior Tranche. As shownin Fig.1, when defaults occur, the Equity Tranche in-vestors bear the loss first, then the Mezzanine Trancheinvestors if the loss is greater than the first attachmentpoint. Only when the loss is greater than the secondattachment point, will the Senior Tranche investors losemoney. Therefore, Senior Tranche has the priority ofreceiving principle and interest payment, and the bestprotection from risk while having the lowest return.Let K L k and K U k denote the lower and upper attach-ment point for Tranche k , respectively. When defaultsoccur, the buyer of the Tranche k will bear the loss in ex-cess of K L k , and up to K U k − K L k . Let L denote the totalloss for the portfolio and the loss suffered by the holdersof Tranche k be L k . As there are various default scenariosunder some uncertainty distribution, we evaluate the ex-pectation value of the tranche loss E [ L k ] for each Tranche k : E [ L k ] = E [ min [ K U k − K L k , max (0 , L − K L k )]]. Thenwe can get the fair spread for this tranche denoted as r k : r k = E [ L k ] N k = E [ L k ] K U k − K L k (1)where N k is the notional value of Tranche k of the port-folio, which can be calculated by K U k − K L k . To arriveat a fair price of a CDO, the return for investors of eachtranche should be consistent with the expected loss theinvestors would bear. Therefore, such a fair spread isconsidered as the return for this tranche. B. The conditional independence approach
Usually the pool in CDO is a portfolio of correlatedassets. Their default events are not independent, whichcan be modeled using the single-factor Gaussian copula.Meanwhile, through years’ practice on the Gaussianmodel, it is found not to well portray the phenomenain real CDO markets, e.g. , the ‘correlation smile’ . In2005, the Normal Inverse Gaussian (NIG) model was in-troduced to CDO pricing. In fact, price volatilities inderivative markets seldom show perfect Gaussian distri-bution. NIG can flexibly introduce a target skewness andkurtosis which the Gaussian model cannot achieve .Explanation for NIG distribution and its probability den-sity function (pdf) can be seen in Appendix I.For either the Gaussian copula or NIG copula model,both of them can use the conditional independenceapproach originally developed by Vaˇs´ıˇcek for themultivariate distribution problems. Consider a portfo-lio that comprises n assets, each with a default risk X i ,and a correlation γ i with the systematic risk Z . Themutually independent latent variables W i can be used: W i = γ i Z + (cid:112) (1 − γ i ) X i , where γ i s are correlation pa-rameters that can be obtained by calibrating the marketdata. Let p i be the original default probability for asset i that is uncorrelated to Z . Via detailed derivation , thedefault probability under the influence of Z follows: p i ( z ) = F ( F − ( p i ) − √ γ i z √ − γ i ) (2)This Eq.(2) is derived for very general scenarios . F stands for the distribution function of Z , which can beany continuous and strictly increasing distribution func-tion, and in this content, they are Gaussian for the Gaus-sian copula model or NIG for the NIG copula model. F − stands for the inverse of distribution F .Using this conditional independence model, given p i ( z )and λ i , the loss that would incur for asset i when defaulthappens, the expected total loss would be: E [ L ] = (cid:90) ∞−∞ n (cid:88) i =1 λ i p i ( z ) f ( z )d z (3)where f ( z ) is the pdf for an uncertainty model. For Gaus-sian distribution with a variance σ , integrating z from − σ to 3 σ would cover 99.73% of the distribution. II. QUANTUM CIRCUIT CONSTRUCTIONA. Load the correlating default risks
To apply quantum computation for CDO pricing, theprimary task is to load the default risk for each asset ofthe portfolio into the quantum circuit. Either the Gaus-sian or NIG model can be loaded in the circuit following
HH ........ ....H H H ...... L X L Z Load Z distribution ...
S C&R n z n x ........ s PP +− gg + objective qubit FIG. 2:
The quantum circuit framework.
The quantumcircuit firstly uses operator L X to load the assets with non-correlated independent default risks, and then loads z distri-bution and uses operator L Z to address the correlation amongasset default risks, and further sums up the total portfolio lossusing operator S . Then it comes to the comparator operator C and the piecewise linear rotation operator R to calculatethe tranche loss, which is related to P , the probability of theobjective qubit at the state | (cid:105) after rotation. a previous treatment for the conditional independencemodel.We firstly load the uncorrelated default events X i ( i =1 , , ..., n ) using linear Y -rotation gate. The default prob-ability p i for each asset i can be obtained from its histor-ical performance. The operator L X for loading X i pre-pares the state (cid:112) − p i | (cid:105) + (cid:112) p i | (cid:105) , so the probabilityfor state | (cid:105) encodes the default probability p i .Then we need to load the default correlation using op-erator L Z . The linear Y -rotation gate R Y ( z ) for opera-tor L Z satisfies that the probability for state | (cid:105) wouldbe p i ( z ). The expression for p i ( z ) as a function of z andthe correlation-free p i just follows Eq.(2), which derivesthe slope and offset for the rotation gate for operator L Z , i.e. sin − ( (cid:112) p i ( z )) = slope ∗ z + offset . Meanwhile, we use n z qubits to discretize the distribution to 2 n z slots, anduse affine mapping to encode the influence of z valuein the quantum circuit. See derivation of slope and offsetand explanation on affine mapping in Appendix II.We need to note that the Gaussian or NIG distributionfor systematic risk Z has to be loaded before operator L Z . The probability of z , i.e. the y axis of Gaussiandistribution function, can be loaded using the built-incodes of uncertainty model in Qiskit, and we contributethe similar codes for NIG distribution.Furthermore, we set an operator S to sum up the lossdue to all default events in this asset pool. The sum ofloss equals to (cid:80) a i λ i , where λ i is the loss given defaultfor asset i , and a i is 1 if asset i default and is 0 viseversa. The probability for a i = 1 is just p i ( z ) given byoperator L Z . The maximum loss would be (cid:80) λ i when allassets default, so ensuring the maximum loss to be en-coded needs n s qubits that (cid:80) λ i ≤ n s −
1. The operator S uses 2 n s qubits following the previous design of sumoperator . See Fig.2 for the quantum circuit frameworkfor loading default events and summing up the loss. B. Calculate tranche loss
We use the comparator operator C L k ( k =1, 2 and 3)to compare the sum of loss with the fixed lower attach-ment point K L k for each Tranche k . The comparator hasbeen used to compare the underlying asset value with thestriking price for option pricing in a recent work , wherethe detailed quantum comparator circuit has been given.The operator C L k would flip a qubit from | (cid:105) to | (cid:105) if L ( z ), the sum of loss under the systematic risk Z , ishigher than K L k , and would keep | (cid:105) otherwise. Mean-while, the objective qubit will also rotate its state un-der the control of the comparator ancilla qubit. Thepiecewise linear rotation operator R would then rotatecos( g ) | (cid:105) + sin( g ) | (cid:105) into: √ − P | (cid:105) + √ P | (cid:105) (SeeFig.2). g π − c , and c is a scaling factor. The set-ting ensures the probability after rotation remains in themonotonously increasing region when total loss increases. P , the probability at state | (cid:105) , is found to have a rela-tionship with the tranche loss: P = ( 12 − c ) + 2 cK U k − K L k ( E [ L k ]) (4)where E [ L k ] is the expectation of loss for a certaintranche k , for instance, the loss for Equity Tranche whensetting K L and K U . See detailed derivation for Eq.(4)in Appendix III.Then it comes to the issue how to read the value of P . We know that measuring an arbitrary state willmake it collapse to the orthogonal basic state. Instead,Quantum Amplitude Estimation (QAE) would be ableto map the amplitude to be estimated ( P in this case)to the discretized value using m additional qubits, andhence read P . The QAE circuit involving controlled ro-tation gates and inverse quantum Fourier transform isexplained in Appendix IV. QAE has been demonstratedas a good alternative to Monte Carlo simulation for fi-nance pricing . In this work, QAE that estimates P allows us to obtain the CDO tranche loss and return. III. RESULT ANALYSIS
We consider an example to show the pricing for CDOtranches. As listed in Table I, the CDO pool has fourassets, each showing a default probability p i , a sensitivityto the systematic risk γ i and a loss given default λ i .The CDO is divided into three tranches: the Equity,Mezzanine and Senior Tranches. Values for the lowerattachment point K L k and upper attachment point K U k for three tranches are provided in Table II.For this task, we need n x = 4 qubits to represent thefour assets in operator L X , and n z = 4 qubits in operator L Z to make 2 = 16 slots for the uncertainty distributionof systematic risk Z . We implement Gaussian (Fig.3a)and NIG (Fig.3b) distribution for Z .For NIG distribution, by setting the parameters givenin Appendix I, it shows a skewness of 1 and kurtosis FIG. 3:
Systematic risk distribution and the tranche loss functions ( a-b ) The probabilities of 2 n z different z valuesusing n z qubits follow the Gaussian distribution(mean=0, variance=1) in ( a ) and the NIG distribution (skewness=1, kurtosis=6,mean=0, variance=1) in ( b ). For both distributions, the range is from -3*variance to 3*variance. ( c-e ) The tranche loss as afunction of cumulative loss for ( c ) Equity Tranche, ( d ) Mezzanine Tranche, and ( e ) Senior Tranche. In the white box in ( c-e ),the first, second and third array respectively show the breakpoints, slopes and offsets for this tranche. T is the attachmentpoint between Equity and Mezzanine Tranche, while T is that between Mezzanine and Senior Tranche. of 6, which are consistent with a real CDO market .Comparing with Gaussian distribution, this is narrowerand centered to the left. TABLE I:
The relevant parameters for each asset.
Asset i λ i p i γ i The attachment points for each tranche.
Tranche Name Lower K L k Upper K U k Equity 0 1Mezzanine 1 2Senior 2 7
The step after loading distribution is to calculate thecumulative loss. The maximal loss is (cid:80) λ i = 7 for thisportfolio. Therefore, we can use n s = 3 qubits to encodethe total loss in the weighted sum operator S .The pricing of the tranche loss is similar to the calloption pricing, where there is a linear ‘payoff function’that goes up from zero after the option striking price or,for the CDO tranche, the attachment point. The trancheloss as a function of the total cumulative loss is given inFig.3c-e for this specific example. See Appendix V forhow to set the input parameters ( e.g. , ‘breakpoint’) forthe piecewise linear rotation function. We then use QAE to estimate P and convert it tothe tranche loss according to Eq.(4). We use the QASMcloud backend that is in the Noisy Intermediate-ScaleQuantum (NISQ) environment. Fig.4 demonstrates thecalculated tranche loss for an NIG distribution (Fig.3b)using quantum computation with m = 4 qubits in QAEand the classical Monte Carlo method. The results fromthe two approaches match well. When z follows Gaussiandistribution (Fig.3a), consistent results have also beenobtained, as shown in Fig. A2 in the appendix.With the calculated tranche loss, we can price the CDOtranche return according to Eq.(1). The notional value N for the Equity Tranche, Mezzanine Tranche and Se-nior Tranche is 1, 1, and 5, respectively, by calculating K U − K L for each tranche. We get the tranche loss fromquantum computation and calculate the tranche returnfor Equity, Mezzanine and Senior Tranche to be 49.9%,48.1% and 1.8%. The low return for the Senior Trancheis consistent with the practice in reality. Such a low valueis firstly due to the last sequence to bear the loss, and sec-ondly owing to its large notional value, which is normallyabove 80% of the sum for the three tranches. See morediscussion on tranche return in practice in Appendix VI. IV. DISCUSSION AND CONCLUSION
The CDO is a relatively advanced and complex struc-tured finance product, and the credit market plays a sig-nificant role in the finance industry. Therefore, despitethere were some disputes on CDOs during the 2008 fi-nancial crisis, CDOs are still widely studied products inquantitative finance, and are being improved with var-ious financial models. In this work, we implement thenormal inverse Gaussian model that is now regarded asadvantageous over the Gaussian model. There is also thevariance gamma model that was first applied to optionpricing and later found to be a good model for CDOpricing . Such improved models can also be calculatedvia quantum computation. -0.5 0.0 0.5 1.0 1.5 P r obab ili t y Equity Tranche Loss -0.5 0.0 0.5 1.0 1.5
Mezzanine Tranche Loss P r obab ili t y -5.0 -2.5 0.0 2.5 5.0 P r obab ili t y Senior Tranche Loss a b c FIG. 4:
CDO tranche loss with Z under the NIG dis-tribution. The calculated loss for ( a ) Equity Tranche, ( b )Mezzanine Tranche, and ( c ) Senior Tranche. Z follows theNIG distribution depicted in Fig.3b. Blue bars indicate quan-tum computation results with m = 4 for QAE, and red dashedlines indicate Monte Carlo results with repeating 10000 times. Note that the quantum adaption of generative adver-sarial network has now been considered as an effectiveway to load any distribution in quantum circuits andcan be applied to more finance models. Besides, the pa-rameter shift rule has been raised to solve the issue ofencoding gradients in quantum circuits, which facilitatesthe mapping of machine learning techniques in quantumalgorithms. Furthermore, the trendy variational quan-tum algorithms that are suitable for NISQ environment,and the alternative approach using quantum annealing ,may work on a large variety of optimization tasks in fi- nance. In all, there’s much room to explore for quantumcomputation in finance applications. Acknowledgments
H.T. thanks Prof. Stephen Schaefer’s previous helpfor studies on fixed income and interest rate deriva-tive at London Business School. The authors thankJian-Wei Pan for helpful discussions. This research wassupported by National Key R&D Program of China(2019YFA0308700, 2017YFA0303700), National NaturalScience Foundation of China (61734005, 11761141014,11690033, 11904229), Science and Technology Commis-sion of Shanghai Municipality (STCSM) (17JC1400403),and Shanghai Municipal Education Commission (SMEC)(2017-01-07-00-02- E00049). X.-M.J. acknowledges addi-tional support from a Shanghai talent program.
Author Contributions.
H.T. and X.-M.J. conceivedand supervised the project. H.T. and A.P. designed thescheme. A.P. wrote the Qiskit code. H.T. did MonteCarlo simulation. H.T., A.P., L.F.Q., T.Y.W., J.G. andX.M.J. analyzed the data and presented the figures. H.T.wrote the paper, including the appendix, with input fromall the other authors.
Competing Interests.
The au-thors declare no competing interests.
Data Availabil-ity.
The data that support the plots within this paperand other findings of this study are available from thecorresponding author upon reasonable request. ∗ Electronic address: [email protected] † Electronic address: [email protected]. Or´us, R., Mugel, S., & Lizaso, E. Quantum computing forfinance: Overview and prospects.
Rev. in Phys. , 100028(2019).2. Baaquie, B. E. Quantum finance: Path integrals andHamiltonians for options and interest rates.
CambridgeUniversity Press (2007).3. Zhang, C., & Huang, L. A quantum model for the stockmarket.
Physica A: Statistical Mechanics and its Applica-tions , 5769-5775 (2010).4. Meng, X., Zhang, J. W., & Guo, H. Quantum Brownianmotion model for the stock market.
Physica A: StatisticalMechanics and its Applications , 281-288 (2016).5. Brassard, G., Hoyer, P., Mosca, M., & A. Tapp, A. Quan-tum Amplitude Amplification and Estimation.
Contempo-rary Mathematics , 53-74 (2002).6. Lloyd, S., Mohseni, M., & Rebentrost, P. Quantum prin-cipal component analysis.
Nat. Phys. , 631-633 (2014).7. Lloyd, S.,& Weedbrook, C. Quantum Generative Adver-sarial Learning. Phys. Rev. Lett. ,040502 (2018).8. Peruzzo, A., McClean, J., Shadbolt, P., Yung, M.-H.,Zhou, X.-Q., Love, P. J., Aspuru-Guzik, A., & O’Brien,J. L. A variational eigenvalue solver on a quantum proces-sor.
Nat. Commun. , 4213 (2014). 9. Farhi, E., Goldstone, J., & Gutmann, S. A QuantumApproximate Optimization Algorithm. arXiv Preprint arXiv:1411.4028 (2014).10. Rebentrost, P., Gupt,B., & Bromley, T. B. Quantum com-putational finance: Monte Carlo pricing of financial deriva-tives. Phys. Rev. A , 022321 (2018).11. Stamatopoulos, N., Egger, D. J., Sun, Y., Zoufal, C., Iten,R., Shen, N., & Woerner, S. Option Pricing using QuantumComputers. Quantum , 291 (2020).12. Egger, D. J., Guti´errez, R. G., Mestre, J. C., & Woerner,S. Credit Risk Analysis using Quantum Computers. arXivPreprint arXiv:1907.03044 (2019).13. Woerner, S., & Egger, D. J. Quantum risk analysis. npjQuantum Info. , 15 (2019).14. Martin, A., Candelas, B., Rod´ıguez-Rozas, A., Mart´ın-Guerrero, J. D., Chen, X., Lamata, L., Or´us, R., Solano,E., & Sanz, M., Towards Pricing Financial Derivativeswith an IBM Quantum Computer. arXiv Preprint arXiv:1904.05803 (2019).15. Zoufal, C., Lucchi, A., & Woerner, S. Quantum Genera-tive Adversarial Network for Learning and Loading Ran-dom Distributions. npj Quantum Info. , 103 (2019).16. Goodfellow, I. J., Pouget-Abadie, J., Mirza, M., Xu, B.,Warde-Farley, D., Ozair, S., Courville, A., & Bengio, Y.Generative Adversarial Nets. Conference proceedings of advances in neural information processing systems
Options futures and other derivatives.
PearsonEducation India (2003).18. Tuckman, B., & Serrat, A.
Fixed income securities:tools for today’s markets, 3rd Edition.
John Wiley & Sons(2012).19. Chacko, G., Sj¨oman, A., Motohashi, H., & Dessain,V
Credit Derivatives, Revised Edition: A Primer onCredit Risk, Modeling, and Instruments.
Pearson Educa-tion (2016).20. Black, F. & Scholes, M. The pricing of options and cor-porate liabilities.
J. Political Econ. , 637654 (1973).21. Merton, R.C. Theory of rational option pricing. The BellJournal of economics and management science , 141-183(1973).22. Li, D. X. On default correlation: A copula function ap-proach. The Journal of Fixed Income , 43-54 (2000).23. Barndorff-Nielsen, O. Hyperbolic Distributions and Dis-tributions on Hyperbolae. Scandinavian Journal of Statis-tics , 151-157 (1978).24. Barndorff-Nielsen, O. Normal Inverse Gaussian Distribu-tions and Stochastic Volatility Modelling. ScandinavianJournal of Statistics , 1-13 (1997).25. Gu´egan, D., & Houdain, J. Collateralized Debt Obliga-tions pricing and factor models: a new methodology usingNormal Inverse Gaussian distributions. Note de RechercheIDHE-MORA No. 007-2005, ENS Cachan. (2005).26. Kalemanova, A., Schmid, B., & Werner, R. The NormalInverse Gaussian distribution for Synthetic CDO Pricing. The Journal of Derivatives
Spring , 80-93 (2007).27. Schl¨osser, A. Normal Inverse Gaussian Factor CopulaModel.
Pricing and Risk Management of Synthetic CDOs
International Journal of Theoretical &Applied Finance , 5 (2014).29. Vaˇs´ıˇcek, O. Probability of loss on loan portfolio. Technicalreport, KMV Corporation (1987).30. Vaˇs´ıˇcek, O. The distribution of loan portfolio value.
Risk , 160-162 (2002).31. Madan, D. B., Carr, P. P., & Chang, E. C. The VarianceGamma Process and Option Pricing. European FinanceReview , 79-105 (1998).32. Moosbrucker, T. Pricing CDOs with Correlated VarianceGamma Distributions. working paper, Centre for FinancialResearch, University of Cologne. (2006).33. Li, J., Yang, X. D., Peng, X. H., & Sun, C. P. Hybridquantum-classical approach to quantum optimal control. Phys. Rev. Lett. , 150503 (2017).34. Schuld, M., Bergholm, V., Gogolin, C., Izaac, J., & Kil-loran, N. Evaluating analytic gradients on quantum hard-ware.
Phys. Rev. A , 032331 (2019).35. Cohen, J., Khan, A., & Alexander, C. Portfolio Optimiza-tion of 40 Stocks Using the DWave Quantum Annealer. arXiv Preprint arXiv: 2007.01430 (2020). Appendix I The model of Normal Inverse Gaussian Distribution
The Normal Inverse Gaussian (NIG) distribution mixes the normal Gaussian distribution and inverse Gaussiandistributions[Ref A1].The word ‘inverse’ in the name needs to be explained. While normal distribution reflects the location distributionunder Brownian motion at a certain time, the inverse Gaussian distribution shows the time distribution when theBrownian motion moves to a certain location, so inverse suggests an inverse way in viewing location and time.Firstly, for a random variable Y that has inverse Gaussian distribution, its density of function is of the form: f IG ( y ; α, β ) = α √ πβ y − / exp( − ( α − β y) β y ) (A1)Then if a random variable X satisfies the following requirement with parameters α , β , µ and δ , it follows the NormalInverse Gaussian distribution N IG ( x ; α, β, µ, δ ): X | Y = y ∼ N ( µ + βy, y ) Y ∼ IG ( γ, γ ) with γ = (cid:112) α − β (A2)The full expression of the probability density function for NIG distribution is a bit complicated: N IG ( x ; α, β, µ, δ ) = a ( α, β, µ, δ ) q ( x − µδ ) − K ( δαq ( x − µδ )) e βx (A3)where the function q follows: q ( x ) = √ x , and a is the function with variables α , β , µ and δ : a ( α, β, µ, δ ) = π − α exp( δ (cid:112) α − β − βµ ) (A4)with parameters satisfying: 0 < = | β | < α and δ >
0, and K is the first index of the Bessel function of the third kind: K ( x ) = x (cid:90) ∞ exp( − xt ) (cid:112) t − dt (A5)The parameter α is related to steepness, β to symmetry, µ to location and δ to scale. In order to realize the NIGdistribution (mean=0, variance=1, skewness=1, and kurtosis=6), the parameters need to be set as follows[Ref A2]: α =-1.6771, β =0.75, µ =-0.6, and δ = 1 . \ qiskit \ aqua \ components \ uncertainty models). Appendix II Derive the linear rotation parameters for operator L Z and explain affine mapping A linear rotation function would use the state qubit | x (cid:105) to work on the target qubit | (cid:105) : | x (cid:105) | (cid:105) → | x (cid:105) (cos( slope ∗ x + offset ) | (cid:105) +sin( slope ∗ x + offset ) | (cid:105) ) (A6)After the rotation, the probability at the state | (cid:105) would be the probability that we are interested in, and it is the z -tuned default probability p i ( z ) in operator L Z : (cid:112) p i ( z ) = sin( slope ∗ z + offset ) (A7)so sin − (cid:112) p i ( z ) has to be expressed in the form of slope ∗ z + offset to get the slope and offset for the linear rotationquantum gate in operator L Z .Combining the expression for p i ( z ) in Eq.(2) using the conditional independence model, we have:sin − (cid:112) p i ( z ) = sin − (cid:115) F ( F − ( p i ) − √ γ i z √ − γ i ) (A8)where F is the cumulative distribution function (CDF), and we denote ψ = F − ( p i ) √ − γ i , so the above equation can besimply expressed as: sin − (cid:112) p i ( z ) = sin − (cid:115) F ( ψ − √ γ i z √ − γ i ) (A9)The first order Taylor’s theorem can be expressed as: g ( x ) = g ( a ) + g (cid:48) ( a )( x − a ) (A10)where g (cid:48) ( a ) is the first derivative of the function g ( x ).Now let x be ψ − √ γ i z √ − γ i , and a be ψ , then x − a is √ γ i z √ − γ i . We know that for Taylor expansion, x − a has to bea very marginal value. √ γ i z √ − γ i satisfies this. From the example given in the Result Analysis Section, γ is generallybelow 0.2 and z follows a Gaussian or NIG distribution, ranging between [ − ,
1] with a value close to zero at mosttimes. Therefore, Taylor’s theorem for this task is reasonable.In this scenario, function g ( x ) is sin − (cid:112) F ( x ). It is a composite function, so the chain rule should be consideredfor the derivative function g (cid:48) ( a ). sin − ( h ) (cid:48) = √ − h and ( √ h ) (cid:48) = − √ h apply for any variable h , and note that F isthe cumulative distribution function (CDF), the derivative of F would just be the probability density function (pdf)which is normally denoted as f . Now g (cid:48) ( a ) reads: g (cid:48) ( a ) = − f ( ψ )2 (cid:112) − F ( ψ ) (cid:112) F ( ψ ) (A11)Therefore we get: sin − (cid:115) F ( ψ − √ γ i z √ − γ i ) =sin − (cid:112) F ( ψ ) + − f ( ψ )2 (cid:112) − F ( ψ ) (cid:112) F ( ψ ) √ γ i z √ − γ i (A12)One can now very easily get the expression for slope and offset : slope = −√ γ i √ − γ i f ( ψ ) (cid:112) − F ( ψ ) (cid:112) F ( ψ ) (A13) offset = 2arcsin( (cid:112) F ( ψ )) (A14)In the L Z operator, with n z qubits, we can truncate and discretize the distribution to 2 n z slots. For instance,with 3 qubits, z ranges between 0 and 7 that corresponds to -3 σ to 3 σ of the Gaussian distribution; for z = 4 =1 ∗ + 0 ∗ + 1 ∗ −
1, Qubit 1 and Qubit 3 turn on their controlled rotation gate R Y (2 ) and R Y (2 ), whileQubit 2 does not switch on its R Y (2 ). Through such a so-called affine mapping[Ref A3], the influence of z value onthe linear rotation is encoded in the quantum circuit. Appendix III Piecewise linear rotation for objective qubit
We use the comparator operator C L k ( k =1, 2 and 3) to compare the sum of loss with the fixed lower attachmentpoint K L k for each Tranche k . The comparator has been used to compare the underlying asset value with the strikingprice for option pricing in a recent work[Ref A4], where the detailed quantum comparator circuit has been given.The operator C L k would flip a qubit from | (cid:105) to | (cid:105) if L ( z ), the sum of loss under the systematic risk Z , is higherthan K L k , and would keep | (cid:105) otherwise. For simplicity, we just discuss C L generally that can later apply to all C L k sby just setting K L to K L k .Meanwhile, the objective qubit will also rotate its state under the control of the comparator ancilla qubit. Thecomparator C L and a piecewise linear rotation operator R map | ψ (cid:105) | (cid:105) [cos( g ) | (cid:105) + sin( g ) | (cid:105) ] into: (cid:40) | ψ (cid:105) | (cid:105) [cos( g ) | (cid:105) + sin( g ) | (cid:105) ] if L ( z ) ≤ K L | ψ (cid:105) | (cid:105) [cos( g + g z ) | (cid:105) + sin( g + g z ) | (cid:105) ] if L ( z ) > K L (A15)and the probability at state | (cid:105) that can be measured using QAE would be expressed as follows: P = (cid:88) L ( z ) ≤ K L f ( z )sin ( g ) + (cid:88) L ( z ) >K L f ( z )sin ( g + g z ) (A16)where f ( z ) shows the probability of a systematic risk value z distributed following a certain probability densityfunction f . g is set to be π − c . g z is a linear function that contains the total loss L ( z ) and scaling factor c . g z canbe implemented using controlled Y-rotations, and it’s mapped to integer value z ∈ { , . . . , n z − } . Note that thereis an upperbound breakpoint K U as well, so we set another comparator operator C U that encodes K U , and g z finallyreads as: g z = 2 c min ( L ( z ) , K U ) − K L K U − K L (A17)With such settings g z would be in the range { , c } , and by choosing a small scaling parameter c , we can ensuresin( g + g c ) in a monotonously increasing regime. Given that for marginal x : sin ( x + π ) = x + + O ( x ), theexpression for probability P can be further derived: P = (cid:88) L ( z ) ≤ K L f ( z )( 12 − c )+ (cid:88) L ( z ) >K L f ( z )( 12 − c + 2 c min ( L ( z ) , K U ) − K L K U − K L )= ( 12 − c ) + (cid:88) L ( z ) >K L f ( z )(2 c min ( L ( z ) , K U ) − K L K U − K L )= ( 12 − c ) + 2 cK U − K L ( E [ L tranche ]) (A18)where E [ L tranche ] is the expectation of loss for a certain tranche, for instance, setting K L and K U to be K L and K U , we get the loss for the Equity Tranche. Therefore, when we have obtained P from the QAE circuit, we wouldbe able to get the tranche loss and return. Appendix IV Quantum Amplitude Estimation
Given a Boolean function f : to find an x ∈ X where X → { , } such that f ( x ) = 1, we can denote N as thenumber of inputs on which f takes the value 1, and it can be written as N = |{ x ∈ ( X | f ( x ) = 1 }| .If we have a classical probabilistic algorithm P that outputs a guess on input x , the solution to instance x can befound by repeatedly calling P and X . If X ( x, P ( x )) = 1 with probability p >
0, we have to repeat the process p times on average. + ψ - ψψ γ a θ a θ γθ − a Final State ψ S S Initial State a b
FIG. A1:
Theoretical framework for QAE. ( a ) Illustration for angle rotation by operator Q . ( b ) The quantum circuit forquantum phase estimation. H denotes the Hadamard gate. QF T − denotes the inverse quantum Fourier transform. a b c -0.5 0.0 0.5 1.0 1.50.00.20.40.60.81.0 P r obab ili t y Equity Tranche Loss -0.5 0.0 0.5 1.0 1.50.00.20.40.60.81.0 P r obab ili t y Mezzanine Tranche Loss -4 -2 0 2 40.00.20.40.60.81.0 P r obab ili t y Senior Tranche Loss
FIG. A2:
CDO tranche loss with Z under the Gaussian distribution. The calculated loss for ( a ) the Equity Tranche,( b ) the Mezzanine Tranche, and ( c ) the Senior Tranche. The systematic risk Z follows the Gaussian distribution (mean=0,variance=1). Blue bars indicate the results from quantum computation with m = 4 for QAE, and red dashed lines indicate theMonte Carlo results with repeating 10000 times. Suppose given a unitary transformation, A , which is a quantum algorithm, unlike making measurement in the caseof classical P algorithm, this produces a quantum superposition state of the “desired” result that X ( x, P ( x )) = 1 and“undesired” result that X ( x, P ( x )) (cid:54) = 1. Then amplitude estimation is the problem of estimating a , the probabilitythat a measurement of | ψ (cid:105) yields a good solution. It is sufficient to evaluate A and X in an expected number of timesthat is proportional to √ a . To explain amplitude estimation, the quantum state after unitary transformation A can be expressed as a linearcombination of ψ + and the orthogonal ψ − : A| (cid:105) = | ψ (cid:105) = − i √ (cid:0) e i θ a | ψ + (cid:105) + e − i θ a | ψ − (cid:105) (cid:1) (A19)so the success probability “ a ” is converted to the solution of angle θ a that decides the eigenvalue for unitary trans-formation A .Apply an operator Q to A : Q = AS A † S ψ (A20)where S = 1 − | (cid:105) (cid:104) | , and S ψ = 1 − | ψ (cid:105) | (cid:105) (cid:104) | (cid:104) ψ | . As illustrated in Fig.A1a, an initial state first rotates along ψ by the operator S ψ , and then rotates along ψ + by the operator S . Therefore, the angle between the arbitrarily setinitial state and the final state after the operator Q becomes 2 θ a . In this case when the initial state is A , operator Q just shifts from | ψ (cid:105) for 2 θ a .The task of finding the eigenvalue for quantum state | ψ (cid:105) of the unitary transformation A can be fulfilled by QuantumPhase Estimation that requires another register with m additional qubits. As shown in Fig.A1b, the phase estimationquantum circuit comprises Hadamard gates, controlled-rotation operators and an inverse quantum Fourier transform( QFT − ) operation.Firstly, the Hamamard gates prepare the m qubits in the uniform superposition: | (cid:105) ⊗ m | ψ (cid:105) → √ m m − (cid:88) j =0 | j (cid:105) | ψ (cid:105) (A21)As has been mentioned, the operator Q essentially causes a Y-rotation of angle 2 θ a , i.e. , Q = R y (2 θ a ). In this phaseestimation circuit, the many controlled-rotation operators Q j satisfies: Q j = R y (2 jθ a ), and they turn the above Eq.(A21) into: 1 √ m m − (cid:88) j =0 e iθ a j | j (cid:105) | ψ (cid:105) (A22)1Applying the QFT − operation, we can reverse the action on vector | j (cid:105) to that on | θ a (cid:105) : QFT − (cid:16) √ m m − (cid:88) j =0 e iθ a j | j (cid:105) | ψ (cid:105) (cid:17) = 1 π | θ a (cid:105) | ψ (cid:105) (A23)By taking measurement on the register of m qubits. we can get the approximation of θ a . This is done by obtainingthe measured integer y → { , , , . . . m − } . Taking M = 2 m , then θ a can be approximated as (cid:101) θ a = yπ/M , whichyields (cid:101) a , the approximation of the aforementioned probability a : (cid:101) a = sin (cid:16) yπM (cid:17) ∈ [0 ,
1] (A24)satisfying the following inequality: | a − (cid:101) a | ≤ πM + π M = O ( M − ) (A25)with probability at least M . Comparing with the O ( M − ) convergence rate of the classical Monte Carlo method,the quantum amplitude estimation method converges faster with a quadratic speed-up. Appendix V Set input parameters for the built-in piecewise linear rotation function in Qiskit
For the C & R part of the quantum circuit shown in Fig.2 of the main text, we can use the built-in code named‘PwlObjective’ for piecewise linear rotation function that includes the comparator C , and the piecewise linear rotator R . The built-in function uses the ‘breakpoints’ array to record the attachment points, and uses the ‘slopes’ and‘offsets’ arrays in which slope k and offset k correspond to these for the line segment between breakpoint k − k . Note the offset is the y -axis value for the starting point of the line segment, instead of the intercept byextending the line segment to the y axis. The breakpoints, slopes and offsets for the tranche loss function are shownin each figure in Fig. 3c-e, which can be very straightforwardly calculated. These are used as the input parametersfor the built-in piecewise linear rotation function. Appendix VI Discussion on tranche return in reality
It’s worth noting that the returns for Equity and Mezzanine Tranche in the case study of the main text are a bittoo high, comparing to the custom returns that would be around 15-25% and 5-15% for the Equity and MezzanineTranche, respectively[Ref A5]. It’s partially because that default probabilities p i s are a bit high. One more reason isthat we ignore the recovery rate of the asset in order to focus on the essential structure. The recovery rate η , which isgenerally set as 40%, means that when asset defaults, some values can be recovered by ways like selling real estates toget funds to compensate investors. Then the maximum loss would equal to the total notional value multiplies (1- η ).In this example, the loss given default λ to λ would become 1.2, 1.2, 0.6 and 1.2, while the tranche attachmentpoints keep unchanged. This would bring down the tranche loss. Supplementary References [A1] Schl¨osser, A. Normal Inverse Gaussian Factor Copula Model.
Normal Inverse Gaussian Factor Copula Model. arXiv Preprint , arXiv:1907.03044 (2019).[A4] Stamatopoulos, N., Egger, D. J., Sun, Y., Zoufal, C., Iten, R., Shen, N., & Woerner, S. Option Pricing usingQuantum Computers.