Pricing Bermudan Interest Rate Swaptions via Parallel Simulation under the Extended Multi-factor LIBOR Market Model

Bermudan swaptions are options on interest rate swaps which can be exercised on one or more dates before the final maturity of the swap. Because the exercise boundary between the continuation area and stopping area is inherently complex and multi-dimensional for interest rate products, there is an inherent “tug of war” between the pursuit of calibration and pricing precision, tractability, and implementation efficiency. After reviewing the main ideas and implementation techniques underlying both single- and multi-factor models, we offer our own approach based on dimension reduction via Markovian projection. Specifically, on the theoretical side, we provide a reinterpretation and extension of the classic result due to Gyöngy which covers non-probabilistic, discounted, distributions relevant in option pricing. Thus, we show that for purposes of swaption pricing, a potentially complex and multidimensional process for the underlying swap rate can be collapsed to a one-dimensional one. The empirical contribution of the paper consists in demonstrating that even though we only match the marginal distributions of the two processes, Bermudan swaptions prices calculated using such an approach appear well-behaved and closely aligned to counterparts from more sophisticated models.

Bermudan Swaptions in the LIBOR Market Model

We present a parallel algorithm and its multi-threaded implementation for computing lower and upper bound prices of multi-asset Bermudan options. Our baseline sequential algorithm follows Longstaff and Schwartz’s least-squares Monte Carlo method in computing the lower bound and Andersen and Broadie’s simulation-based procedure with sub-optimality checking for the upper bound. The parallel implementation uses POSIX Threads for thread manipulation and Intel’s MKL functions for random number generation and linear algebra operations. Tests were made on Intel x86 multi-core processors using the same option examples as the previous work, and the runtimes of the same computations were reduced from minutes to a few seconds.

We present an experimental study of parallel biconnected components algorithms employing several fundamental parallel primitives, e.g., prefix sum, list ranking, sorting, connectivity, spanning tree, and tree computations. Previous experimental studies of these primitives demonstrate reasonable parallel speedups. However, when these algorithms are used as subroutines to solve higher-level problems, there are two factors that hinder fast parallel implementations. One is parallel overhead, i.e., the large constant factors hidden in the asymptotic bounds; the other is the discrepancy among the data structures used in the primitives that brings non-negligible conversion cost. We present various optimization techniques and a new parallel algorithm that significantly improve the performance of finding biconnected components of a graph on symmetric multiprocessors (SMPs). Finding biconnected components has application in fault-tolerant network design, and is also used in graph planarity testing. Our parallel implementation achieves speedups up to 4 using 12 processors on a Sun E4500 for large, sparse graphs, and the source code is freely-available at our Web site.

The single-site dynamics are a canonical class of Markov chains for sampling from high-dimensional probability distributions, e.g. the ones represented by graphical models. We give a simple and generic parallel algorithm that can faithfully simulate single-site dynamics. When the chain asymptotically satisfies the $\ell_p$-Dobrushin's condition, specifically, when the Dobrushin's influence matrix has constantly bounded $\ell_p$-induced operator norm for an arbitrary $p\in[1,\infty]$, the parallel simulation of $N$ steps of single-site updates succeeds within $O\left({N}/{n}+\log n\right)$ depth of parallel computing using $\tilde{O}(m)$ processors, where $n$ is the number of sites and $m$ is the size of graphical model. Since the Dobrushin's condition is almost always satisfied asymptotically by mixing chains, this parallel simulation algorithm essentially transforms single-site dynamics with optimal $O(n\log n)$ mixing time to RNC algorithms for sampling. In particular we obtain RNC samplers, for the Ising models on general graphs in the uniqueness regime, and for satisfying solutions of CNF formulas in a local lemma regime. With non-adaptive simulated annealing, these RNC samplers can be transformed routinely to RNC algorithms for approximate counting. A key step in our parallel simulation algorithm, is a so-called procedure, which tries to simultaneously couple all distributions over the same sample space. We construct such a universal coupling, that for every pair of distributions the coupled probability is at least their Jaccard similarity. We also prove this is optimal in the worst case. The universal coupling and its applications are of independent interests.

In this paper we empirically compare a wide range of different term structure models when it comes to the pricing and, in particular, hedging of caps and swaptions. We analyze the influence of the number of factors on the hedging and pricing results, and investigate which type of rate or derivative price data should be used to estimate the model parameters to obtain the best hedging and pricing results. We use on interest rates, and cap and swaption prices from 1995 to 1999. The empirical results on the hedging of caps and swaptions show that, if the number of hedge instruments is equal to the number of factors, the multi-factor models outperform one-factor models in hedging caps and swaptions. However, if one uses a large set of hedge instruments, one-factor models perform as well as multi-factor models. In terms of pricing, we find that models with two or three factors imply better out-of-sample predictions of cap and swaption prices than one-factor models. Also, estimation on the basis of current derivative prices leads to more accurate out-of-sample prediction of cap and swaption prices than estimation on the basis of interest rate data.

Linear algebra operations play an important role in scientific computing and data analysis. With increasing data volume and complexity in the "Big Data" era, linear algebra operations are important tools to process massive datasets. On one hand, the advent of modern high-performance computing architectures with increasing computing power has greatly enhanced our capability to deal with a large volume of data. One the other hand, many classical, deterministic numerical linear algebra algorithms have difficulty to scale to handle large data sets. Monte Carlo methods, which are based on statistical sampling, exhibit many attractive properties in dealing with large volume of datasets, including fast approximated results, memory efficiency, reduced data accesses, natural parallelism, and inherent fault tolerance. In this dissertation, we present new Monte Carlo methods to accommodate a set of fundamental and ubiquitous large-scale linear algebra operations, including solving large-scale linear systems, constructing low-rank matrix approximation, and approximating the extreme eigenvalues/ eigenvectors, across modern distributed and parallel computing architectures. First of all, we revisit the classical Ulam-von Neumann Monte Carlo algorithm and derive the necessary and sufficient condition for its convergence. To support a broad family of linear systems, we develop Krylov subspace Monte Carlo solvers that go beyond the use of Neumann series. New algorithms used in the Krylov subspace Monte Carlo solvers include (1) a Breakdown-Free Block Conjugate Gradient algorithm to address the potential rank deficiency problem occurred in block Krylov subspace methods; (2) a Block Conjugate Gradient for Least Squares algorithm to stably approximate the least squares solutions of general linear systems; (3) a BCGLS algorithm with deflation to gain convergence acceleration; and (4) a Monte Carlo Generalized Minimal Residual algorithm based on sampling matrix-vector products to provide fast approximation of solutions. Secondly, we design a rank-revealing randomized Singular Value Decomposition (R3SVD) algorithm for adaptively constructing low-rank matrix approximations to satisfy application-specific accuracy. Thirdly, we study the block power method on Markov Chain Monte Carlo transition matrices and find that the convergence is actually depending on the number of independent vectors in the block. Correspondingly, we develop a sliding window power method to find stationary distribution, which has demonstrated success in modeling stochastic luminal Calcium release site. Fourthly, we take advantage of hybrid CPU-GPU computing platforms to accelerate the performance of the Breakdown-Free Block Conjugate Gradient algorithm and the randomized Singular Value Decomposition algorithm. Finally, we design a Gaussian variant of Freivalds’ algorithm to efficiently verify the correctness of matrix-matrix multiplication while avoiding undetectable fault patterns encountered in deterministic algorithms.

The support of range, segment and rectangle queries are fundamental problems in computational geometry, and have extensive applications in many domains. Despite significant theoretical work on these problems, efficient implementations can be complicated, and most implementations do not have useful theoretical bounds. In this paper, we focus on simple and efficient parallel algorithms and implementations for range, segment and rectangle queries, which have worst-case bounds in theory and good performance in practice, both sequentially and in parallel. We propose to use a framework based on the abstract data type augmented map, to model the problems. Based on the augmented map interface, we develop both multi-level tree structures and sweepline algorithms supporting range, segment and rectangle queries in two dimensions. For the sweepline algorithms, we also propose a parallel paradigm and show corresponding cost bounds. Theoretically, the construction algorithms of all of our data structures are work-efficient and highly parallelized.We have implemented all the data structures described in the paper, ten in all, using a parallel augmented map library. Based on the library, each data structure only requires about 100 lines of C++ code. We test their performance on large data sets (up to 108 elements) and a machine with 72-cores (144 hyperthreads). The parallel construction achieves 32-68x speedup, and the speedup numbers for queries are up to 126-fold. Sequentially, each of our implementations outperforms the CGAL library by at least 2x in both construction and queries. Our sequential implementation has approximately the same construction time as the R-tree in the Boost library, but has significantly better query performance (1.6-1200x). We believe this paper provides the most comprehensive experimental study of data structures on range, segment and rectangle queries, both in parallel and sequential setting.

Unlocking massively parallel spectral proper orthogonal decompositions in the PySPOD package

A new method to retrieve aerosol optical thickness from satellite images on a parallel system

Parallel simulation has been an active research area for more than a decade, and several parallel simulation algorithms have been proposed. To evaluate parallel simulation environments, there is a need for a common benchmark suite. Such benchmarks should allow the designer of simulation kernels to: (i) evaluate how efficiently the simulation kernel runs on specific architectures; and (ii) evaluate how simulation problems scale on the kernel. A vast majority of benchmarks suggested in the literature focus on the latter problem. The first requirement is however equally important, as such benchmarks are needed to implement efficient simulation systems. In this report, we advocate an incremental benchmark methodology which is primarily intended for performance tuning. The benchmark suite is based on a small set of ping models with which it is possible to effectively isolate and estimate various overheads, contention and latencies encountered in simulation kernels. The use of the benchmark suite is illustrated by its application to performance tuning and evaluation of a time warp kernel.

Purpose. Modeling systems and programming platforms provide ample opportunities for the use of statistical tools in research activities. Since the normal distribution is one of the most common distribution laws, the criterion for checking the sample for normality is in high demand among statistical assessment tools, among which the Epps-Pulley test has the status as one of the most powerful tests to check the deviation of the distribution from the normal one. There are a number of implementations of this test in the R and Python languages. However, this test is not implemented in one of the most popular Matlab modeling software. Thus, the purpose of this study is to develop a software implementation of the Epps-Pulley criterion in the Matlab environment and verify the correctness of the performed calculations.Materials and Methods. We implemented the calculation of Epps-Pulley statistics by two methods – classical, using cycles, and matrix-vector, using linear algebra operations. The classical method requires calculating the intermediate values necessary to obtain the criterion statistics using two independent cycles, the second cycle being a double one, in which one cycle is nested into the other. The matrix-vector method requires fewer code by performing calculations using linear algebra operations on matrices and vectors. We obtained critical statistical values for the sample size from 8 to 1000 observations with two-dimensional linear interpolation of tabular values. We used an approximation by a beta function of the third kind for a sample of over 1000 elements.Results. An assessment of the computational efficiency of the methods showed that the cyclic approach is about three times higher than the matrix-vector approach in terms of consumed time, which is presumably due to the processing of insignificant elements in triangular matrices when performing component-by-component operations. The correctness of the software implementation of the Epps-Pulley criterion was tested on several examples, which confirmed the compliance of the calculated values of the criterion statistics, as well as the critical values of statistics, with known data. We carried out a criterion statistical evaluation based on the empirical values of the error of the first kind. We obtained the error values correspondence to the specified significance levels. We performed comparative estimates of the Epps-Pulley test with the Anders-Darling and Shapiro-Wilk tests in terms of the criterion empirical power and tabulated the evaluation results. We published the software implementation of the Epps-Pulley test on the MATLAB Central Internet resource and for free use.Conclusion. We developed software implementation of the Epps-Pulley criterion as a new research tool that was previously unavailable in the Matlab modeling environment. We used the time spent on calculations to make a reasonable choice of the calculation algorithm for the criterion statistics. We confirmed correctness of the calculation algorithms by a set of selective checks and statistical estimates that showed the compliance with well-known theoretical provisions.

Deep MVA: Deep Learning for Margin Valuation Adjustment of Callable Products

In the framework of a multifactor LIBOR market model (LMM), this paper presents a new approach for finding an approximate distribution of constant maturity swap (CMS) rates under the terminal martingale measure. With this approach, we derive an analytical pricing formula for CMS range notes, which is both intuitive and tractable. Many exotic CMS rate derivatives are widely traded in the marketplace or embedded in structure notes.

We present a semi-static replication algorithm for Bermudan swaptions under an affine, multi-factor term structure model. In contrast to dynamic replication, which needs to be continuously updated as the market moves, a semi-static replication needs to be rebalanced on just a finite number of instances. We show that the exotic derivative can be decomposed into a portfolio of vanilla discount bond options, which mirrors its value as the market moves and can be priced in closed form. This paves the way toward the efficient numerical simulation of xVA, market, and credit risk metrics for which forward valuation is the key ingredient. The static portfolio composition is obtained by regressing the target option’s value using an interpretable, artificial neural network. Leveraging the universal approximation power of neural networks, we prove that the replication error can be arbitrarily small for a sufficiently large portfolio. A direct, a lower bound, and an upper bound estimator for the Bermudan swaption price are inferred from the replication algorithm. Additionally, closed-form error margins to the price statistics are determined. We practically study the accuracy and convergence of the method through several numerical experiments. The results indicate that the semi-static replication approaches the LSM benchmark with basis point accuracy and provides tight, efficient error bounds. For in-model simulations, the semi-static replication outperforms a traditional dynamic hedge.