Directeur De

Raphael Douady

La Présidente

Dominique Guegan

Professeur M Rama Cont

Professeur Rapporteur M Andrew Mullhaupt

Rapporteur M Pierre

Del Moral

M Ricadro Perez-Marco

Professeur Jean-Paul Laurent

Particle Methods in Finance

Liste des figures

1.1 HMC sampling from Bivariate Gaussian distribution . . . . . . . . . . . . . . 2.1 Left : Positive exposure profiles. Right : EPE(Blue), ENE(Red) . . . . . . . . . 2.2 The paths of CIR process (left) and survival function (right) . . . . . . . . . . 2.3 OU Jump-diffusion process with exponentially tempered α-stable jumps (left), Survival probability (right) with α = 0.8 . . . . . . . . . . . . . . . . . . . . . . 2.4 OU Jump-diffusion process with exponentially tempered α-stable jumps (left), Survival probability (right) with α = 1.5 . . . . . . . . . . . . . . . . . . . . . . 2.5 PMC simulation of OU Jump-diffusion process with exponentially tempered α-stable jumps (left), Survival probability (right) with α = 0.8 . . . . . . . . . 2.6 PMC Euler, Euler jump-adapted, Milstein jump-adapted strong errors, α = 0.8, #(particles)= 1000. . . . . . . . . . . . . . .

Résumé

This thesis consists of two parts, namely rare event simulation and a homotopy transport for stochastic volatility model estimation.

Particle methods, that generalize hidden Markov models, are widely used in different fields such as signal processing, biology, rare events estimation, finance and etc. There are a number of approaches that are based on Monte Carlo methods that allow to approximate a target density such as Markov Chain Monte Carlo (MCMC), sequential Monte Carlo (SMC). We apply SMC algorithms to estimate default probabilities in a stable process based intensity process to compute a credit value adjustment (CVA) with a wrong way risk (WWR). We propose a novel approach to estimate rare events, which is based on the generation of Markov Chains by simulating the Hamiltonian system. We demonstrate the properties, that allows us to have ergodic Markov Chain and show the performance of our approach on the example that we encounter in option pricing.

In the second part, we aim at numerically estimating a stochastic volatility model, and consider it in the context of a transportation problem, when we would like to find "an optimal transport map" that pushes forward the measure. In a filtering context, we understand it as the transportation of particles from a prior to a posterior distribution in pseudotime. We also proposed to reweight transported particles, so as we can direct to the area, where particles with high weights are concentrated. We showed on the example of Stein-Stein stochastic volatility model the application of our method and illustrated the bias and variance.

Keywords : Hamiltonian flow Monte Carlo, Particle Monte Carlo, Sequential Monte Carlo, Monte Carlo, rare events, option pricing, stochastic volatility, optimal transport

Table des matières

Motivation

Throughout last 20 years development of computation power allowed to use sophisticated Monte Carlo methods in signal-processing, rare events estimation, computational biology, queuing theory, computational statistics and etc. In finance we have to deal with large dimensionality of problems, where other techniques due to some constraints, such as curse dimensionality, computational burden makes us to look for alternative numerical techniques. Particle methods are a broad class of interacting type Monte Carlo algorithms for simulating from a sequence of probability distributions satisfying a nonlinear evolution equation. These flows of probability measures can always be interpreted as the distributions of the random states of a Markov process whose transition probabilities depends on the distributions of the current random states ( [6], [5]).

Particle methods have a variety of applications in finance : rare events and stochastic volatility estimation. For instance, portfolio managers need to estimate a large portfolio loss for risk management. In banking industry, banks need to compute default probabilities, so as to compute different value adjustments and to comply with financial regulation. In insurance industry, companies are interested in estimating ruin probabilities in a given time horizon. All the above mentioned cases are the examples, where rare event simulation is applied.

Stochastic volatility models are widely used in the financial industry. In fact, while realized volatility could be hedged away by trading other options, stochastic volatility models are needed to model the dynamics of implied volatilities, which will provide their user with simple break-even accounting conditions for the Profit and Loss of a hedged position [1].

The current thesis consists of two main parts, namely computing rare events via simulation and general stochastic volatility estimation.

Credit Risk Estimation

Since the credit crisis 2007-2009, the importance of counterparty credit risk for regulators increased dramatically. According to Basel III regulation [33], banks are required to hold a regulatory capital based on CVA charges against each of their counterparties. The is already a number of article on CVA valuation, the most common Credit Value Adjustment (CVA) formula ( [5], [22], [22], [START_REF] Pykhtin | A Guide to Modelling Counterparty Credit Risk[END_REF] and [START_REF] Prisco | Modeling Stochastic Counterparty Credit Exposures for Derivatives Portfolios[END_REF]) is given by : CVA = (1 -R)

∫ T 0 D t E[V + t |τ = t ]d PD(t ) (1.1)
where R is a recovery rate, D t -risk-free discount rate, PD(t ) probability of default up to time t and V t -the value of an underlying asset.

Two problems related to CVA computation that usually arise : incorporation of wrong way risk (WWR) into the value of CVA and a high computational burden. WWR is a statistical dependence between the exposure and a counterparty's credit risk (CCR). Another challenge related to CVA is the computation of default probability. There are two main approaches in default probability computation : structural models (first-passage approach) [3] and reduced form credit models ( [19], [26] and [28]).

Structural models

The first passage approach is modelled like in a barrier option pricing, and formulated in the following way. Assume that the default barrier b is constant valued. Then the default time τ is a continuous random variable (r.v.) values in ]0, +∞[ and given by τ = inf{t > 0 : V t < b} (1.2) then the first passage probability is given by : P(T) = P(M T < b) = P(inf

s≤T ( µs + σ D W s ) < log(b/V 0 )) (1.3)
where M T = inf s≤T V s .

Default intensity

The main difference between reduced form and structural models is in the fact that, in the latest one does not need any economic model of firm's default, i.e. defaults are exogenous.

D t = 1 τ<t = { 1 if τ ≤ t 0 otherwise (1.4)
Observe that default intensity is an increasing process, such that the conditional probability at time t that the firm defaults at time s ≥ t is at least as big as the process D t .

A process with such property is called submartingale. The Doob-Meyer decomposition theorem says that we can isolate upward trend from D. This very important result says that there exist an increasing process D τ starting at zero such that D -A τ becomes a martingale. The unique process A τ counteracts the upward trend in D, and it is called a compensator. It describes the cumulative and conditional likelihood of default, which is parameterized by non-negative process λ.

A τ t = ∫ t τ 0 λ s d s = ∫ t 0 λ s 1 τ>s d s (1.5)
where λ t describes conditional default rate for small interval ∆t and τ > t . λ∆t approximates the probability that default occurs in the interval ]t , t + ∆t].

Rare events

Rare events simulation is an important field of research in computational and numerical probability. It has a wide range of applications starting from catastrophe theory to finance, as an example one can consider pricing of a barrier option pricing or a credit risk estimation.

Let us consider the problem of rare event estimation, where the probability z = P(A) is very small. The crude Monte Carlo (CMC) estimates z through the proportion of times on the simulated points attain rare event area A over M independent trials :

b z M = 1 M M ∑ m=1 1 ξ (m) ∈A (1.6)
with variance σ = √ z(1z). In "easy problems" we can use the central limit theorem (CLT), which asserts that

1 p M (z M -z) ----→ M→∞ σZ (1.7)
where Z ∼ N (0, 1) is a standard Gaussian distribution.

In a rare-events simulation we are not very concerned about the absolute error, but instead we would like to measure via relative error(RE) the precision of our estimation with respect to some "true" value of quantity we are computing. Relative error exposes the problems one encounters by using CMC in rare event setting :

RE = Z √ z(1 -z) M 1/2 z = Z √ 1 -z Mz ∝ Z √ Mz ---→ z→0 ∞ (1.8)
The approximation (1.8) shows that CMC requires n ≫ 1/z. Lets show this on a simple example, assume that we would like to estimate rare event probability z = 10 -7 , RE is targeted at 0.1 with a 95% confidence interval, then we have 1.96 p 10 -7 M ≤ 0.1

(1.9) from the equation above we see that we need at least M ≥ 3.84 × 10 -9 samples to have a 10% RE.

Assume that we have a sequence of random variables {ξ(a) (m) } M m=1 . Define the rare event set as :

A a = {x ∈ R M , f (x) > a} (1.10)

then we can define a probability of interest and its estimator as

z(a) = E[1 ξ(a)∈A a ] and b z(a) = 1 M M ∑ m=1 1 ξ (m) (a)∈A a (1.11)
We can show that if our estimator has a bounded RE, then the desired precision is independent of that the rarity of the set A a . (1.12)

From the equation (1.12) we that the set will not depend on the "rarity" of the set A a , if we impose the following condition :

Pricing in Partial Observation Models

Stochastic volatility models are the one of examples of the use of partial observation models (hidden Markov models) in finance, i.e. the situation, when one can observe the prices but not the dynamics of the stochastic volatility. In the last chapter we will show the application of different numerical approaches, such as particle filtering and homotopy transport on the example of a barrier option pricing in Stein-Stein stochastic volatility models. For example, in [1] author showed that the price of a barrier option is mostly dependent on the dynamics of the at-the-money skew conditional on price hitting the barrier. A stochastic volatility model for barrier options would need to provide a direct handle on this precise feature of the dynamics of the volatility surface so as to appropriately reflect its Profit and Loss impact in the option price. On of the first applications in option pricing could be found in [6]. Other applications and analysis of hidden Markov models could be found in [8] and [START_REF] Macdonald | Hidden Markov models for time series : an introduction using R[END_REF].

Sampling Methods

Monte Carlo

The application of CMC in a rare event setting was already demonstrated in previous sections. These methods are generally used to compute integral that can not be calculated analytically or very difficult to compute. For example, we would like to calculte the following expectation :

E P [h(X)] = ∫ X h(x)P(d x) (1.14)
All we need is to generate M independent and identically distributed (i.i.d.) random variables {ξ (m) } M m=1 according to the law P. Using the law of large numbers (LLN) we have an unbiased estimator

1 M M ∑ m=1 δ ξ (m) (x)h(x) a.s. ----→ M→∞ E P [h(X)]
(1.15)

As we saw in the previous section, often we can not efficiently sample from the measure P. Those problems arises in rare events setting and sampling from fat-tailed distributions, then we have to use advanced Monte Carlo techniques, such as an importance sampling (IS), control variates, a stratified sampling and etc.

Importance Sampling

One of the ways to deal with rare events probabilities is IS Monte Carlo method. The idea is to change a measure from P, where samples rarely reach rare event sets, to the measure Q, so that we can sample at low cost in the new measure .

E P [h(X)] = ∫ X h(x)P(d x) = ∫ X h(x)ω(x)Q(d x) = E Q [h(X)ω(X)] (1.16) M ∑ m=1 δ ξ (m) (x)ω(x)h(x) a.s. ----→ M→∞ E Q [h(X)ω(X)] (1.17)
The disadvatage of this method is in the fact, that we do not know the explicit form of the Radon-Nykodim derivative d P dQ = ω. This method is unfeasible in many examples, if we exclude some simple cases.

Markov Chain Monte Carlo

Markov Chain Monte Carlo (MCMC) method allows to approximate a target measure π, by constructing an ergodic Markov Chain, that admits the target measure as a stationary law. We do not need to know the law explicitly as in the case of the IS, but instead we construct a Markov kernel K , that leaves the target measure π invariant. There are two very popular MCMC sampling techniques : Metropolis-Hastings(MH) and Gibbs algorithms. One can show that the Gibbs algorithm is a special case of MH. The idea behind Metropolis-Hastings algorithm is to propose new set of candidates, and accept them, so that the transition kernel K is left invariant with respect to the target, where Markov transition kernel is given by :

K (x, d y) = α(x, y)Q(x, d y) + ( 1 - ∫ α(x, z)Q(x, d z) ) δ x (d y) (1.18) CHAPITRE 1. INTRODUCTION
and MH transition kernel K is reversible, i.e. πK = π. By recursion we can show that the chain {X n } n≥0 follows the law π n = π 0 K n . Using reversibility property of the kernel K , we can have convergence of the chain X n to the target measure π, given that the kernel K has contraction :

π n -π = (π 0 -π)K n ∀n > 0 (1.19)
If the chain X l is irreducible, or in other words it admits a unique invariant measure, then the paths of the chain satisfies ergordic theory :

lim n→∞ 1 n n H +n ∑ l =n H +1 h(X l ) = E π [h(X)] (1.20)
where n H is the number of iterations needed for "burn-out period", so that the chain leaves its initial law, and converges to the invariant measure π. Usually one needs to make high number of iterations, so that MH algorithm will approach the target density.

The performance of MCMC algorithms, like in other Monte Carlo sampling algorithms depends on the experience and ability to tune and optimize it. In general, the quality of final samples depends on the kernel K and how fast it is able to explore the state space. In chapter 3 we will show, that the transition kernel represented by a metropolized Hamiltonian dynamics allows to explore it fast.

Algorithm 1 : MH MCMC algorithm

1 Initialization : N -#(time steps), π 0 -initial measure 2 X 0 ∼ π 0 ; 3 for n = 1,...,N do 4 Generate x * from Q(•, X n-1 ) ; 5
Compute the weight :

a(X n , x * ) = 1 ∧ π(x * )q(x * , X n ) π(X n )q(X n , x * ) assuming that a(x, y) = 0 if π(x)q(x, y) = 0 ; 6 Draw u ∼ U nif(0, 1) ; 7 if u < a then 8 Set X n+1 = x * ;

Particle Methods and Sequential Monte Carlo

Interacting particle system(IPS) is very power method, that is applied in non-linear filtering problems, rare event estimation and many other applications. IPS allows to overcome the problem that we face in importance sampling Monte Carlo technique, where we needed to know explicitly the importance sampling measure in order to sample from it. On contrary, the importance measure is approximated by a collection of trajectories of a Markov process {X n } n≥0 weighted by a set of potential functions ω n . IPS is related to the unnormalized Feynman-Kac models. In the next section we give a brief overview of these stochastic models, for more details we refer to monographs [6] and [5].

Feynman-Kac Approximations

We give some key notations that are used to describe Feynman-Kac models, that we use in the next chapters. First, we define normalised and unnormalized Feynman-Kac measures (η n , γ n ) for any bounded function f by the formulae :

γ N ( f ) = E[ f (X N ) N-1 ∏ n=0 ω n (X n )] and η N ( f ) = γ N ( f ) γ N (1 X ) (1.21)
Since non-negative measures (η n ) n≥0 satisfy for any bounded function f the recursive linear equation :

γ n ( f ) = γ n-1 (Q n ( f )) we have γ N ( f ) = ∫ ... ∫ f (X N )η 0 (d x 0 ) ∏ Q n (x n , d x n+1 ) (1.22)
where {Q n } n≥0 are unnormalised transition kernels, for example, we can choose them of the following form

Q n (x n , d x n+1 ) = ω n (x n )k n (x n , d x n+1
), where {k n } n≥0 is a sequence of elementary transitions.

The following definition was given in [6], for a linear semigroup Q m,n , 0 ≤ m ≤ n associated with a measure γ n and defined by : Q m,n = Q m+1 , ..., Q n . For any bounded function f and (m, n) ∈ N :

γ n ( f ) = γ m Q m ...Q n-1 f and η n ( f ) = γ m Q m ...Q n-1 ( f ) γ m Q m ...Q n-1 (1 X ) = η m Q m ...Q n-1 ( f ) η m Q m ...Q n-1 (1 X ) (1.23)
with a convention Q m , ..., Q n = Id , if m > n by definition of a normalised measure η n and a semigroup Q m,n we readily obtain

η n+1 ( f ) = η n Q n ( f ) η n Q n (1 X ) (1.24)
Let us consider the example of rare events estimation, assume that we have a sequence of rare event sets {A p } 0≤p≤n , X is a random variable on the probability space (Ω, P, F ). One of the application in barrier option pricing or credit risk estimation is computation of conditional expectation

η N ( f ) = E[ f (X N )|X p ∈ A p ] (1.25)
Feynman-Kac interpretation for any bounded function f and potential function

ω n (x) = 1 A n (x), using the fact that 1 A n 1 A n+1 = 1 A n (x), since A n+1 ⊂ A n η N ( f ) = 1 P(A N ) ∫ 1 A N (x)h(x)P(d x) = E[ f (X N ) ∏ N-1 l =1 ω l (X l )] E[ ∏ N-1 l =1 ω l (X l )] = E[ f (X N )1 A N (X)] E[1 A N (X)] (1.26)

Hidden Markov Models

As we mentioned in the motivation section, hidden Markov Models (HMM) can be formulated in the context of a filtering problem. Under the filtering problem one assumes that we have a couple of processes in discrete time (X n , Y n ) n≥0 , where (X n ) n≥0 is a sequence of hidden variables and (Y n ) n≥0 is partially observed data.

X n+1 |X n = x n ∼ k n (x n , d x n+1 ) (1.27) Y n |X n = x n ∼ ρ n (x n , d y n ) (1.28)
where k n is an elementary transition kernel, and ρ n is a likelihood function. Using Bayes formula for a bounded function f we can derive

E[ f (X 0 , ..., X n )|Y 0 , ..., Y n-1 ] = ∫ ... ∫ f (X 0 , ..., X n ) ∏ ρ p (X p , Y p )ν(d x 0 )k 0 (X 0 , d x 1 )...k n-1 (X n-1 , d x n ) ∫ ... ∫ ∏ ρ p (X p , Y p )ν(d x 0 )k 0 (X 0 , d x 1 )...k n-1 (X n-1 , d x n )
(1.29) Motivating example for studying and using HMM in finance is stochastic volatility model. It is defined by the following system :

d Y t = µ(X t )d t + σ Y (X t )d W t d X t = β(X t )d t + σ X (X t )d B t (1.30)
Using Euler discretization scheme we have

Y n = Y n-1 + µ(X n )∆t + σ Y (X n ) p
∆tZ (1) n

X n+1 = X n + β(X n )∆t + σ X (X n ) p ∆tZ (2) n (1.31)
where Z (1) and Z (2) are two independent standard Gaussian r.v., ∆t is a discretization step.

In this case, the unnormalised transition kernel has the following form

Q n (x n , d x n+1 ) = ρ n (x n , d y n )k n (x n , x n+1 ).
For the sake of illustration, assume that we have constant vola-

tilities σ X (X n ) = σ X , σ Y (X n ) = σ Y ,
and other parameters are given by ∆t = 1, β(X n ) = 0, µ(X n ) = X n and Y n-1 = 0, then the likelihood is given by

ρ n (x n , d y n ) = 1 p 2πσ Y exp ( - (y n -x n ) 2 2σ 2 Y ) (1.32)
and an elementary Markov transition kernel

k n (x n , d x n+1 ) = 1 p 2πσ X exp ( - (x n+1 -x n ) 2 2σ 2 X ) (1.33)
If we exclude toy examples, there does not exist the exact simulation method for such type of models. Sequential Monte Carlo based particle algorithms allow to approximate normalised measures η n .

Sequential Monte Carlo

Assume that on the measurable space (X, X ), there exists an unnormalized transition

Q n : X n × X n → R + , which is absolutely continuous with respect to a kernel K n : X n × X n+1 → R + , i.e. Q n (x n , •) ≪ K n (x n ,
•) for all x n ∈ X n , then we can define an importance weight function ω n :

ω n (x n , x n+1 ) = d Q n (x n , •) d K n (x n , •) (x n+1 ) (1.34)
Sequential Monte Carlo (SMC) produces a fixed number of weighted samples (particles).

At each time instant we generate a couple of particles and their corresponding weights {ξ m n , ω (m) n } M m=1 , where M is a fixed number of particles. If at time n ∈ N, the set of weighted particles {ξ (m) n , ω (m) n } M m=1 approximates measure η n , then using IS method, we can also approximate η n+1 by a couple of the set of particles and their corresponding weights. {ξ (m) n+1 , ω (m) n+1 } M m=1 . This procedure is called a sequential importance sampling (SIS), which CHAPITRE 1. INTRODUCTION allows us recursively construct the measure η N from the initial measure η 0 ∈ P (X) and the set of importance weights {ω (m) n } M m=1 .

Algorithm 2 : Sequential Importance Sampling algorithm

1 Initialization : M -#(simulations), N -#(time steps), η 0 -initial measure

2 for m = 1,...,M do 3 ξ (m) 0 ∼ η 0 ; 4 ξ (m) 1 ∼ K 0 (ξ (m) 0 , •) ; 5 ω (m) 1 = ω (m) 0 (ξ (m) 0 , ξ (m) 1 ) 6 end 7 for n = 1,...,N do 8 for m = 1,...,M do 9 Generate ξ (m) n from K n (•, ξ (m) n-1 ) and set b ξ (m) n = ( b ξ (m) n , ξ (m) n-1 ) ; 10
Compute the weight :

ω (m) n ( b ξ (m) n ) ; 11 end 12 end
In some cases, we can observe the weight degeneracy, when the variance of importance weights increase over time, as the fact that most of weights are negligible. To overcome this problem, in [24] authors proposed to resample particles at each iteration, which was called sequential importance resampling (SIRS) or bootstrap algorithm. There was an extensive research to optimize the bootstrap using, for example, effective sample size (ESS), we refer to [27], [16] and [14] for details.

In [6] and in [5] it was shown, that SMC estimators are unbiased, i.e. we can prove that sup

|| f ||≤1 ¯¯E [ 1 M M ∑ m=1 f (ξ m N ) -η N ( f ) ]¯¯≤ c(N) M (1.35)
where c(N) is some positive constant, whose values depend on time horizon N.

Algorithm 3 : Sequential Importance Resampling algorithm

1 Initialization : M -#(simulations), N -#(time steps), η 0 -initial measure

2 for m = 1,...,M do 3 ξ (m) 0 ∼ η 0 ; 4 ξ (m) 1 ∼ K 0 (ξ (m) 0 , •) ; 5 ω (m) 1 = ω (m) 0 (ξ (m) 0 , ξ (m) 1 ) 6 end 7 for n = 1,...,N do 8 for m = 1,...,M do 9 if n < N then 10 Resample using probability weight : W n ( b ξ (m) n ) = ω n ( b ξ (m) n ) 1 M ∑ M j =1 ω n ( b ξ ( j ) n ) ; 11 end 12 Generate ξ (m) n from K n (•, ξ (m) n-1 ) and set b ξ (m) n = ( b ξ (m) n , ξ (m) n-1 ) ; 13
Compute the weight :

ω (m) n ( b ξ (m) n ) ; 14 end 15 end

Hamiltonian Flow Monte Carlo

Hamiltonian Flow Monte Carlo (HFMC) methods came from statistical physics, where the computation of macroscopic properties requires sampling phase space configurations distributed according to some probability measure. One of its first application in sampling from high-dimensional distribution was proposed in [15], and its differenti extensions and converging properties in [2], [12] and [6]. Sampling on Riemann manifolds using generalised HFMC was introduced in [14].

We consider the system of particles described by a position and momentum X and P respectively, that are modelled by a Hamiltonian energy H (X, P). In statistical physics, the macroscopic properties can be obtained by averaging of some function A with respect to a probability measure P describing the state of the system of particles :

E P [A(X)] = ∫ X A(x, p)P(d x, d p) (1.36)
In most cases we approximate the above expectation numerically, such that as number of iteration increases, the microscopic sequences converge to the target distribution P :

lim N→∞ 1 N N-1 ∑ i =0 A(X n , P n ) = E P [A(X)] P -a.s. (1.37)
One of the simplest examples are so called canonical measures that have the following form :

P(d x, d p) = 1 Z e -βH (x,p) d xd p, Z = ∫ X e -βH (x,p) d xd p (1.38)
The measure P is called Gibbs measure. Since we can separate potential and kinetic energies in the Hamiltonian, we can consider positions sampling, by using projection operator.Define the projection operator as pr o j •P(d x, d p) = P(d x), and consider the canonical measure, ν, which is a projection of P.

ν(d x) = 1 Z e -βΨ(x) d x, Z = ∫ X e -βΨ(x) d x (1.39)
If the Hamiltonian is separable, then the measure P takes tensorized form, and each element of momenta follows independent Gaussian distributions. The main difficulty, that one encounters in computational statistics, biology is sampling of potential energy. The approaches to sample new configurations are as follows :

-"Brute force" sampling methods, such as rejection sampling -MCMC techniques, when we accept or reject new proposals using MH algorithm -Markovian stochastic dynamics, when we generate new samples using Langevin or generalized Langevin equations -Deterministic methods on the extended state space

Sampling with HFMC

In chapter 3 we apply HFMC algorithms to estimate rare events. There are number of challenges to compute ensemble averages in eq. (1.37), for example the dynamics of X n is not ergodic with respect to the Gibbs measure, i.e. paths are not sampled from ν. The system can contain high energy barriers that will prevent fast sampling.

CHAPITRE 1. INTRODUCTION

A large deviation theory defines rare event, when the dynamics has to overcome a potential barrier ∆Ψ, such that the exit time scales like τ ∝ exp(-β∆Ψ) (1.40) We can interpet τ in the sense that sampling from the Gibbs measure takes exponentially long time. There is the following law of large numbers.

Theorem 1.4.1 [30] Given τ sufficiently small, let the numerical flow Ξ τ be symmetric and symplectic. Then

X n+1 = (pr o j • Ξ)(X n , P n ), P n ∝ e -βH (X n ,•) (1.41) P(X n+1 = x n+1 ) = 1 ∧ e -β∆H n+1
(1.42)

defines an irreducible Markov process {X 0 , X 1 , ...} ⊂ X with unique invariant probability measure π and the property

1 N N ∑ n=1 f (X n ) → ∫ E f d π a.s., ∀X 0 ∈ X (1.43)
One can simply understand the irreducibility property as the probability of reaching any point on the configuration space is nonzero, i.e.

P(X n+1 ∈ B(x τ )|X n = x 0 ) > 0 (1.44)
holds true ∀x 0 , x τ ∈ X and any Borel set B = B(X).

In algorithm

for l f = 1,...L -1 do P (m) H ((l f + 1 2 )δ) = P (m) H (l f ) -δ 2 ∂Ψ ∂x H (x (m) H (l f )) x (m) H ((l f + 1)δ) = x (m) H (l f ) + δP (m) H ((l f + 1 2 )δ)I -1 M P (m) H ((l f + 1)δ)) = P (m) H ((l f + 1 2 )δ) -δ 2 ∂Ψ ∂x H (x (m) H ((l f + 1)δ))
end Calculate acceptance probability and set x * = x (m) H (L), P * = P (m) H (L) : Then the potential energy function and partial derivatives are given by

a = 1 ∧ e (-H (x * ,P * )+H (x (m) H ,P (m) H ))∆t (1.45) Draw u ∼ U nif(0, 1) ; if u < a then Set X (m) n+1 = x * ; else Reject,
Ψ(x) = -log(p(x)) = x T Σx 2 , ∂Ψ ∂x i = x i (1.47)
Figure 1.1 shows samples from bivariate Gaussian distribution. We see that HFMC fastly explores the state space and samples from a target distribution.

Filtering by Optimal Transport

In previous sections we showed how particle filters can be applied to approximate likelihood function by a set of weighted particles. It was shown in [1], [14], [15] that if the dimension of filtering problem increases, then the collapse of particles grows super -exponentially. In [29], authors introduced a new approach to update system's measurement, where the homotopy is formed, that gradually transforms an initial prior into a posterior density as scale parameter λ increases from 0 to 1. The idea is to choose a parameterized density and minimize the deviation between this density and homotopy, according to the measure of deviation. The approach can be interpreted as the optimal transport problem, which is defined as :

inf T E[||T (X) -X|| 2 ] s. th. Q = T # P (1.48)
That means that we try to approximate some importance measure Q by an optimal transport map T , that minimize the deviation between particles transported using the homotopy and posterior distribution. Instead of approximating the likelihood function by a set of particles sampled through importance resampling algorithm, we find a transport map that moves random variables using the following homotopy :

ψ(X t ,λ |Y t ) = 1 Z λ p(X t ,λ |Y t -1 )ρ(Y t |X t ,λ ) λ (1.49)
The homotopy continuously deforms a prior distribution p(X t |Y t -1 ) into an unnormalized posterior p(X t |Y t -1 )ρ(Y t |X t , ) as λ approaches to 1 :

ψ(X t ,λ |Y t ) λ(0→1) -----→ ψ(X t |Y t ) (1.50)
The transportation of particles in n X dimensional state space is performed by a particle flow. The idea is to find a flow of particles that correspond to the flow of probability density defined by a homotopy. We suppose that flow of particles for Bayes rule follows the following dynamics in pseudo time λ :

d X t ,λ = g (X t ,λ )d λ + η(X t ,λ )d W λ (1.51)
where g (X t ,λ ) = dX t ,λ d λ .

In chapter 4 we will show that the flow g (x, λ) can be found as a solution to Fokker-Plank partial differential equations, given that its diffusion matrix sums to zero, the flow is given by :

g (X t ,λ ) = [ ∂ 2 Ψ(X t ,λ ) ∂X 2 t ,λ ] -1 [ ∂L(X t ,λ ) ∂X t ,λ
] T (1.52) Algorithm 5 : Homotopy Transport Algorithm 1 Initialization : i = 1,...,n X -#(simulations), t = 1,...,N -#(time steps)

2 Draw {X (i ) 0 } n X i =1 from the prior p 0 (x). 3 Set {ω (i ) 0 } n X i =1 =
d λ = - [ ∂ 2 Ψ(X (i ) t ,λ ) ∂X 2 t ,λ ]-1 [ ∂L(X (i ) t ,λ ) ∂X t ,λ
] T ;

16

Transport particles according to its flow : X (i ) t ,λ = X (i )

Xt = 1 n X ∑ n X i =1 X (i ) t ,λ 20 end 21 end

Organinzation of the thesis and main contributions

Current thesis consists of three articles, uploaded to HAL and Arxiv that we plan to submit. A brief description and main contributions are presented below.

Chapter 2. (Article)

Particle L évy tempered α-stable credit risk model simulation with application to CVA with WWR computation in multi-factor BGM setting (R. Douady, Sh. Miryusupov) This chapter presents a Levy driven α-stable stochastic intensity model to estimate default probability in CVA computation. We use the fact that under certain regular conditions the L évy process's small jump part could be approximated by a Brownian motion in order to correlate market and credit risks, hence to take into account the WWR. Since default probabilities are rare events, we adapt IPS methods to simulate L évy stochastic intensity. We use three discretization schemes : Euler with constant time steps, Euler and Milstein jump-adapted time steps.

Chapter 3. (Article)

Hamiltonian flow Simulation of Rare Events (R. Douady, Sh. Miryusupov) In this chapter we present an algorithm to estimate rare events using Hamiltonian dynamics. In this approach we generate Markov Chains that converges to an invariant measure. This approach allows to have rare-event estimates that have small variability. We compare our approach with IPS on the example of Barrier option.

Chapter 4. (Article) Optimal Transport Filtering with Particle Reweighing in Finance (R. Douady, Sh. Miryusupov) For the last portion of this thesis, we will move to partiallyobserved models in order to estimate stochastic volatility. We show the adaptation of particle filter to estimate the hidden parameter, which is not observed. The oprimal transport allows the transportation of particles from prior distribution into posterior by using homotopy that gradually transforms the likelihood function, when we move in pseudo time from 0 to 1. In order to improve we propose reweighted particle transportation.

Abstract

Since the beginning of crisis in 2007 banks observed anomaly in the behaviour of the credit value adjustments (CVA) that was due to the Wrong Way Risk (WWR). The WWR is known as a statistical dependence between exposure and a counterparty credit risk (CCR).

Most of approaches of dealing with the WWR are computationally intensive that makes it hard to use in the banks.

In this article we model CVA with the WWR, where a stochastic intensity follows a Levydriven Ornstern-Uhlenbeck dynamics with jump marks that have α-stable distribution. We use the fact that small jumps of a L évy process under certain regulatory conditions could be approximated by a Brownian motion with a drift. A L évy process decomposition into a Brownian motion and a compound Poisson process allows naturally embed the WWR into the CVA values through the correlation of the Brownian motion in the dynamics of exposure and the one in the Levy processs Gaussian approximation.

To reduce the samples' variability, we used mean field interacting particle system (IPS), that allowed to reduce the variability of our CVA estimates on the one hand, and reduce the number of simulated paths on the other hand. This approach also allowed to reduce errors of sample paths using Euler with constant time step, Euler and Milstein scheme with jump-adapted time steps to several order of magnitudes compared to a crude Monte Carlo.

Our results show that the WWR risk has a huge impact on the values of CVA. We hope that our results will be a message for banks and regulators in CVA computation.

Introduction

After the crisis banks observed anomaly in the behaviour of CVA that was due to the WWR. The WWR is known as a statistical dependence between exposure and counterparty credit risk (CCR). Most of approaches of dealing with the WWR are computationally intensive that makes it hard to use in the banks. The WWR is a negative statistical dependency between exposure and a counterparty's credit quality. The CVA estimation with the WWR poses a major challenge with respect to both computational burden and tractability, consequently it is not clearly accounted in Basel III regulation.

There are different approaches that have been proposed to assess the WWR, for example in [13] they add jump into the exposure process. In [5] authors proposed the change of measure, that models the presence of the WWR. [15] proposed adjusted default probability in the dependent CVA formula and in [22] used copula method to model the dependence between default time and exposures. We focus on the dependence between the counterparty default and general market risk factors and in particular we use the properties of α-stable process to correlation the approximation of small jumps in the L évy intensity and exposure processes.

In the following sections we focus on two issues : CVA with WWR estimation and IPS interpretation of α-stable process. One the one hand we propose α-stable intensity process for default probability estimation, on the other hand we would like to see the impact of WWR/RWR on the values of CVA.

Not surprisingly, we find that WWR depends on the correlation between counterparty default risk and the credit spread of the underlying asset, in particular, the regulatory ratio of 1.4 appears to be underestimated, when the correlation is negative. This paper contributes as follows : first, we construct tempered-α stable stochastic intensity model for probability of default estimation ; second, we apply interacting particle system to estimate default probability in the rare events estimation context ; third, the IPS framework that we develop for simulation of α -stable process demonstrated reduced variability of its estimates compared to a crude Monte Carlo estimators ; forth, we apply IPS to estimate CVA with right/wrong way risk. Our results show that the WWR has a big impact on the values of CVA, taking into account that BCBSs regulatory CCR capital charges assume that WWR risk increasing by a constant factor value of CVA, our results show that this assumption has to be reviewed.

Problem Formulation and Definitions

The general formula for CVA computation is given by :

CVA := LGDE[B(0, τ)V + τ 1 τ≤T ] (2.1)
where LGD = (1 -R), and LGD=loss given default, R is a recovery rate, B(0, t ) is a discount factor for maturity t . V t is a market exposure at time t , x + = max(x, 0) and τ is counterparty's default time. The expectation is under the risk neutral measure Q.

Definition

The counterparty positive exposure of a derivative security V + t , is the nonnegative part of the difference between a security value S t minus a collateral C t :

V + t = max(S t -C t , 0) (2.2)
The expected positive exposure at each time instant t :

EPE t = E[V + t |F 0 ] (2.3) CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING
We assume that C t = 0.

Definition The default distribution function up to time t is given by :

P(τ ≤ t ) := E[1 τ≤t ] (2.4)
The density of the first default time is given by :

f τ (t ) := ∂ ∂t P(τ ≤ t ) (2.5)
Given that the default time τ and discounted positive exposures B t V + t are independent, CVA at time t can be computed as :

CVA i nd = LGD ∫ T 0 E[B t V + t ] f τ (t )d t (2.6)
BCBS in [33] defined the following formula for a potential mark-to-market losses associated with a deterioration in the credit worthiness of a counterparty :

CVA ≈ LGD mkt T ∑ i =1 max(0; e - s i -1 t i -1
LGD mkte

- s i t i LGD mkt ) | {z }
De f aul t pr obabi l i t y t er m

( EE i -1 D i -1 + EE i D i 2 ) | {z } Exposur e t er m | {z } Uni l at er al CVA (2.7)
It is a good question, whether (2.7) makes any sense, since you have two different probability measures in one formula, on the one hand CDS spreads, that are calculated using risk-neutral probability measure, on the other hand you have expected exposure, which is computed using historical risk measure.

Next question, how the WWR is treated within banking regulation ? In [16] and [15] it is indicated that Basel II defines α = 1.4, but gives banks to estimate their own α, subject to a floor of 1.2. So CVA WWR is approximated using the following formula :

CVA WWR,BIII = a × CVA (2.8)
In this article we model CVA with WWR, where stochastic intensity has Levy-driven Ornstein-Uhlenbeck dynamics. Since for a different level of α, the jump magnitude changes, we can see its impact on the CVA and, in particular, how it impacts CVA with WWR. The specific form of Levy process decomposition allows us to correlate two Brownian motions, the first one in the stochastic intensity and another one in the exposure process.

To reduce the variance and the number of generated path of MC simulations, we will use particle methods. To validate this approach we will verify the mean absolute deviation(MAD) and || • || L 2 norm of simulated paths' errors using Euler scheme with constant time steps, Euler with jump-adapted and Milstein with jump-adapted time steps using MC and PMC algorithms.

Market Risk and Exposure Profiles

BGM Model

We give some basic notions on BGM model, and remind the way one could compute the exposure profiles for multifactor BGM model [26], [4] and on the example of two factor model we will demonstrate expected exposure and expected positive exposure.

L(t , T i ) = 1 δ ( 1 P(t , T i ) -1 ) (2.9)
The forward rate form T i -1 to T i , set at time t , where t ≤ T i -1 ≤ T i is given by :

L i (t ) := L(t , T i , T i -1 ) = 1 δ ( P(t , T i -1 ) P(t , T i ) -1 ) (2.10)
From this we can derive a money market rate :

P(t , T i ) = P(t , T i -1 ) 1 + δL i (t ) = i ∏ j =1 1 1 + δL i (t ) (2.11)
Since the forward rate is a martingale under the masure Q T i , we have

E Q T i [L(t , T i ) -L i (t )] = 0 (2.12)
and thus follows the dynamics

d L i (t ) = L i (t )ζ i (t )d W i t (2.13)
where W i t is a Brownian motion, σ i is a bounded deterministic function. The forward measure dynamics, under a forward-adjusted measure Q T i is given by :

d L i (t ) = i ∑ j =g (t ) δ j L j (t )L i (t )ζ i (t )ζ j (t ) 1 + δ j L j (t ) d t + L i (t )ζ i (t )d W i t ,q (2.14)
The interest rate swap value with the strike K is given by :

S b (t i ) = b-1 ∑ k=i P(t i , t k+1 )δ k (L k (t i ) -K) (2.15)
A multifactor BGM ([4], [1]) is given by the following formula :

d L i (t ) = i ∑ j =g (t ) δ j L j (t )L i (t ) ∑ d q=1 ζ i ,q (t )ζ j ,q (t ) 1 + δ j L j (t ) d t + L i (t ) d ∑ q=1 ζ i ,q (t )d W i t ,q (2.16)
where d is the number of factors.

Expected Exposure

We show the exposure profiles computation on the example of a two factor BGM model for LIBOR rates, by simulating exposure profiles of LIBOR swaps with 5 year maturity, which coincides with the last exercise date and lockout after 3 months and tenor spacing of three months.

Using Euler's discretization scheme [29], we have :

L i (t n+1 ) = L i (t n ) + µ i (L(t n ), t n )L i (t n )∆t + L i (t n ) p ∆t 2 ∑ q=1 ζ i ,q (t n )ϵ n+1,q (2.
µ i (L(t n ), t n ) = i ∑ j =g (t n ) δ j L j (t n ) ∑ 2 q=1 ζ i ,q (t n )ζ j ,q (t n ) 1 + δ j L j (t n ) (2.

18)

Then for m = 1,...,M, at time t l the positive (PE) and negative (NE) exposures are given by V (m),+ i

(t l ) = max(S (m) l (t l ), 0) V (m),- i (t l ) = min(S (m) l (t l ), 0) (2.19)
The expected positive exposure(EPE) at each time instant t l :

EPE t l = 1 M M ∑ m=1 V (m),+ l (t l ) (2.20)
In figure 2.1. you can see the exposure profiles of 1000 different scenarios.

Multi-curve Interest Rate Modelling

Before the crisis of 2007-2009 Euribor and Eonia overnight index swap (OIS) rates (a fixed swap rate that is paid/received in exchange for a floating rate tied to the Euro overnight rate Eonia), when both rates were assumed to be risk-free and followed each other closely. But during the crisis the situation radically changed and the spread between Euribor and Eonia OIS rates increased for all tenures. Since Euribor/Libor is not risk-three, many financial institutions moved away from the traditional discounting, instead Eonia is considered now as risk-free rate.

The previous assumption that one curve could be used for discounting and funding is not true anymore :

L i (t ) ̸ = E Q T i [L(t , T i )]
(2.21)

Credit Risk

In this section we develop a L évy driven α-stable stochastic intensity model. For the sake of completeness, we give some notions on L évy processes, for more detailed treatment, we refer to [3], [17] and [14]. In the next sections, we provide some well-known properties and definitions on α-stable processes, then we show sufficient conditions for CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING a L évy process to be decomposed into the sum of a drift, a Brownian motion and a compound Poisson process, that will be very important in imposing correlation structure of CVA with WWR.

In credit risk modelling there are two main approaches to model default probabilities. First one is the first passage approach, which is very similar to the first cross of a given barrier, the technique is widely used in a barrier option pricing and structural models developed by Merton. Second approach is a stochastic intensity model, which will be used for default probabilities computation in the next sections.

The first jump of the process can be simulated, conditional on the filtration F as (1), and by inversion we define the default time as :

Λ(τ) =: ξ = E xpo
τ = Λ(ξ) -1 .
The survival probability is given by :

P(τ > t ) = E[e - ∫ t 0 λ u du ] (2.22)
The density of the first default time is given by :

f τ (t ) = E[λ t e - ∫ t 0 λ u du ] (2.23)
As a simple example consider CIR process that describes the dynamics of the stochastic intensity λ t

d λ t = k(ϕ -λ t )d t + σ λ √ λ t d W t (2.24)
The figure 1.7, shows the path of a CIR stochastic intensity and a survival function.

L évy driven OU Process

Given that X is a L évy process X = {X t : t ≥ 0} in R d , then it is uniquely determined by its characteristic function in the L évy-Khintchine form :

E[e i 〈u,X t 〉 ] = exp ( t [ i 〈µ, u〉 + ∫ (e i 〈u,x〉 -1 -i 〈u, x〉1 ||x||≤1 )ν(d x) ]) (2.25)
Let us consider the following L évy driven OU process :

d λ t = k(ϕ -λ t )d t + ∫ h t (λ t -, x)ϱ(d x,
i : Ω → E, i = 1,...,ϱ(T)}.
If we assume that h t (λ t -, x) = x, then using Ito's formula, we can find the solution of (2.26) :

λ t = λ 0 e -kt + ϕ(1 -e -kt ) + ∫ t 0 ∫ e k(s-t ) h s (λ s-, x)ϱ(d x, d s) (2.27)
Equivalently we can rewrite it as :

λ t = λ 0 e -kt + ϕ(1 -e -kt ) + ϱ(t ) ∑ i =1 h(λ τ-, ξ i ) (2.28) α-Stable Processes Definition X is α-stable if, ∀n ∈ N, i.i.d. r.v. of X, ∑ n i =1 X i = bX + c for some constants b = b(n) > 0 and c = c(n) ∈ R d . Definition A stable X is called α-stable, α ∈]0, 2], whenever X * r d = r 1/α X + c and some constant c = c(r ) ∈ R d (2.29)
Theorem 2.4.1 Let X be a non-trivial L évy process in R with generating triplet (µ, σ, ν).

Then X is α-stable for some α > 0 iff exactly one of these conditions holds :

1. α = 2 and ν = 0.

2. α ∈]0, 2[,µ = 0, and

ν(d x) = (c 1 |x| -1-α 1 x<0 + c 2 x -1-α 1 x>0 )d x on R for some c 1 , c 2 > 0.
(2.30)

Definition A L évy process X with L évy triplet (µ, σ, ν) is said to be tempered α-stable if the L évy measure ν is given by

ν(d x) = ||x|| -1-α υ(x)d x (2.31)

Gaussian Approximation of L évy Driven Process

One of the main challenges to simulate L évy processes roots in the fact that there is no general way to simulate its innovations. In this subsection we present the approach developed in [27] and [28], where necessary and sufficient conditions were established under which the small jumps of a one-dimensional and multidimensional L évi process could be approximated by a Brownian motion. As we mentioned in the first section, this allow us to decompose a L évy process into a Brownian motion and a compound Poisson process. Indeed, we will show that simulated paths of the Levy process resemble a Brownian motion for large values of stability index α and a compound Poisson process for smaller values. Suppose that for ϵ ∈]0, 1] we can decompose the L évy measure ν as :

ν = ν ϵ + ν ϵ (2.32) where ∫ ||x|| 2 ν ϵ (d x) < ∞ and ν ϵ (R d ) < ∞ (2.
X = X 1,ϵ + P ϵ + µ ϵ (2.34)
where X 1,ϵ is a L évy process determined by :

E[e i 〈u,X 1,ϵ t 〉 ] = exp ( t [ i 〈µ, u〉 + ∫ (e i 〈u,x〉 -1 -i 〈u, x〉1 ||x||≤1 )ν ϵ (d x)
]) (2.35)

P ϵ = {P ϵ t : t ≥ 0} is a compound Poisson process with jump measure ν ϵ and µ ϵ = {t µ ϵ : t ≥ 0} is a drift, i.e. µ ϵ = µ + ∫ ||x||≥1 xν ϵ (d x) - ∫ ||x||≤1 xν ϵ (d x) (2.36)
X 1,ϵ is a "small-jump" part of the process X. Next theorem says that, the small-jump part X 1,ϵ can be approximated by p σ ϵ W, where W ∈ R d is a Brownian motion, independent of P ϵ . So, we have :

X D = µ ϵ + p σ ϵ W + P ϵ := X ϵ (2.37)
Next results from [2](Theorem 2.1) gives necessary and sufficient conditions that guarantees Gaussian approximation of small jumps.

Theorem 2.4.2 σ(ϵ) -1 X 1,ϵ D -→ W as ϵ → 0, if and only if for each κ > 0, σ(κσ(ϵ) ∧ ϵ) ∼ σ(ϵ), as ϵ → 0 (2.
lim ϵ→0 σ(ϵ) ϵ = ∞ (2.39)
Example Consider a symmetric α-stable process X that satisfies conditions in theorem 4.1 with L évy measure given by :

ν(d x) = ( c 1 |x| 1+α 1 x<0 + c 2 x 1+α 1 x>0 ) d x (2.40)
where

c 1 + c 2 > 0, then σ(ϵ) 2 = σ 2 ϵ = ∫ |x|≤ϵ x 2 ν(d x) = c 1 + c 2 2 -α ϵ 2-α (2.41)
Now we can check the condition (2.39) :

σ(ϵ) ϵ = √ c 1 + c 2 2 -α ϵ -α/2 ---→ ϵ→0 ∞ (2.42)
So, this condition is satisfied. Example Another example, that we apply in the next sections, is exponentially tempered α-stable process with L évy triplet (0, 0,ν), where L évy measure ν is given by :

ν(d x) = ( e |-x| |x| 1+α 1 ]-∞,0] + e -x x 1+α 1 [0,∞[ ) d x (2.43) µ ϵ = 0 - ∫ ϵ<|x|<1 xν(d x) (2.44)
The variance of small jumps σ(ϵ

) 2 = σ 2 ϵ . σ(ϵ) 2 = σ 2 ϵ = ∫ |x|<ϵ x 2 ν(d x) = ∫ 0 -ϵ x 2 e x (-x) 1+α d x + ∫ ϵ 0 x 2 e -x x 1+α d x = 2 ∫ ϵ 0 x 1-α e -x d x (2.45)
Next, using Karamata's integral theorem, which says that the integrals of regularly varying functions are again regularly varying [31]. For α > -1,

∫ ϵ 0 x 1-α e -x d x ∼ ϵ 2-α 2 -α for ϵ → ∞ (2.46)
So, we have

σ(ϵ) ∼ √ 2 2 -α ϵ 1-α/2
(2.47) and

σ(ϵ) ϵ = √ 2 2 -α ϵ -1/2-α/2 ---→ ϵ→0 ∞ (2.48)
Consequently, condition (2.39) is satisfied.

Jump-Diffusion Approximation of OU L évy Process

A L évy process decomposition allows us to approximate L évy driven SDE by a jump-diffusion SDE. In this section we get "almost exact" solution that allows to compare errors of different numerical schemes under Monte Carlo and IPS approximations. Consider the following OU SDE :

d λ t = µ(λ t )d t + h(λ t -)d X t (2.49)
where X is a L évy process with a triplet (0, 0,ν). Let µ(λ t ) = k(ϕ-λ t ), σ(λ s ) = 0 and h(λ s-) = 1, then the solution of (2.49) is

λ t = λ 0 e -kt + ϕ(1 -e -kt )d t + ∫ t 0 e k(s-t ) d X s (2.50) X t = ∫ t 0 ∫ ||x||≤1 xe ϱ(d x, d s) + ∫ t 0 ∫ ||x||≥1 xϱ(d x, d s) (2.51)
The approximation of a L évy process X t is given by

X ϵ t = ∫ t 0 ∫ ||x||≤ϵ σ ϵ e ϱ ϵ (d x, d s) + ∫ t 0 ∫ ||x||≥ϵ xϱ ϵ (d x, d s) (2.52)
and the approximated solution of eq. (2.49) is

λ ϵ t = λ 0 e -kt + ϕ(1 -e -kt )d t + ∫ t 0 e k(s-t ) d X ϵ s (2.
λ ϵ t = λ 0 e -kt + ( ϕ(1 -e -kt ) + µ ϵ ) d t + ∫ t 0 e k(s-t ) p σ ϵ d W s + ∫ t 0 ∫ ||x||≥ϵ e k(s-t ) xϱ ϵ (d x, d s) (2.54)
The drift of the L évy approximation µ ϵ = 0, i.e.

µ ϵ = 0 - ∫ ϵ<|x|≤1 xν(d x) = - ∫ -ϵ -1 x e x (-x) 1+α d x + ∫ 1 ϵ x e -x
x 1+α d x = 0

(2.55)

Simulation of α-Stable Processes and Numerical Schemes

In previous sections we saw that any L évy process X t with a generating triplet (µ, 0,ν) can be decomposed into the sum of two independent Levy processes :

X t = X 1,ϵ t + P ϵ t (2.56)
The Gaussian approximation of "small-jump" part X 1,ϵ t allows us to get the following approximation of L évy process in equation (2.56) :

X t D = P ϵ t + σ ϵ W t (2.57)
To get the best approximation of X, we need to choose ϵ experimentally and find an efficient way for generating a compound Poisson process P ϵ as ϵ → 0. That can be done using a series representation. Rosi ńsky showed that L évy process X t with generating triplet (µ, 0,ν) can be represented almost surely as a uniform in t convergent series of the following form :

X t = ∞ ∑ j =1 [ H(Γ j , ξ j )1 u j ≤t -t µ j ] , 0 ≤ t ≤ T (2.58)
where {Γ j } are arrival times in a Poisson process with rate one, {ξ j } are i.i.d. jump marks taking values in some Euclidean space, {u j } are uniform i.i.d. r.v. on the interval [0, T] and the sequences {Γ j }, {ξ j } and {u j } are independent from each other. The function H is a jointly measurable real-valued function, such that the map r 7 → |H(r, v)| is nonincreasing and µ j ∈ R. According to [27], the choice of H and ξ is not unique, but they have to satisfy to condition :

∫ R f (x)ν(d x) = T -1 ∫ ∞ 0 E[ f (H(r, ξ))]d r (2.59)
for any nonnegative Borel function f with f (0) = 0. For the case of P ϵ we have the following series representation :

P ϵ,τ t = ∑ j =Γ j ≤τ [ H(Γ j , ξ j )1 u j ≤t -t μ(τ) ] , 0 ≤ t ≤ T (2.60)
where

μ = ∫ s 0 ∫ ||x||≤1 xϑ(r, d x)d r, s ≥ 0 (2.61)
for B ∈ B(R), the measure ϑ is defined by : Series representation of a L évy processes provides a uniform approximation of a L évy process along the sample paths. Next step is the simulation of the L évy process based on this representation. A detailed treatment of stable processes can be found in [28], we use a series representation in the case of tempered-α stable processes. Let Γ j , U j , Υ j , Ξ j and ξ j be sequences of i.i.d. random variables, such that P(ξ j = ±1) = 0.5, Υ j are exponentially distributed with rate one, Ξ j are uniformly distributed on ]0, 1[, then we have

ϑ(r, B) = P(H(r, ξ i ) ∈ B), r > 0, and ν(B) = ∫ ∞ 0 ϑ(r, B)d
P ϵ t = ∞ ∑ j =1 ξ j ( [ αΓ j 2κT ] -1 α ∧ Υ j Ξ 1 α j ) 1 U j ≤t , 0 ≤ t ≤ T (2.63)
For the compound Poisson process P ϵ,τ , we simulate it using the following series :

P ϵ,τ t = ∑ j =Γ j ≤τ ξ j ( [ αΓ j 2κT ] -1 α ∧ Υ j Ξ 1 α j ) 1 U j ≤t , 0 ≤ t ≤ T (2.64)
Details are given in the Algorithm 6.

Algorithm 6 : α-stable process generation using series representation

αι(n) 2κT ] -1 α ∧ Υ n Ξ 1 α n ) ; 22 end 23 end

Numerical Schemes

In this section we show the way SDE in eq.( 2.49) could be solved numerically using three discretization schemes : Euler with a constant time step, Euler and Milstein with

d λ t = k(ϕ -λ t )d t + σ ϵ d W t + ∫ ||x||≥ϵ xϱ ϵ (d x, d s) (2.65)
The simplest way to discretize it is the Euler scheme.

λ t n+1 = λ t n + k(ϕ -λ t n )∆t + σ ϵ ∆W n + ∫ t n+1 t n ∫ |x|>ϵ xϱ(d x, d s) = = λ t n + k(ϕ -λ t n )∆t + σ ϵ ∆W n + ϱ(t n+1 ) ∑ k=ϱ(t n )+1 ξ k (2.66)
where ∆t = t n+1 -t n and ∆W n = W t n+1 -W t n ∼ N (0, ∆t) is n-th increment of the Brownian motion process.

A jump-adapted Euler scheme was introduced in [21]. Define λ t n+1 -= lim s↑t n+1 λ s as the left limit of the process λ at time t n+1 .

Consider the following equidistant discretization of the interval [0, T] : 0 = t 0 , ..., < t n = T and let the jump times to be denoted as τ 1 , ..., τ ϱ(T) , then the jump-adapted discretization be {t 0 , t 1 , ..., t ϱ(T) } = {t 0 , t 1 , ..., t N } ∪ {τ 0 , τ 1 , ..., τ ϱ(T) )}.

The process before jump is given by

λ t n+1 -= λ t n + µ∆ t n + σ∆W t n (2.67)
and the process after the jump

λ t n+1 = λ t n+1 -+ ∫ h(λ t n+1 -, x)ϱ(d x, {t n+1 }) (2.68)
Note that the diffusion approximation is given by (2.67), and the impact of jump is simulated in the eq. (2.68). We can represent the two equations in a simple form :

{ λ t n+1 = λ t n+1 -+ h(λ t n+1 -, ξ ϱ(t n+1 ) ) if jumps at t n+1 λ t n+1 = λ t n+1 - otherwise (2.69) where ∫ h(λ t n+1 -, x)ϱ(d x, {t n+1 }) = 0, if t n+1 is not a jump time.
We can achieve a strong order -1 convergence, if we use a Milstein scheme. The process after the jump has the same form as in (2.68).

λ t n+1 -= λ t n + µ∆ t n + σ∆W t n + σ 2 2 ( (∆W t n ) 2 -∆t n ) ( 2

WWR Modelling

One of the most interesting cases of CVA valuation is a dependence between the EPE and the default probability of a counterparty. When the dependence is negative, i.e. exposure tends to increase, when the probability of default of our counterparty increases or the credit quality worsens, we have the WWR. When the dependence is positive, we have Right Way Risk. There are a number of different approaches to correlate the exposure and default probability. One can use copulas, or correlate two Brownian motions, i.e. 〈d W t , d W λ t 〉 = ρd t , where ρ is a correlation parameter. Another approach proposed in ( [10], [23], [24] and [25]), models WWR by adding jumps in the exposure process (jump-at-default) modelling to reflect the effect of drastic drops in credit ratings to the underlying risk factors.

The following proposition, inspired by [30], allows to incorporate WWR into CVA.

Proposition 2.5.1 Given the expected positive exposure V + τ , the default time τ, the CVA with wrong-way risk is given by

CVA WWR = LGDE[B τ V + τ 1 τ<T ] = LGD ∫ T 0 E[B t V + t f τ (t )]d t (2.71)
Proof By using the tower property of the conditional expectation and Fubini theorem, we have the following :

LGDE

[B τ V + τ 1 τ<T ] = LGDE[E[B τ V + τ 1 τ<T |G t ]] = LGDE [∫ T 0 E[B τ V + τ |G t , τ = t ] f τ (t )d t ] = (2.72) = LGDE [∫ T 0 B t V + t f τ (t )d t ] = LGD ∫ T 0 E[B t V + t f τ (t )]d t (2.73)
As we mentioned in previous sections, α stable process X t allows Gaussian approximation of small jumps, so that we can express it as :

X t D = σ ϵ W t + P ϵ t (2.

Particle Interpretations

Interacting Particle System for Rare Events Estimation

Interacting particle system is efficient method to calculate rare-event probabilities. For example, usually we formulate the problem in the following way. Given a sequence of rare event sets A n ⊂ ... ⊂ A 1 , we would like to compute the following expectation after some Radon-Nikodym transformation, that turns a rare event into a regular one :

E P [g (λ N ) N ∏ n=1 1 λ n ∈A n ] = E Q [g (λ N ) d P d Q N ∏ n=1 1 λ n ∈A n ] (2.76)
In practice we do not have an explicit form of the likelihood ratio d P d Q , so it becomes unfeasible unless one considers some simple examples. One of the solutions was interacting particle system(IPS), which was proposed by Del Moral and Garnier [8], where they CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING proposed to generate particles(samples) in two steps, i.e. particle mutation and selection. The idea is to approximate the ratio of P with respect to some importance measure Q by choosing a weight function ω that approximates the Radon-Nikodym derivative dP d Q . If we assume that P and Q have density function p and q respectively, for M particles {λ (m) n } M m=1 we can define the weight function iteratively by :

N ∏ n=1 ω n (λ (m) n ) ∝ d P d Q (λ (m) N ) = N ∏ n=1 p n (λ (m) n , λ (m) n-1 ) q n (λ (m) n , λ (m) n-1 )
(2.77)

Since two measures P and Q form a Markov chain, the Radon-Nykodim derivative can be decomposed into the product of ratio of the transition density p n (•, λ (m) n-1 ) to the transition density q n (•, λ (m) n-1 )with respect to the measure Q. The normalized importance weight function is given by :

W n (λ (m) n ) = ω n (λ (m) n ) 1 M ∑ M j =1 ω n (λ ( j ) n ) (2.78)
The IPS estimate of an expectation (2.76) has the following form :

C IPS = E [ f (λ N ) ∏ N-1 n=1 ω n (λ n )1 λ n ∈A n ] E [ ∏ N-1 n=1 ω n (λ n )1 λ n ∈A n ] (2.79)
In our experiments we choose a potential function (an unnormalized importance weight) of the form :

ω n (λ (m) N ) = N ∏ n=1 e δ(λ (m) n -λ (m) n-1 ) (2.80)
where δ is an exponential tilting parameter.

Mean Field IPS Approximations

In this subsection we give a brief overview of interacting particle methods, that is very important to understand the mechanism behind the sequential importance resampling scheme that we will use in the next sections. For a thorough survey we refer to [5], [6].

We can express eq.(2.76) through the linear semigroups associated to the sequence of non-negative measures γ n for any test function g that have the following form :

C IPS = γ N (g ) = η N (g )γ N (1)
(2.81)

η n represent normalized measures. The flow of Feynman-Kac measures (γ N , η N ) for any measurable and bounded function g is given by

γ N (g ) = E[g (λ N ) N-1 ∏ n=0 ω n (λ n )] and η n (g ) = γ n (g ) γ n (1) (2.82)
where the sequence of measures γ n satisfy a linear recursion

γ n (g ) = γ n-1 (Q n (g )) (2.83) with integral operator Q n : Q n (g )(x) = E[g (λ n )|λ n-1 = x] = ∫ Q n (x, d y)g (y) (2.
Q n (x, d y) = ω n-1 (x)M n (x, d y) (2.85)
The sequence of measures η n satisfy the nonlinear equation :

η n = Φ n (η n ) = Ψ n-1 (η n-1 )M n (2.86)
where Φ n : E n → E n+1 and the operator

Ψ n : E n → E n is Ψ n (η n )(x) = ω n (x) η n (ω n ) η n (x) (2.87)
Observe, that eq.( 2.86) represents very important class of McKean models. The main problem of Monte Carlo simulation of this type of models stems in the fact that it can not be computed exactly, since its Markov transitions M depends on transitions that are unknown( [7], [5]). The IPS allows to approximate the solution by a discrete-time process.

Assume that the Markov transition kernels M n have a sequence of densities M n , and there exist some reference probability measures κ n , such that M n ≪ κ n , i.e. M n are absolutely continuous with respect to κ n :

M n (x, d y) = M n (x, y)κ n (d y) (2.88) Radon-Nikodym derivative given that Q n (x, •) ≪ κ n : d Q n (x, •) d κ n (y) = ω n-1 (x)M n (x, y) (2.89)
The following result shows the approximation of importance measure Q n by a flow of non-negative measures η n . The lemma is using techniques, that were adapted from [5] and [6].

Lemma 2.6.1 Assume that sup x,y∈E n M n (x, z) M n (y, z) := h n (z) s.th. ||h n || < ∞ and sup x,y∈E n ω n (x) ω n (y) < ∞, ∀n ≥ 0 (2.90)
then the integral operator (2.84) is given by

Q n (g )(x) = η n-1 (ω n-1 ) ∫ ω n-1 (x)M n (x, y) η n-1 (ω n-1 M n (•, y)) g (y)η n (d y) (2.

91)

Proof Using eq.( 2.86) we have :

η n (d y) = Ψ n-1 (η n-1 )M n (d y) = η n-1 (ω n-1 M n (•, y))κ n (d y) η n-1 (ω n-1 ) (2.92)
So the importance density takes the following form :

Q n (x, d y) = η n-1 (ω n-1 ) ω n-1 (x)M n (x, y) η n-1 (ω n-1 M n (•, y)) η n (d y) (2.93)
The statement is obtained by applying eq.(2.84).

CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING IPS could be understood as a stochastic linearization of the nonlinear recursive equation(2.86). The algorithm involves construction of E M

n valued M-tuple Markov Chain X n = {X (i ) n } M i =1 ∈ E M n with transitions given by P(X n ∈ d (x (1) , ...,

x (M) )|X n-1 ) = M ∏ i =1 Φ n (η M n )(d x (i ) ) (2.94)
to approximate measure η n :

lim n→∞ η M n (g ) = η n (g ) where η M n = 1 M M ∑ i =1 δ X (i ) (2.95)
The above could be proven using the law of large numbers under certain regularity conditions. Now we can obtain the approximation b Q n of the importance measure Q n by replacing the measures η n by their M-tuple approximation η M n , which is the definition of mean field IPS approximation scheme :

b Q n (g )(x) = η M n-1 (ω n-1 ) ∫ ω n-1 (x)M n (x, y) ∫ η M n-1 (d z)(ω n-1 (z)M n (z, y) g (y)η M n (d y) = = η M n-1 (ω n-1 ) M ∑ i =1 ω n-1 (x)M n (x, X (i ) ) ∑ M j =1 η M n-1 (d z)(ω n-1 (z)M n (z, X ( j ) )
g (X (i ) ) (2.96)

Particle Monte Carlo Convergence Properties

There is a lot of research on the convergence of particle methods. Now, given M particles {X (i ) } M i =1 , we cite from [8] the following result :

1 M M ∑ i =1 g (x)δ X (i ) N (x) ----→ M→∞ E[g (λ N )] a.s. (2.97)
where δ is a Dirac function. The following proposition shows that the random sequence {γ N n (1)} k n=1 converges in probability to a deterministic sequence {γ n (1)} k n=1 as its number M goes to infinity. Proposition 2.6.2 The Particle estimator is unbiased estimator, so that for any p ≥ 1 and for all bounded and measurable functions g on some measurable space

(E k , E k ) with ||g || ≤ 1, we have E [ γ N k (g ) ] = γ k (g ) (2.
Λ(t N ) = ∫ t N 0 λ u d u
for the time interval 0 = t 0 < ...t N = t . From large deviation analysis, we know that P(ξ > Λ(t N )) → 0, as Λ(t N ) → ∞, i.e.

P(τ > t

N ) = P(ξ > Λ(t N )) = E P [1 ξ>Λ(t N ) ] = E P [1 ξ>Λ(t ) e αξ-αξ ] ≤ e -(αΛ(t N )-Γ(α))
with Γ(α) = log(E P [e αξ ]) (2.103)

Unnormalized particle approximations are given by :

E P [1 ξ>Λ(t N ) ] = E Q [1 ξ>Λ(t N ) e -(αΛ(t N )) ] = E Q [1 ξ>Λ(t N ) N ∏ l =1
e α(Λ(t l )-Λ(t l -1 )) ] (2.104) Now we can give particle interpretation of the default probability. Assuming that initial values are equal, i.e. p 0 = µ 0 and for any bounded measurable function g : E → R, under new measure Q we have the following representation of the original expectation under measure P :

f τ (t N ) = E Q [λ t N e - ∫ t N 0 λ u d u ∏ N-1 l =1 ωl (λ l )] E Q [ ∏ N-1 l =1 ω l (λ l )]
(2.105)

In figure 5, we can see the simulated paths of three discretization schemes with constant and jump-adapted time steps to compare with almost-exact solution.

Numerical Analysis and Discussion

Numerical Schemes Error Estimates

We split errors into two categories : ||•|| L 2 and ||•|| L 1 errors. Strong errors are given in terms of || • || L 2 norm. Given the estimates b X l l t , where l l = {MC Eu , MC Eu-j ad p , MC Ml s-j ad p } set of numerical schemes, of the process X t , we define a strong error estimates as :

er r s = ( E [ | b X t -X t | 2 ]) 1/2 ≈ v u u t 1 M M ∑ i =1 | b X (i ) t -X (i ) t | 2 (2.106) || • || L 1 errors are defined as : er r w = E [ | b X t -X t | ] ≈ 1 M M ∑ i =1 | b X (i ) t -X (i ) t | (2.107)
Table 1 and 2 shows ||•|| L 2 and ||•|| L 1 errors respectively of 5000 simulated Monte Carlo paths and 1000 particles. Level of stability is chosen to be α = 0.8, when most of jumps are large. We observe that errors are sufficiently minimized when we use IPS. The same demonstrates figure 7. Figure 6 shows that the particles variability and errors are converging between -7 and -8 for three numerical schemes as we increase the number of time steps.

CVA Computation

In this subsection we present the computation methodology for CVA estimation. We consider two cases of CVA modelling :

-CVA i nd , when there is no correlation between default probability and the exposure process. -CVA w wr , when we assume a negative correlation between market risk and credit risk. First we outline two factor BGM model that will be used for expected exposure V t simulation. Second, depending on the case of presence or absence of the WWR, we simulate the Levy driven OU stochastic intensity process.

Exposure Profiles Specification

We consider a two-factor BGM to model LIBOR interest rates :

d L i (t ) = i ∑ j =g (t ) δ j L j (t )L i (t ) ∑ 2 q=1 ζ i ,q (t )ζ j ,q (t ) 1 + δ j L j (t ) d t + L n (t ) 2 ∑ q=1 ζ i ,q (t )d W t (2.108)
Euler discretization scheme gives us :

L i (t n+1 ) = L i (t n ) + µ n (L i (t n ), t n )L i (t n )∆t + L i (t n ) p ∆tσ i (t n )ϵ n+1 (2.109)
where ∆t = t n+1t m and drift µ n is given by

µ i (L i (t n ), t n ) = i ∑ j =g (t ) δ j L j (t )L i (t ) ∑ 2 q=1 ζ i ,q (t )ζ j ,q (t ) 1 + δ j L j (t ) (2
i nd = LGD ∫ T 0 E[B t V + t ] f (t )d t = LGD N ∑ l =1 ∫ T l +1 T l E[B t V + t ] f (t )d t = (2.111) = LGD N ∑ l =1 B l M ∑ m=1 V (m),+ t l ∫ T l +1 T l f τ (t )d t = (2.112) = LGD N ∑ i =1 B l M ∑ m=1 V (m),+ t l M j ∑ j =1 λ ( j ) l e - ∑ l k=1 λ ( j ) k ∆k ∆t (2.113)
where ∆t = T l +1 -T l , B l = B(0, T l +1 ) and ∆k = t kt k-1 .

In the case of CVA with WWR we specify different levels of correlation ρ between two Brownian motions in (2.75). Monte Carlo estimator of CVA with WWR is

CVA WWR = LGD ∫ T 0 B t E[V + t f τ (t )]d t = LGD N ∑ l =1 ∫ T l +1 T l B t E[B t V + t f τ (t )]d t = (2.114) = LGD N ∑ l =1 ∫ T l +1 T l B t M ∑ m=1 V (m),+ t l f (m) τ (t )

CVA. Particle Monte Carlo

Another unbiased estimator CVA P of the expectation in eq.(2.6) can be obtained by means of IPS.

The exposure profiles are simulated similarly as in the previous section. The simulation of a survival probability, as it is described in algorithm 2, consists of two steps : sampling (mutation) and survival (selection). Define M initial particles {X (m) } M m=1 . The selection step consists of randomly choosing a set of M particles from the selection transition, which, for example, could be chosen as :

S n = ω n (X (m) n )1 X (m) n (x) + (1 -ω n (X (m) n )) ∑ i ω n (X (m) n ) ∑ j ω n (X ( j ) n ) 1 X (m) n (x) (2.116)
Define selected particles as λ n . The mutation step from λ n to λ n+1 : every survived particle λ (i ) n propagates randomly to a new particle λ (i ) n+1 = y randomly chosen from the transition kernel M (λ (i ) n , d y). At the final time N, the CVA estimate is

CVA P,i nd = LGD N ∑ l =1 1 M M ∑ i =1 B l V (m),+ l × [( N-1 ∏ l =1 1 M M ∑ m=1 ω l (λ (m) l ) ) × M ∑ m=1 f τ (λ (m) N ) ] (2.117)
CVA with WWR estimate is given by CVA

P,WWR = LGD N ∑ l =1 1 M M ∑ i =1 B l V (m),+ l f τ (λ (m) l ) × ( N-1 ∏ l =1 1 M M ∑ m=1 ω l (λ (m) l ) ) (2.118)
We assume that parameters of the underlying stochastic processes are already calibrated. We simulate 20 times CVA values :

C VA = M s ∑ l =1 CVA l (2.119)
where M s is the number of repeated simulations.

Numerical Results

We simulate CVA values using two numerical methods : first, EPE using standard MC approach, second, EPE using particle methods ; We compare CVA without WWR and CVA with WWR on the one hand, and compare with Basel III recommendations on the other.

Table 5. shows that Basel III's approximation is overvalued we use the formula (2.118) for 10 4 and 10 5 scenarios.

Next, we are interested in CVA values with WWR, when we correlate two Brownian motions, the one is the exposure and another one in Gaussian small jump approximation of Levy process. Table 6. shows CVA with RWR values for different levels initial intensity and correlation.

Conclusion

In this paper we adapted IPS approach to the simulation of α-stable process. Using mean field algorithms we reduced ||•|| L 2 and ||•|| L 1 errors for three discretization scheme : Euler with constant time step, Euler and Milstein with jump-adapted time steps compared to a vanilla Monte Carlo algorithm. IPS allowed reduce the number of simulation paths and reduce the variability of estimates.

We also developed α-stable intensity process to estimate probability of default. We showed how WWR could be naturally embedded into the CVA values. Importantly, our results demonstrated that the Basel III approximation of CVA with WWR far from being adequate. As an extension, in the next article we will use OIS spreads in the exposure valuation in the case of other products, such as Bermudan swaptions.

Chapitre 3 Hamiltonian flow Simulation of Rare Events

« "Die Welt ist alles, was der Fall ist" » Ludwig Wittgenstein « "The world is all that is the case" » Ludwig Wittgenstein

Introduction

Hamiltonian flow based Monte Carlo simulations originates from physics have been used in many applications in statistics, engineering for a number of years. However these methods are not widely used in the estimation of rare events and in the financial option pricing practise.

This paper proposes Hamiltonian flow Monte Carlo technique for an efficient estimation of the rare event probability. Similarly to an importance sampling technique this method involves a change of probability measure. The random variables are sampled according to a modified probability measure that differs from a reference measure.

A rare event is the probability P( f (x t ) > a t ) for large values of x. One way to deal with this problem is to change an original measure, so that a k is not too large in the new measure. Define the following set :

A t = {x ∈ R M , f (x t ) > a t } (3.1) E P [1 A t ] = E Q [1 A t L t ] (3.2)
where L = dP dQ is a Radon-Nykodim derivative. Hamiltonian approach in the Monte Carlo context was developed by Duane et al. [11], R. Neal [15] where they proposed an algorithm for sampling probability distributions with continuous state spaces. The advantage of Hamiltonian based Monte Carlo is in the fact that we can extend the state space by including a momentum variable that will force in our context to move long distances in the state space in a single update step. We use this property of Hamiltonian dynamics to explore rarely-visited areas of the state space and efficiently estimate rare-event probability. Algorithm consists of two parts : simulation of Hamiltonian dynamics and Metropolis-Hastings test, that removes the bias and allows large jumps in the state space. We will show the performance of the algorithm on the Down-Out Barrier option technique with low level of a barrier.

HFMC can be considered within the optimal transportation problem, which was posed back in the 18-th century. Like HFMC, other simulation based approaches such as particle methods [5], [8] or the transportation using a homotopy [10] allow to move a set of particles from the measure P to the measure Q, by minimizing the transportation cost. We will show how the rare events estimation can be computed using interacting particle systems [6].

The paper is organized as follows. Section 2 introduces state of the art approach to estimate rare events : MC and PMC and formulates the problem within the context of Barrier option pricing. Section 3. describes Hamiltonian flow Monte Carlo algorithm. Section 4 presents the law of large numbers and convergence for HFMC. An adaptation of Hamiltonian flow Algorithm in the case of a Barrier option, numerical results and discussion is presented in Section 5. Section 6 concludes.

Monte Carlo and Interacting Particle System

Problem Formulation

Barrier option pricing is one of the cases when we encounter to the case of rare events. Consider a sequence of random variables {X n } n≥0 , which in the financial context can be interpreted as asset prices, which forms Markov Chain on the space R n S . Given some stochastic process {X t } t ≥0 , for any test function f , we would like to compute the following expectation : CHAPITRE 3. HAMILTONIAN FLOW SIMULATION OF RARE EVENTS

C = E P [ f (X T )1 X t ∈[0,T] ∈A t ] (3.3)
One of the most popular ways to deal with this problem is importance sampling, when we replace the original statistical measure P by an importance measure P δ . Then for X n = (X 1 , ..., X n ). :

C = ∫ f (X n t ) d P d P δ (X n t )d P δ (X n t ) n t ∏ n=1 1 X n ∈A n = E P δ [ f (X n t ) d P d P δ n t ∏ n=1 1 X n ∈A n ] (3.4)
where the rare event set A n is given by :

A n = {X n ∈ R n S , f (X n ) > a n } (3.5)
In practise we don't have an explicit form of the likelihood ratio d P dP δ , so it becomes unfeasable unless one considers very simple toy examples. One of the solutions was interacting particle system(IPS), which was proposed by Del Moral and Garnier [8], where they proposed to generate particles(samples) in two steps, i.e. particle mutation and selection. The idea is to approximate the ratio of P with respect to some importance measure P δ by choosing a weight function ω n that approximates the Radon-Nikodym derivative d P dP δ . If we assume that P and P δ have density function p and e p respectively, for n S particles {X (m) n } n S m=1 we can define the weight function iteratively by :

n t ∏ n=1 ω n (X (m) n ) ∝ d P d P δ (X (m) n t ) = n t ∏ n=1 p n (X (m) n , X (m) n-1 ) p δ n (X (m) n , X (m) n-1 ) (3.6)
Since two measures P and P δ form Markov chain, the Radon-Nykodim derivative can be decomposed into the product of ratio of the transition density p n (•, X (m) n-1 ) to the transition density with respect to the measure P δ .

The normalized importance weight function is given by :

W n (X (m) n ) = ω n (X (m) n ) 1 n S ∑ n S j =1 ω n (X ( j ) n ) (3.7)
The IPS estimate of an expectation (3.3) will have the following form :

b C IPS = E [ f (X n t ) ∏ n t -1 n=1 ω n (X n )1 X n ∈A n ] E [ ∏ n t -1 n=1 ω n (X n )1 X n ∈A n ] (3.8)
In our experiments we choose a potential function(an unnormalized importance weight) of the form :

ω n (X (m) n t ) = n t ∏ n=1 e δ(X (m) n -X (m) n-1 ) (3.9)
where δ is an exponential tilting parameter. One issue with this approach is in the fact that an optimal choice of tilting parameter δ has to be judiciously chosen by running simulations and, in fact, it is fixed across all time steps n = 1,...,n t .

To approximate Radon-Nikodym derivative dP dQ we will generate Markov Chain that will converge to an ergodic distribtion. Let us introduce a random process X u in a pseudotime u and consider the following SDE, which is a gradient flow distrubed by a noise :

d X u = -∇Ψ(X u )d u + 2 √ β -1 d W u (3.13)
where Ψ(X) := -log(p(X)) is a potential. Under the assumption of ergodicity, the autocorrelated path X u asymptotically, i.e. u → ∞ draws samples from a stationary distribution :

π(X) = 1 Z exp(Ψ(X)) (3.14)
where a normalizing constant Z :

Z = ∫ R n e -βΨ(x) d x (3.15)
This can be seen as a unique solution of the following Fokker-Plank equation, given that Ψ satisfies to some growth condition :

∂p(t, x) ∂t = d i v(∇(Ψ(x)ρ)) + β -1 ∆p (3.16)
When we mentioned ergodicity, we meant, that for a class of regular functions ϕ : R X → R and x 0 a.s., the Markov Chain satisfies :

1 L L ∑ l =1 ϕ(x l ) → ∫ R X ϕ(x)π(d x) = E π [ϕ(X)] (3.17)
Observe that eq. (3.13) is a reversible process, which is interesting from theoretical point of view, but in practise the speed of convergence is not optimal. One of the ways to improve the convergence is to add a divergence-free drift b and consider the following modified SDE :

d X u = (-∇Ψ(X u ) + b(X t ))d u + 2 √ β -1 d W u (3.18)
in order to satisfy detailed balance condition, we assume that ∇(be -Φ ) = 0.

Another way to improve the convergence is to consider a generalized Langevin SDE :

Ẍγ u = -∇Ψ(X γ u ) - Ẋγ u + √ 2β -1 Ẇu (3.19)
We can rewrite it as :

{ d X u = P u d u d P u = -∇Ψ(X u )d u -P u d u + √ 2 β d W u (3.20)
where the pair (X, P) is a kinetic process with X is the position and P = d X du is the velocity, that acts as an instantaneous memory.

The invariant function of the Markov process {x, P}, if it exists, is given by :

π 0 (x, P) = 1 Z e -βH (x,P) , Z = ∫ R 2 e -βH (x,P) d Pd x (3.21)
where

H (x, P) = 1 2 PM -1 P + Ψ(x) (3.22)
equal to one, this is because candidates are proposed though simulation of Hamiltonian flow. Generate X (m) n from prior e p(X (m) 0 , •) ;

5

Simulate initial momentum P (m) 1 ∼ N (0, I M ), set x (m) H = X (m) n and run Hamiltonian flow :

6 for l f = 1,...L -1 do 7 P (m) H ((l f + 1 2 )δ) = P (m) H (l f ) -δ 2 ∂Ψ ∂x H (x (m) H (l f )) x (m) H ((l f + 1)δ) = x (m) H (l f ) + δP (m) H ((l f + 1 2 )δ)I -1 M P (m) H ((l f + 1)δ)) = P (m) H ((l f + 1 2 )δ) -δ 2 ∂Ψ ∂x H (x (m) H ((l f + 1)δ)) 8 end 9
Calculate acceptance probability and set x * = x (m) H (L), P * = P (m) H (L) :

a = 1 ∧ e (-H (x * ,P * )+H (x (m) H ,P (m) H ))∆t ; (3.45) Draw u ∼ U nif(0, 1) ; 10 if u < a then 11 Set X (m) n+1 = x * ;
C HFMC = 1 n S n S ∑ m=1 ( f (X (m) n t ) n t ∏ n=1 e (-H (X (m) n+1 ,P (m) n+1 )+H (X (m) n ,P (m) n )∆t 1 X (m) n ,X (m) n+1 ∈A n )

Convergence Analysis

IPS convergence

IPS convergence, and in particular the asymptotic behaviour as number of particles n S → ∞ was thoroughly studied in [6].

The following result given in [5] allows a non asymptotic control of variance of the rare event probability.

E [ ¯CIPC C -1 ¯2] ≤ 4 n S n t ∑ s=1 b δ (n t ) s p k (3.47)
where b δ

(n t ) s = ∏ s≤k<s+n t e δ k

HFMC Convergence

LLN and Convergence Rate

Birkhoff ergodic theorem allows us have law of large numbers(LLN) like convergence. So, we are interested in a sigma-algebra G of invariant events, in particular when G is trivial.

From lemma 3.3 we know that Markov chain generated by HFMC is irreducible, and we can see that Markov Chain that we get from the rare event kernel M satisfies irreducibility conditions due to the fact that the transition density is always positive. Applying the results by [17], we have : Proposition 3.4.3 [17] Suppose that M is a π-irreducible Metropolis kernel. Then M is a Harris reccurent. Proposition 3.4.4 [17] If M is positive Harris and aperiodic then for every initial distribution λ :

|| ∫ λ(d x)(M ) l (x, •) -π|| TV → 0, l → ∞ (3.48) for π almost all x.
where || • || TV is a total variation distance.

Geometric Ergodicity

To establish central limit theorem (CLT), we need a geometric ergodicity of the chain.

Definition A subset C of that state space (R n S , B(R n S )) is petite if there exists a non-zero positive measure ν on the state space and subsampling distribution q such that

K q (x, A) ≥ ν(A), ∀A ∈ B(R n S ) and x ∈ C (3.49) Definition A subset C of that state space (R n S , B(R n S )
) is small if there exists a non-zero positive measure ν on the state space and real-valued number l ∈ R such that

K l (x, A) ≥ ν(A), ∀A ∈ B(R n S ) and x ∈ C (3.50)
Observe, that every small set is petite.

Theorem 3.4.5 Suppose for an irreducible, aperiodic Markov chain having transition probability kernel K and a state space R n S , there exists a petite set C areal valued function V, satisfying v ≥ 1, and constants b < ∞ and λ < 1 such that

K V(x) ≤ λV(x) + b1 C (x), ∀x ∈ R n S (3.51)
holds. Then the chain is geometrically ergodic.

The function V is called a geometric drift. Take the expectation of the both sides of (3.51) and using the invariance of measure π with respect to the kernel K :

E π [V(X)] ≤ bπ(C) (1 -λ) (3.52)
In other words, for λ ∈ (0, 1] a function satisfying (3.51) is always π-integrable.

Proposition 3.4.6 Assume that there exist λ ∈ [0, 1) and b ∈ R + such that

K V ≤ λV + b (3.53)
and

lim sup K (x, R(x) ∩ B(x)) = 0 (3.54)
It was shown in [12] that under certain conditions, HFMC kernel is geometrically ergodic.

Applications and Numerical Results

We will test our algorithm on down-out(DOC) Barrier option pricing, and compare its estimate with a standard Monte Carlo and particle Monte Carlo methods. Lets consider a toy example and assume that our asset follows the following SDE :

d X t = µX t d t + σX t d W t (3.55)
where µ is a drift, σ is a constant volatility parameter. European DOC call Barrier option is a usual call option contract that pays a payoff max(S T -K, 0), provided that the asset price S has not fallen below a barrier B during the lifetime of the option. If the pricing process ever reaches the barrier B, then the option becomes worthless. We use Euler-Muruyama disretization scheme and we use the following notation X t n := X n , so for a time discretization : 0 = t 0 , t 1 , ..., t n t = T, the solution of the SDE in (3.55) :

X n = X n-1 e (µ-0.5σ 2 )∆t +σ∆t ϵ n (3.56)
The DOC barrier call option price of a discretely monitored barrier at maturity T is :

C = e -(r -q)T E[g (X n S ) n t ∏ n=1 1 X t ∈[t n-1 ,t n ] ∈A n ] (3.57)
where r, q are respectively an interest and a dividend rates, g (x) = (x -K) + is a payoff function and the set A n in the case of a DOC barrier call option :

A n = inf t n-1 ≤t ≤t n {t : X t > B}
We use continuity correction that was proposed in [1] : B = B exp -0.5826σ∆t . HFMC estimator to compute DOC call option is given by : b

C HFMC = e -(r -q)T 1 n S n S ∑ m=1 ( g (X (m) n t ) n t ∏ n=1 e (-H (X (m) n+1 ,P (m) n+1 )+H (X (m) n ,P (m) n )∆t 1 X (m) n ,X (m) n+1 ∈A n ) (3.58) CHAPITRE 3.

HAMILTONIAN FLOW SIMULATION OF RARE EVENTS

Monte Carlo estimate is given by : b

C MC = e -(r -q)T 1 n S n S ∑ m=1 ( g (X (m) n t ) n t ∏ n=1 1 X (m) n ∈A n ) (3.59)
The IPS estimator is given by : b

C IPS = e -(r -q)T 1 n S n S ∑ m=1 ( g ( b X (m) n t ) n t ∏ t =1 W n-1 (X m n-1 )1 b X n ∈A n )
In the context of a rare event, we chose the barrier level at 65, with an initial price X 0 = 100, Strike K = 100, interest rate r = 0.1, volatility σ = 0.3, T = 0.5 and zero dividends q = 0. In the table 1 and 2, Hamiltonian flow MC, MC and IPS are presented. We used 50000 and 75000 particles with 750 equally spaced time steps in Table 1 and Table 2.

It is very important to choose the number and the size of leapfrog steps. We chose them such that the acceptance probability α is bigger than 0.8.

We compare each approach by estimating the standard deviations, root mean squared error (RMSE), bias, relative mean squared error(RRMSE), time required to compute each estimate and the figure of merit (FOM). We run 20 MC experiments. The RMSE estimator is given by :

RMSE = v u u t 1 M s M s ∑ l =1 ||C -b C l || 2 (3.60)
where C is price computed analytically, b C l are Monte Carlo estimates and M s is the number of Monte Carlo experiments.

The RRMSE is computed using the following formula :

RRMSE = RMSE b C (3.61)
To measure the efficiency of each method presented in the article, we will use the figure of merit(FOM) :

FOM = 1 R 2 × CPU t (3.62)
where CPU t is CPU time need to compute the estimator and R is a relative error, which is the measure of a statistical precision :

R = St .d ev C ∝ 1 p n S (3.63) where C = ∑ M s l =1 b C l .
We run 20 independent Monte Carlo experiments for each estimate. Since IPS and the simulation of Hamiltonian dynamics requires more time to compute an estimate, we use the figure of merit to compare three approaches. From the table 1 and 2 we can observe that HFMC demonstrates standard deviations, bias and relative RMSE.

[H] TABLEAU 3.1 -DOC Barrier option estimates statistics. B = 65,X 0 = 100, K = 100, r = 0.1,σ = 0.

Conclusion and Further Research

We proposed an importance sampling algorithm based on the simulation of Hamiltonian system, that generates Markov Chain that follows along the gradient of the target distributions over large distances of the state space, while producing low-variance samples.

From the simulated results we saw that HFMC allows efficiently estimate rare event probabilities, which we tested on the case of DOC Barrier options. Its estimates show lower variance and bias than that of MC and IPS.

It will interesting to adapt a stochastic gradient Hamiltonian Monte Carlo algorithm [6], when one can avoid computing the gradient at each simulations. Taking into account the big data problem and the necessity of online estimations, we can get sufficient improvements. Another extension is the adaptation to the Riemann Manifold Hamiltonian Monte Carlo [14], when we can create a statistical manifold and tune HFMC by computing explicitely the mass matrix M in the kinetic energy of the algorithm.

In the next article we will show the performance of mixed IPS and Hamiltonian flow Monte Carlo. It will allow faster explore the state space on the one hand, and push trajectories into rare event area on the other hand. By resampling we can reduce the correlation between generated from Hamiltonian system Markov chains.

Introduction

Optimal transport problem that was formulated by Monge in the XVIII century, then reformulated fifty years ago by Kantorovich and it has recently been rediscovered by C. Villani [12]. This problem then was applied in different contexts, for example in option pricing [11].

Particle methods which were extensively researched by P. Del Moral in [6], [5] and [7] allow to find so-called "optimal transport". For this purpose a set of discrete weighted samples, i.e. particles, is used to approximate an importance measure, and then to predict a posterior distribution by propagating the set of particles until we get an estimate.

Another approach has been proposed by Daum's et al. [2], [3] that allows the reduction of the number of particles we need in order to get a tolerable level of errors in the filtering problem. The main idea behind this method is the evolution in homotopy parameter λ (a "pseudotime") from prior to the target density. They introduced a particle flow, in which particles are gradually transported without the necessity to randomly sample from any distribution. This approach as an optimal transport problem allows optimally move the set of particles according to Bayes' rule. In other words, the particles are progressively transported according to their flow. One can in this way reduce the number of needed samples, since the variance and bias of the estimator is lower and as a result reduce the computational burden in both the estimation and the prediction steps.

In this paper we adapt homotopy transport in Stein-Stein stochastic volatility model [9] to price a European option and extend Daum's et al. method by reweighing the generated particles' trajectories that allows to efficiently transport the particles from a prior transition density to a posterior one under the measurement impact. The idea of transportation and reweighing mechanism is to transport particles through the sequence of densities that move the least during the synthetic time until they reach the posterior distribution. By regenerating particles according to their weight at each time step we are able to direct the flow and further minimize the variance of the estimates. The transportation of particles can be understood as a geodesic flow in a convex subset of a Euclidean space.

We show that homotopy transport allows to significantly reduce the variance compared to a particle filtering technique. Path reweighing allows further reduce both the variance and the bias of estimators.

The rest of the article is organized as follows. Section 2 formulates the problem of computing the expectation when we have partially observed variables and shows the solution using particle filter method. Section 3 formulates the problem defined in section 2 in the context of optimal transport and presents the homotopy transport approach to solve the problem. Section 4 shows the mixture of homotopy tranport and path reweighing approach and, actually, extends the method proposed in section 3. Section 5 provides numerical results. Section 6 concludes.

Particle Filtering

Problem formulation

Many problems arises in financial applications when one has to compute expectations with partially observed information. A simple example is an option pricing with hidden volatility dynamics. Assume that we denote by {Y t } t ≥0 ∈ R n Y asset returns, which are observed from the dynamics of prices, while the hidden factor {X t } t ≥0 ∈ R n X is unobservable. Let (Ω, F , P) be a probability space and the set of observed data at time t be (F t ) REWEIGHING IN FINANCE is a filtration generated by a process (Y t ) t ≥0 .

The classical problem, where particle filtering is applied, is to extract a sequence of hidden variables X t . It is formalized in the following way, given an initial R n X -dimensional random variables x 0 with distribution P x 0 , then for t ∈ N :

{ X t = f (X t -1 , ϵ t ) Y t = h(X t , Y t -1 , η t ) (4.1)
where the first equation is the hidden process, with ϵ t : Ω → R n X are i.i.d random variables, the map f : R n X → R n X is B(R n X ) -measurable. The second equation is called a measurement model with η t : Ω → R n Y are i.i.d. random variables and the map h : R n

X × R n Y → R n Y is B(R n X ) ⊗ B(R n Y ) -measurable.
Given above stochastic dynamical system, we would like to compute the following conditional expectation :

E[z(X t )|F t ] = 1 Z ∫ ν(d x 0:t )ρ t (x 0:t , Y 1:t )z(x t ) (4.2)
with a distribution of X 0:t :

ν(d x 0:t ) = p 0 (d x 0 ) t ∏ l =0 k l (X l -1 , X l )µ(d x l ) (4.3)
and normalizing constant Z :

Z = ∫ ν(d x 0:t )ρ t (x 0:t , Y 1:t ) (4.4)
where (X t ) t ≥0 forms a Markov Chain in (R n X , B(R n X )) with transition density k t : R n

X × R n X → R n X
+ with respect to the measure µ(d x). The random variables (Y t ) t ≥0 in (R n Y , B(R n Y )) are conditionally independent given (X t ) t ≥0 with transition density (likelihood) ρ t : R n

X × R n Y → R n Y
+ with reference measure γ. Intuitively, one can think that we can use naive Monte Carlo technique to approximate (4.2) :

E[z(X t )|F t ] ≈ ∫ (M n X P)(d (x 0:t )z(x t ) = 1 n X n X ∑ i =1 f (X (i ) t ) (4.5)
where the sampling operator

M n X ν = 1 n X ∑ n X i =1 δ X (i )
, X ∈ R n X and ∀i = 1,...,n X , X (i ) are i.i.d.

draws from ν.

The problem with naive Monte Carlo Sampling lies in the fact that we don't know how to sample from conditional distribution P(d (x 0:t )) = P(X 0:t ∈ d x 0:t |Y 1:t ). Moreover, computation of normalization constant Z is a big challenge.

There is a lot of research made to tackle this problem, for example [6], where the problem is transformed from a partially observed to a fully observed, by introducing a so called filtering distribution, that links observed and latent variables and recursively updates it.

Definition Conditional probability (filtering distribution) Ξ t = P(X t ∈ •|Y 1 , ..., Y t ) with prior X 0 ∼ p 0 can be computed sequentially :

Ξ t z = ∫ Ξ t -1 (d x t -1 )k t (X t -1 , x)µ(d x)ρ t (x, Y t )z(x) ∫ Ξ µ t -1 (d x t -1 )k t (X t -1 , x)µ(d x)ρ t (x, Y t ) (4.6) with Ξ 0 = p 0 and Ξ t z = ∫ P(x 0:t ∈ d x 0:t |Y 1 , ..., Y t )z(x 0:t )
Denote the corresponding values of the hidden process as (X 0 , ..., X t ) and the values of the measurement process as (Y 0 , ..., Y t ). If there exists an absolutely continuous probability measure P ≪ Q, than for t = 0,...,N we have :

E P [z(X t )|F t ] = E Q [z(X t ) d P d Q (X 0:t )|F t ] (4.7)
An importance measure Q can be chosen arbitrarily as soon as the continuity of the measure is preserved. But usually in a sequential importance sampling literature it is common to see the approximation of Q, given that there exists an absolutely continuous importance kernel e K t , such that for K ≪ e K t as :

Q(B) = M ∑ i =1 ω (i ) t e K t (X (i ) t -1 , A i ), B ∈ B(R n X ) (4.8)
where

A i = {X t ∈ R n X |1 B (X (i ) t -1 , X t ) = 1}, (ω (i ) t -1 ) M i =1
is the weight function, and for i = 1,...,n X , (X (i ) 0 , ..., X (i ) t ) are independent trajectory realizations. Now assume that the prior and sampling kernels K t and e K t have densities k t and e k t with respect to the measure µ, ∀t = 1,...,T. For 0 < ... < t , the Radon-Nikodym derivative in (4.7) is :

d P d Q (X 0:t ) = 1 Z ρ 1 (X 1 , Y 1 ) k 1 (X 0 , X 1 ) e k 1 (X 0 , X 1 ) ...ρ t (X t , Y t ) k t (X t -1 , X t ) e k t (X t -1 , X t ) (4.9)
where the importance measure is given by : Q(d x 0:t ) = p 0 (d x 0 ) e k 1 (x 0 , x 1 )µ(d x 1 )... e k t (x t -1 , x t )µ(d x t ) (4.10)

Observe, that we still can not compute a normalization constant Z, otherwise to compute the filtering distribution Ξ t will not be a problem, so we will need to apply normalized operator M N x to approximate filtering distribution :

E P [z(X 0:t )|F t ] ≈ ∫ M n X P(d x 0:t )z(x t ) = n X ∑ i =1 e ω (i ) t z(X (i ) t )δ X (i ) t (d x t ) (4.11)
where the normalized importance weight function :

b ω (i ) t (X (i ) t ) = ω (i ) t (X (i ) t ) ∑ M i =1 ω ( j ) t (X ( j ) t ) (4.12)
and an unnormalized weight is given by :

ω (i ) t (X (i ) t ) = t ∏ l =1 ρ l (X (i ) l , Y l ) k l (X (i ) t -1 , X (i ) l ) e k l (X (i ) l -1 , X (i ) l ) (4.13)
Observe that importance weights { b

ω (i ) t } n X i =1 are positive and ∑ n X i =1 b ω (i ) t = 1.
Since Particle filters showed weight degeneracey as number of time steps increased, [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF] proposed a resampling step to the algorithm, which can be described by the following nonlinear equation :

Ξ t = Φ t Ξ µ t -1 with Ξ 0 = p 0 (4.14)
where the nonlinear operator Φ t is given by :

(Φ t ν)z = ∫ ν(d x t -1 )k t (X t -1 , x)µ(d x)ρ t (x, Y t )z(x) ∫ ν(d x t -1 )k t (X t -1 , x)µ(d x)ρ t (x, Y t ) (4.15)
The action of the operator Φ can be schematically described as :

Ξ t -1 Mut at i on -------→ M Ξ t -1 Rewei g hi ng ---------→ Ω t M Ξ t -1 (4.16)
where the mutation operator M is given by

(M ν)(z) = ∫ ν(d x t -1 )p(x t -1 , x)µ(d x)z(x) (4.17)
and the reweighing operator Ω t has the form

Ω t (ν)z = ∫ ν(d x)ρ t (x, Y t ) f (x) ∫ ν(d x)g (x, Y t ) (4.18)
After the reweighing step we get the following approximation of the filtering distribu-

tion Ξ t -1 : b Ξ t -1 = n X ∑ i =1 e ω (i ) t -1 δ X (i ) t -1 (4.19)
where {X

(i ) t -1 } n X i =1 ∼ M b Ξ t -2
. We see from above equations that n X particles are sampled from an empirical distribution b Ξ t , i.e. it is itself defined through n X particles.

Let us give the intuition behind the reweighing step. The idea behind it is in the fact, that at this step particles with low weights have lower probability to be sampled compared with particles with high importance weights. Consequently, in this step particles with low weights will be neglected, while particles with high weights will be sampled more frequently.

Particle Filtering Algorithm

The algorithm allows to approximate Ξ

Ξ t = b Φ t b Ξ t -1 with b Ξ 0 = p 0 (4.20)
where b

Φ t := Ω t M n X M . It consists of three steps : b Ξ t -1 Mut at i on -------→ M b Ξ t -1 Sampl i ng -------→ M n X M b Ξ t -1 Rewei g hi ng ---------→ Ω t M n X M b Ξ t -1 (4.21)
At time t = 0, we generate M i.i.d. random variables from the prior distribution. For t = 1,...,N -1 we propagate X t ∈ R n X according to the dynamics of the hidden process, update the measurement, to get a couple of random vectors (X t +1 , Y t +1 ) in the first step. Resample particles according to their probability weights b ω t +1 (X t +1 ) and set resampled particles b X t . At the final time step t compute the estimate of (4.6) :

b C PF = 1 n X n X ∑ i =1 z( b X (i ) t ) b ω (i ) t -1 (X (i ) t -1 )δ X (i ) t -1 (
(i ) t -1 } n X i =1 ∼ M b Ξ t -2 .
Algorithm 10 : PF Algorithm

1 Initialization : i = 1,...,n X -#(simulations), t = 1,...,T -#(time steps) 2 Draw {X (i ) 0 } n X i =1 from the prior p 0 (x). Set {ω (i ) 0 } n X i =1 = 1 n X ;
3 for t = 1,...,N do 4 for i = 1,...,n X do 5

Propagate particles using state equation X

(i ) t = f (X (i ) t -1 , Y (i ) t -1 , ϵ t ) ; 6 Measurement update : Y t = h(X (i ) t , Y (i ) t -1 , η t ) ; 7
Compute effective sample size M e f f :

1

∑ M i =1 (ω (i ) t ) 2 ; 8 if M e f f < M or k < N then 9 Resample using weight b ω (i ) t (X (i ) t ) ω (i ) t (X (i ) t ) 1 n X ∑ M j =1 ω ( j ) t X ( j ) t ) 10 end 11 end 12
Set resampled particles as b X (i ) t 13 end Despite the advantage of sampling from highly non-linear and non-gaussian filtering distributions, we need to mention its limitations. In fact, today we have to deal with highdimensional data, as it was shown in [1], [14], [15], the collapse of weights occurs unless the sample size grows super-exponentially. Homotopy transport allows us to sample efficiently in high-dimensional framework, while avoiding the explosion of the sample size.

Homotopy Transport

The classical optimal transport problem is to find over all maps T : R n X → R n X , such that for X ∼ P, T (X) ∼ Q and T ∈ C 1 ; which optimizes the following criterion :

inf T E[||T (X) -X|| 2 ] s.t. Q = T ♯ P (4.23)
In other words, we would like to find a continuous transformation that minimizes the distance between measure P and measure Q among all these that pusheforward a prior measure P towards a measure Q. In the context of filtering problem we would like to find a transformation T , that transport particles from a sampling measure P to Q :

E Q [z(X t ) d P d Q (X 0:t )|F t ] = E P [z(T (X t ))|F t ] (4.24)
One can solve this problem using variational methods [8].

For the sake of exposition we represent posterior distribution, presented in the form of a normalized importance weight in the following way :

ψ(X t |Y t ) = 1 Z t p(X t |Y t -1 )ρ(Y t |X t ) (4.25)
where Y t = (Y 0 , ..., Y t ), the prior is p(X t |Y t -1 ), the likelihood is ρ(Y t |X t ) and Z t is a normalization factor :

Z t = ∫ p(X t |Y t -1 )ρ(Y t |X t )d X t .
Actually, the equation (4.25) is equivalent to the normalized importance weight in the eq. (4.12). Now, if we consider a continuous map T : R n X → R n X , then :

ψ(T (X t )|Y t ) = 1 Z t p(T (X t )|Y t -1 )ρ(Y t |T (X t )) (4.26)
Homotopy gradually modifies the prior density into the posterior density, as a scaling parameter λ ∈ [0, 1] increases from 0 to 1. In other words, by iterating we will transport homotopy ψ(X t ,λ |Y t ) to a true posterior ψ(X t |Y t ), while minimizing the cost of transport. There are several conditions that homotopy has to satisfy. First, at λ 0 we should have our prior, i.e. ψ(x t ,λ 0 |Y t ) = p(X t ) and at some point λ 0→1 , we will get approximation of our posterior density. Define a new set of density functions :

ψ(X t ,λ |Y t ) := ψ(X t |Y t ), p(X t ,λ |Y t -1 ) := p(X t |Y t -1 ), ρ(Y t |X t ,λ ) λ := ρ(Y t |X t ,λ ) and Z λ := ∫ p(X t ,λ |Y t -1 )ρ(Y t |X t ,λ ) λ d x λ
, so that homotopy is defined as :

ψ(X t ,λ |Y t ) = 1 Z λ p(X t ,λ |Y t -1 ) | {z } pr i or ρ(Y t |X t ,λ ) λ | {z } l i kel i hood (4.27)
In order to simplify the calculation we take the logarithm of homotopy :

Ψ(X t ,λ |Y t ) = G(X t ,λ ) + λL(X t ,λ ) -log Z λ (4.28) where Ψ(X t ,λ ) = log ψ(X t ,λ |Y t ), G(X t ,λ ) = log p(X t ,λ |Y t -1 ), L(X t ,λ ) = log ρ(Y t |X t ,λ
). The dynamics of homotopy transport in the artificial time λ is known as l og -homotopy [2]. In some sense, the dynamics of transport will be given by the flow movement in the aritficial time λ, so we will look for a flow d x dλ that rules the movement of particles following log-homotopy.

If we assume that in pseudo-time λ, the flow d x d λ follows the following SDE :

d X t ,λ = g (X t ,λ )d λ + η(X t ,λ )d W λ (4.29)
where W λ is a vector field that pushes forward particles from prior to posterior distribution.

We impose the following assumptions :

I/ The
∂ψ(X t ,λ ) ∂λ = -t r [ ∂ ∂X t ,λ (g (X t ,λ )ψ(X t ,λ )) ] + 1 2 t r [ ∂ ∂X t ,λ Q(X t ,λ ) ∂ψ(X t ,λ ) ∂X t ,λ ] (4.30)
where Q(X t ,λ ) = η(X t ,λ )η T (X t ,λ ) is the diffusion tensor of the process, and t r (•) is a trace operator. The forward Kolmogorov equation is used to relate the flow of particles

d X t ,λ d λ
with the evolution of log-homotopy as λ 0→1 , under the diffusion process.

∂ψ(X t ,λ ) ∂λ

= -t r [ ψ(X t ,λ ) ∂g (X t ,λ ) ∂X t ,λ + g (X t ,λ ) T ∂ψ(X t ,λ ) ∂X t ,λ ] + 1 2 d i v [ ∂ ∂X t ,λ Q(X t ,λ ) ∂ψ(X t ,λ ) λ ] = = -ψ(X t ,λ )t r [ ∂g (X t ,λ ) ∂X t ,λ ] -g (X t ,λ ) T ∂ψ(X t ,λ ) ∂X t ,λ + 1 2 d i v [ ∂ ∂X t ,λ Q(X t ,λ ) ∂ψ(X t ,λ ) ∂X t ,λ ] (4.31)
where d i v(•) is a divergence operator. On the other hand if we take the derivative of equation (4.28) with respect to λ, we have :

∂Ψ(X t ,λ ) ∂λ = L(X t ,λ ) - ∂ ∂λ log Z λ (4.32)
Since Ψ(X t ,λ ) is a composition of two functions, we will need to use the chain rule :

∂Ψ(X t ,λ ) ∂λ = 1 ψ(X t ,λ )
∂ψ(X t ,λ ) ∂λ

∂ψ(X t ,λ ) ∂λ = ψ(X t ,λ ) [ L(X t ,λ ) - ∂ ∂λ log Z λ ] (4.34)
Observe that (4.31) and (4.34) are identical, so by equating and dividing on ψ(X t ,λ ) we get :

L(X t ,λ ) - ∂ ∂λ log Z λ = -g (X t ,λ ) T 1 ψ(X t ,λ ) ∂ψ(X t ,λ ) ∂X t ,λ - -t r [ ∂g (X t ,λ ) ∂X t ,λ ] + 1 2ψ(X t ,λ ) d i v [ ∂ ∂X t ,λ Q(X t ,λ ) ∂ψ(X t ,λ ) ∂X t ,λ ] (4.35)
In [4], authors propose to take the derivative of (4.35) with respect to X t ,λ in order to find explicitely the equation of flow on the one hand, and to get rid of the normalization constant Z λ that lead to instabilities on the other hand.

∂L(X

t ,λ ) ∂X t ,λ = -g (X t ,λ ) T ∂ 2 Ψ(X t ,λ ) ∂X 2 t ,λ - ∂Ψ(X t ,λ ) ∂X t ,λ ∂g (X t ,λ ) ∂X t ,λ - ∂ ∂X t ,λ t r [ ∂g (X t ,λ ) ∂X t ,λ ] + + ∂ ∂X t ,λ ( 1 2ψ(X t ,λ ) d i v [ Q(X t ,λ ) ∂ψ(X t ,λ ) ∂X t ,λ ]) (4.36)
Observe that we get a highly nonlinear PDE. We use the solution found in [2] and [3], which states that if we can find a vector field g (X t ,λ ) and diffusion tensor Q(X t ,λ ), such that sum of the last three terms in (4.36) are equal to zero. The PDE, then simplifies to :

∂L(X t ,λ ) ∂X t ,λ = -g (X t ,λ ) T ∂ 2 Ψ(X t ,λ ) ∂X 2 t ,λ (4.37)
Using the assumption III, i.e. the Hessian matrix

∂ 2 Ψ(X t ,λ ) ∂X 2 t ,λ
is non-singular, we get explicitely the flow g (X t ,λ ) :

g (X t ,λ ) = - [ ∂ 2 Ψ(X t ,λ ) ∂X 2 t ,λ ] -1 [ ∂L(X t ,λ ) ∂X t ,λ
] T (4.38)

Homotopy Transport Algorithm

Sampling from the prior. First we generate M i.i.d random variables X (i ) t from the prior density p 0 (x), initialize pseudo-time λ and set the state variables that will be transported as : X (i ) t ,λ = X (i ) t |t -1 .

Transportation Stage. For t = 2,...,N -1, compute the derivative with respect to X t ,λ of the measurement function. If h is non-linear, a second order Taylor expansion at X t ,λ allows speeding up the calculation by linearizing the first derivative. After that, update the pseudo time by setting : λ = λ + ∆λ.

Compute the flow g (X (i ) t ,λ ). Note, that the first Hessian can be derived by twice differentiating a log-homotopy equation (4.28) :

∂ 2 Ψ(X (i ) t ,λ ) ∂X 2 t ,λ = ∂ 2 G(X (i ) t ,λ ) ∂X 2 t ,λ + λ ∂ 2 L(X (i ) t ,λ ) ∂X 2 t ,λ (4.39)
The first term in (4.39)

∂ 2 G(X (i ) t ,λ ) ∂X 2 t ,λ
is estimated by using a sample covariance matrix of t patricles generated form the prior distribution :

∂ 2 G(X t ,λ ) ∂X 2 t ,λ ≈ -b S -1 M x (4.40)
Compute the transportation of particles from the measure P to the measure Q :

X (i ) t ,λ = X (i ) t ,λ + ∆λg (X (i ) t ,λ ) (4.41)
And finally update the state parameter :

Xt = 1 n X n X ∑ i =1 X (i ) t ,λ (4.42)

Maturity.

At the final time interval ]N -1,N] compute the estimator of (4.24) :

b C HT = 1 n X n X ∑ i =1 z( X(i) t ) ( 4
(i ) t = f (X (i ) t -1 , Y (i ) t -1 , ϵ t ) ; 7 Measurement update : Y t = h(X (i ) t , Y (i ) t -1 , η t ) ; 8 Initialize pseudo-time λ = 0 ; 9 Set X (i ) t ,λ = X (i ) t

Homotopy Transport with Particle Reweighing

Taking into account the difficulties one faces in non-Gaussian and high-dimensional problems, the idea of a particle transport without any use of sampling techniques is very elucidating. The next question that arises is whether we can direct the transportation by choosing those particles that have higher probability of reaching rarely visited areas of the state space ? We propose a mixture of homotopy particle transport with a particle reweighing at each time step. The numerical test that we performed on the toy example of a Stein-Stein stochastic volatility model showes that we significantly reduce the variance and bias of our estimator.

The algorithm consists of two steps : first we transport particles according to its flow, and second, we choose those particles that have higher probability of faster exploring the state space.

E e Q [z(X t ) d P d e Q (X 0:t )|F t ] = E P [z(T (X t ))|F t ] = E Q [ z(T (X t )) d P d Q (X 0:t )|F t ] (4.44)
where T is a flow of particles under the pseudotime λ decribed in the section 4.3.1.

By setting X t = (X 0 , ..., X t ), we can express our Radon-Nikodym derivative in a product form :

d P d Q (X 0:t ) = d P d e Q × d e Q d Q (X 0:t ) (4.45)
where the first Radon-Nikodym derivative denotes the transport of particles from a mesure P to a measure e Q, then we choose the particles that have high probability of reaching rare corners of the state space, using d e Q d Q that allows us to reassess the weights of the particles.

As in the section 2, an importance measure Q that will play a resampling to choose the trajectories with higher weight, given that there exists an importance kernel e K t , such that K t ≪ e K t , can be defined as :

Q(B) = n X ∑ i =1 ω (i ) t e K t (X (i ) t , A i ), B ∈ B(R n X ) (4.46)
where the set A i = {T (X t +1 ) ∈ R n X |1 B (X (i ) t , T (X t +1 )) = 1}. Assuming, that the prior and sampling kernels K t and e K t have densities k t and e k t respectively, then the Radon-Nikodym derivative is

d e Q d Q (X 0:t ) = t ∏ l =0
ρ l (T (X l ), Y l ) ω l -1 (X l -1 )k l (X l -1 , T (X l ))

ω l -1 (X l -1 ) e k l (X l -1 , T (X l )) (4.47) such that ω t (X t ) = ω (i ) t (X (i ) t ) if X t = X (i ) t , and ω t (X t ) = 1 otherwise.

The an unnormalized weight is given by :

ω (i ) t (T (X (i ) t )) = t ∏ l =1 ρ l (T (X (i ) l ), Y l ) k l (X l -1 , T (X (i ) l )) e k l (X l -1 , T (X (i ) l )) (4.48)
So, now we have homotopy transport with particle reweighing estimator :

b C TRW = 1 n X n X ∑ i =1
z(T (X (i ) t )) b ω (i ) t -1 (T (X (i ) t -1 )) (4.49)

PF-Enhanced Homotopy Transport Algorithm

The algorithm can be described by the following scheme, ∀i = 1,...,n X :

X (i ) t Sampl i ng -------→ X (i ) t +1
Tr anspor t at i on ------------→ T (X (i ) t +1 ) = X(i)

t +1
Rewei g hi ng ---------→ Φ( X(i

) t +1 ) = b X (i ) t +1 (4.50)
where Φ is an operator that denotes the resampling mechanism of particles. If we assume that there is a continuous kernel e K t , such that K t ≪ e K t with densities k t and e k t respectively, then we can define a weight function ω (i ) t :

ω (i ) t ( X(i) t ) = t ∏ l =1 ρ l ( X(i) l , Y l ) k l ( b X (i ) l -1 , X(i) l ) e k l ( b X (i ) l -1 , X(i) l ) (4.51)

Detailed Algorithm

Sampling from the prior. As in the section 4.3.1, we start with M particles sampled from the prior distribution p 0 , initialize pseudo-time λ and set the state variables that will be transported as : X (i ) t ,λ = X (i ) t |n-1 . Transportation Stage. Follow steps 6-8 of the Algorithm 2 in the section 4.

ω (i ) t ( X(i) t ) = ω (i ) t ( X(i) t ) ∑ n X i =1 ω ( j ) t ( X (
C TRW = 1 n X n X ∑ i =1 z( b X (i ) t ) b ω (i ) t -

Numerical Applications and Results

As a toy example, we decided to test the algorithms on a Stein-Stein stochastic volatility model. We set log-returns as Y t = log(S t ), then the model takes the following form :

{ d Y t = (µ - X 2 t 2 )d t + X t d B t d X t = κ(θ -X t )d t + σd W t (4.54)
where X t is a volatility process, Y t the dynamics of log-returns, µ is a drift, θ is a longterm variance, κ -the rate of reversion, σ is the volatility of volatility, and B t and W t are two independent Brownian motions, in the sense that 〈d B t , d W t 〉 = 0.

Using the above presented stochastic volatility model, we would like to compute estimates for a European option. For a given interest rate r , maturity T, strike price K, and for a function z(•, x) = max(x -K, 0), the call price of the option is given by : C(X t , Y t ) = B t ,T E P [z(X T , Y T )|F t ] (4.55)

where F t = σ{(Y 0 , ..., Y t )}. We chose Euler-Muruyama discretization scheme, which gives :

{ Y t -Y t -1 = (µ - X 2 t -1 2 )∆t + X t -1 p ∆tϵ t X t -X t -1 = κ(θ -X t -1 )∆t + σ p ∆tη t (4.56)
where ∆t is a discretization size, ϵ t and η t are independent Gaussian variates, N (0, 1).

We compare each approach by estimating the standard deviations, the root mean squared error (RMSE), the bias, the relative mean squared error(RRMSE), the time required to compute each estimate and the figure of merit (FOM). We run 20 Monte Carlo experiments. For l = 1,...,M s the RMSE estimator is given by : RMSE = To measure the efficiency of each method presented in the article, we will use the figure of merit(FOM) [10] :

FOM = 1 R 2 × CPU t (4.60)
where CPU t is the CPU time need to compute the estimator and R is a relative error, which is the measure of statistical precision : b C l We used 20 000 and 40 000 simulations over 64 time intervals for our MC experiments. Table 1. shows that homotopy and reweighted(RW)-homotopy algorithms shows less statistical errors then traditional PF. If we compare homotopy and RW-homotopy, we can see that FOM says that the first is more efficient the latest, due to the fact that we need more time to reweight the paths. Meanwhile RW-homotopy shows less erros and st. deviations.

[H]

Conclusions and Further Research

The estimation of latent variables has a lot of applications in engineering and finance. We provide homotopy based algorithm and its extension with reweighted trajectories that permits to solve the optimal transportation problem.

Numerical results that we applied in European option pricing with stochastic volatility demonstrated the efficiency of the proposed algorithms with respect to error, bias and other stastics. Both algorithms ourperformed particle filtering. The path-reweighing allowed to reduce standard deviations, and in some cases the bias and the RMSE compared to the homotopy transport algorithm.

From our experiments we can observe the following : -Homotopy transport is fast algorithm, which is spectacularily demonstrated in the figure of merit statistics. -Efficiency of homotopy transport algorithm increases as the number of particles increases. -Implementation of homotopy transport requires less effort than a vanilla Monte Carlo algorithm. -Homotopy transport proved to be unbiased estimator.

-Homotopy with path reweighing proved to reduce the bias when the number of particles is small compared to homotopy transport without reweighing. While reweighted homotopy transport approach showed the reduced RMSE and Bias in low-dimensions, the mixture of homotopy transport and bootstrap resampling, the importance weight can converge to unity in high-dimensional problems( [15]) and potentially worsen. In the next article we plan to use the combination of MCMC and homotopy transport methods on real data. It will be also interesting

∂ 2 ψ ∂x 2 = u ′ v -v ′ u v 2 (4.79)

P(

  |b z(a)z(a)| z(a)

  and set X (m) n+1 = x (m) example, consider sampling from a bivariate Gaussian distribution p(x) = N (µ, Σ) with mean and covariance matrix

FIGURE 1 . 1 -

 11 FIGURE 1.1 -HMC sampling from Bivariate Gaussian distribution

CHAPITRE 2 .

 2 PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING

CHAPITRE 2 .

 2 PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING Denote by P(t , T) the discount factor at time t with maturity T. The Libor rate is defined as :

FIGURE 2 . 1 -

 21 FIGURE 2.1 -Left : Positive exposure profiles. Right : EPE(Blue), ENE(Red)

FIGURE 2 . 2 -

 22 FIGURE 2.2 -The paths of CIR process (left) and survival function (right)

CHAPITRE 2 .

 2 PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING

1while s < τ do 4 j = j + 1 ; 5 τ 6 s 7 8 Simulate random signs for symmetricity : 9 r s = r and ; 10 if 16 Simulate 1 λ log(r and ) ; 17 Ξ(j ) = r and ; 18 end 19 for j = 1 : T n do 20 while t j ( j ) < T j do 21

 456789101611712021 Initialization : M -#(simulations), N -#(time steps), j = 0, T n = {T 1 , ..., T n } 2 s = -log(r and ) ; 3 ( j ) = -log(r and ) ; = s + τ( j ) ;Poisson arrival times ι( j ) = s ; Simulate jump times t j ( j ) = r and ; r s < 0.5 then 11 d r ( j ) = 1 ; interarrival times : Υ(j ) = -

CHAPITRE 2 .

 2 PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING jump-adapted time steps. Now consider the following L évy-driven SDE :

CHAPITRE 2 .

 2 PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING

FIGURE 2 . 5 - 8 CHAPITRE 2 .

 2582 FIGURE 2.5 -PMC simulation of OU Jump-diffusion process with exponentially tempered α-stable jumps (left), Survival probability (right) with α = 0.8

Sommaire 3 . 1 3 . 3 3 . 4 3 . 5 3 . 6

 3133343536 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2 Monte Carlo and Interacting Particle System . . . . . . . . . . . . . . . . . 57 3.2.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2.2 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Hamiltonian flow Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.3.1 Markov Chains on a Phase Space . . . . . . . . . . . . . . . . . . . . . 59 3.3.2 Hamiltonian Flow's Integrator and Properties . . . . . . . . . . . . . 61 3.3.3 Hamiltonian Flow Monte Carlo on Rare Events Sets . . . . . . . . . . 62 3.3.4 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Convergence Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.4.1 IPS convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.4.2 HFMC Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Applications and Numerical Results . . . . . . . . . . . . . . . . . . . . . . 67 Conclusion and Further Research . . . . . . . . . . . . . . . . . . . . . . . 69

Algorithm 9 :

 9 Hamiltonian flow Monte Carlo in Rare event setting 1 Initialization : n S -#(simulations), n t -#(time steps) 2 for n = 1,...,n t do 3 for m = 1,...n S do 4

Assumption 3 . 4 . 1 eTheorem 3 . 4 . 2

 341342 δ n := sup x,y ω n (x) ω n (y)< +∞ (3.46) When the assumption (3.4.1) is met for some e δ n , we have the nonasymptotic estimates :

µ t - 1

 1 by the empirical distribution b Ξ µ t -1 compute by the following reccurence equations : b

(4. 33 )

 33 By substituting eq. (4.33) into (4.32) and rearranging the terms :

3 . 1 .

 31 Path Reweighing Stage. Compute the normalized importance weight :

  b

Maturity

  At the time interval ]N -1,N] compute the final Homotopy transport reweighted estimator : b

2 -

 2 b C l || 2 (4.57) where C is the price computed analytically, b C l are Monte Carlo estimates and M s = 20 is the number of Monte Carlo experiments. As a reference price, we used the article by EM Stein [9]. REWEIGHING IN FINANCE Bi as = √ RMSE St .d ev 2 (4.58)where St .d ev are standard deviations of MC estimates. The RRMSE is computed using the following formula :

TABLEAU 4 . 1 -

 41 Stein-Stein Stochastic volatility option price estimates statistics. S 0 = 100, K = 90, r = 0.0953,σ = 0.2, κ = 4, θ = 0.25, V 0 = 0.25, T = 1/2, and dividends d = 0 True price : 16.05, t = 20000,

FIGURE 4 . 1 -

 41 FIGURE 4.1 -Volatily dynamics, PF (Blue), Homotopy (Red), RW-homotopy(Yellow)

CHAPITRE 4 .

 4 OPTIMAL TRANSPORT FILTERING WITH PARTICLE REWEIGHING IN FINANCE [H]

FIGURE 4 . 2 -

 42 FIGURE 4.2 -Zoomed volatilty dynamics. Homotopy (left), RW-homotopy (right)

Table des matières iii Liste des figures v Liste des tableaux vii Chapitre 1 Introduction

« " πάντες αντροποι του είδέναι ορέιγονται ϕύσει. τὰ μετά τὰ ϕισικά ' » Αριστοτέλη « "All men by nature desire to know" » Aristotle

17 end 18 Update state estimate : 19

	t ,λ + ∆λ	dX (i ) t ,λ dλ ;

Chapitre 2 Particle Levy tempered α-stable credit risk model simulation with application to CVA with WWR computation in multi-factor BGM setting

	« "Stabilité première condition du
	bonheur publique. Comment
	saccommode-t-elle avec la
	perfectibilité indéfinie ?" »
	A.S. Pouchkine
	« "Stability -the first condition of
	public happiness. How does it fit
	with indefinite perfectibility ?" »
	A.S.

Puchkin Sommaire 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2 Problem Formulation and Definitions . . . . . . . . . . . . . . . . . . . . 24 2.3 Market Risk and Exposure Profiles . . . . . . . . . . . . . . . . . . . . . . . 25

	2.3.1 BGM Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
	2.3.2 Expected Exposure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
	2.3.3 Multi-curve Interest Rate Modelling . . . . . . . . . . . . . . . . . . . 27
	2

.4 Credit Risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2

.4.1 L évy driven OU Process . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4.2 Simulation of α-Stable Processes and Numerical Schemes . . . . . . 32 2.5 WWR Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.6 Particle Interpretations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.6.1 Interacting Particle System for Rare Events Estimation . . . . . . . . 37 2.6.2 Mean Field IPS Approximations . . . . . . . . . . . . . . . . . . . . . 38 2.6.3 IPS Simulation of the Stochastic Intensity . . . . . . . . . . . . . . . . 41 CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING

  PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING 1,...,ϱ(T)} with intensity l l and jump sizes or marks are represented by the following set : {ξ

d t )

(2.26)

where : ϱ(d x, d t ) is a Poisson random measure that generates a sequence of pairs of random variables {τ i , ξ i , i = 1,2,...,ρ(T)}. The set of jump times is given by {τ

i : Ω → R + , i = CHAPITRE 2.

  33) CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING Assume that corresponding L évy X process could be decomposed as :

  [START_REF] Pykhtin | Pricing Counterparty Risk at the Trade Level and CVA Allocations[END_REF] If we do not consider toy examples, condition (2.38) is difficult to verify. The next proposition from[2] gives a sufficient condition for eq.(2.38).

	Proposition 2.4.3 Condition (2.38) is implied by

  53) CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING Then the approximated version of the L évy-driven SDE have the following form

  .70) In the figures 3 and 4 we see the sample paths of α-stable process with values of α = 0.8 and α = 1.5 and survival probabilities. We simulated so-called almost exact solution given in eq. (2.54) and compared with three discretization schemes : Euler with constant discretization steps ∆t, Euler and Milstein jump-adapted versions. You can observe that for larger values of α, the process resembles Brownian motion and compound Poisson process for low values, that gives the intuition behind the decomposition in eq. (2.56).CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING

FIGURE 2.3 -OU Jump-diffusion process with exponentially tempered α-stable jumps (left), Survival probability (right) with α = 0.8 FIGURE 2.4 -OU Jump-diffusion process with exponentially tempered α-stable jumps (left), Survival probability (right) with α = 1.5

  84) CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING where the importance measure Q n (•, d y)

  TABLEAU 2.1 -||•|| L 2 error estimates of numerical schemes using MC and PMC density approximations, α = 0.8

	CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL
	SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN
	MULTI-FACTOR BGM SETTING
	∆t 2 -6 2.54e-01 0.22 2.62e-01 0.15 2.62e-01 0.14 7.10e-02 0.03 2.16e-01 0.02 MC eu t MC eu j ad p t MC ml s j ad p t PMC eu t PMC eu j ad p t 2 -7 1.55e-01 0.27 1.67e-01 0.16 1.67e-01 0.17 6.05e-02 0.27 8.23e-02 0.16 2 -8 1.85e-01 0.32 3.43e-02 0.22 1.85e-01 0.25 7.88e-02 0.06 1.08e-01 0.05 2 -9 2.74e-01 0.28 2.91e-01 0.19 2.92e-01 0.20 6.63e-02 0.12 1.30e-01 0.10 2 -10 3.19e-01 0.29 3.30e-01 0.20 3.29e-01 0.20 7.14e-02 0.23 7.90e-02 0.21	PMC ml s j ad p 2.15e-01 8.10e-02 1.09e-01 1.30e-01 7.89e-02	t 0.02 0.17 0.05 0.10 0.21
	TABLEAU 2.2 -||•|| L 1 error estimates of numerical schemes using MC and PMC density approxima-tions . α = 0.8
	∆t 2 -6 2.04e-01 0.22 2.09e-01 0.15 2.09e-01 0.14 5.48e-02 0.03 1.656e-01 0.02 1.652e-01 0.02 MC eu t MC eu j ad p t MC ml s j ad p t PMC eu t PMC eu j ad p t t PMC ml s j ad p 6.23e-02 0.03 2 -7 1.03e-01 0.27 1.12e-01 0.16 1.12e-01 0.17 5.25e-02 0.04 6.28e-02 0.03 1.3e-01 0.32 1.36e-01 0.22 1.36e-01 0.25 6.88e-02 0.06 8.32e-02 0.05 8.34e-02 0.05 2 -8 2.1e-01 0.20 5.91e-02 0.12 7.95e-02 0.1 7.95e-02 0.1 2 -9 2.03e-01 0.28 2.1e-01 0.19 2.1e-01 0.20 5.86e-02 0.23 6.32e-02 0.21 6.36e-02 0.21 2 -10 1.95e-01 0.29 2.1e-01 0.20 TABLEAU 2.3 -||•|| L 2 error estimates of numerical schemes using MC and PMC density approxima-tions, α = 1.5 ∆t MC eu t MC eu j ad p t MC ml s j ad p t PMC eu t PMC eu j ad p t PMC ml s j ad p 2.51e-01 0.28 8.1e-01 0.44 8.01e-01 0.49 1.82e-01 0.4 4.34e-01 1.42 4.33e-01 1.36 t 2 -6 3.840e-01 1.27 2 -7 3.548e-01 0.13 6.962e-01 0.33 6.922e-01 0.33 1.694e-01 0.54 3.857e-01 3.21 7.6 2.673e-01 1.9 2 -8 3.203e-01 0.1 4.572e-01 0.25 4.563e-01 0.29 1.966e-01 1.59 2.831e-01 2 -9 3.243e-01 1.05 6.082e-01 2.55 6.011e-01 0.75 2.037e-01 4.01 3.206e-01 10.24 3.232e-01 2.29 2 -10 3.33e-01 0.29 3.77e-01 0.44 4.11e-01 0.42 2.08e-01 15.9 3.16e-01 21.9 3.04e-01 4.47
	∆t 2 -6 2 -7 3.626e-01 0.13 6.964e-01 0.33 6.922e-01 0.33 1.823e-01 0.54 3.861e-01 3.21 MC eu t MC eu j ad p t MC ml s j ad p t PMC eu t PMC eu j ad p t 3.15e-01 0.28 9.89e-01 0.44 9.85e-01 0.49 2.3e-01 0.4 5.23e-01 1.42 7.6 2 -8 3.230e-01 0.1 4.575e-01 0.25 4.563e-01 0.29 2.036e-01 1.5 2.837e-01 2 -9 3.264e-01 1.05 6.083e-01 2.55 6.011e-01 0.75 2.066e-01 4.01 3.212e-01 10.24 3.232e-01 2.29 t PMC ml s j ad p 5.24e-01 1.36 3.840e-01 1.27 2.673e-01 1.9 2 -10 3.34e-01 0.29 3.78e-01 0.44 4.11e-01 0.42 2.10e-01 15.9 3.17e-01 21.9 3.04e-01 4.47

CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING

FIGURE 2.6 -PMC Euler, Euler jump-adapted, Milstein jump-adapted strong errors, α = 0.8, #(particles)= 1000. TABLEAU 2.4 -||•|| L 1 error estimates of numerical schemes using MC and PMC density approximations . α = 1.5

FIGURE 2.7 -MC vs PMC Euler (left) and jump-adapted (right) strong errors. α = 0.8, #(particles)= 1000, #(simulated paths)= 5000 CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING

Table 7 .

 7 gives CVA with WWR values, we can see that Basel III approximation is overvalued for initial intensity values lower than 0.03 and for 0.04. CHAPITRE 2. PARTICLE LEVY TEMPERED α-STABLE CREDIT RISK MODEL SIMULATION WITH APPLICATION TO CVA WITH WWR COMPUTATION IN MULTI-FACTOR BGM SETTING TABLEAU 2.5 -CVA WWR/free for LIBOR SWAP. PMC approach. α = 0.8, #(scenar i os) := #(sc) LIBOR SWAP CVA WWR/free for different values of ρ and λ 0 , α = 0.8. PMC approach.

	λ 0	CVA	CVA WWR CVA WWR(BIII)	CVA	CVA WWR CVA WWR(BIII)
	#(sc)	(st.dev) 10 4	(st.dev) 10 4	(st.dev) 10 4	Θ BIII si m 10 4	(st.dev) 10 5	(st.dev) 10 5	(st.dev) 10 5	Θ BIII si m 10 5
	0	1.02e-04	1.29e-04	1.43e-04	1.11	8.32e-05	1.19e-04	1.16e-04	0.84
		(7.71e-05) (9.93e-05)	-		(7.7e-05)	(1e-04)	-
	0.01	1.03e-04	1.36e-04	1.45e-04	1.07	1.25e-04	1.2e-04	1.75e-04	1.46
		(8.03e-05) (1.02e-04)	-		(6.5e-05)	(1.1e-04)	-
	0.02	9.95e-04	1.06e-04	1.39e-04	1.31	1.24e-04	1.63e-04	1.74e-04	1.07
		(7.55e-05) (1.02e-04)	-		(7.78e-05) (9.12e-05)	-
	0.03	8.62e-05	9.82e-05	1.2e-04	1.22	1.07e-04	1.2e-04	1.5e-04	1.25
		(8.16e-05) (1.07e-04)	-		(8.46e-05) (1.1e-04)	-
	0.04	8.85e-05	1.44e-04	1.24e-04	1.37	1.15e-04	1.46e-04	1.61e-04	1.1
		(7.63e-05) (9.63e-05)	-		(7.79e-05) (9.29e-05)	-
		ρ	λ 0	0	0.01	0.02	0.03	0.04
			0.1	8.6e-05	9.1e-05 1.22e-04 7.52e-05 8.54e-05
			0.2	6.9e-05	8e-05	1.3e-04 8.93e-05 1.14e-04
			0.3	6.8e-05	9.4e-05 7.69e-05 1.13e-04 1.04e-04
			0.4	1e-04	6.2e-05 9.07e-05 9.36e-05 1.01e-04
			0.5	8.2e-05	7.0e-05 9.81e-05 1.15e-04 9.98e-05
		CVA WWR(B III) 1.43e-04 1.45e-04 1.39e-04 1.12e-04 1.97e-04
		CVA	1.02e-04 1.03e-04 9.95e-05 8.62e-05 1.4e-04
	TABLEAU 2.7 -ρ	λ 0	0	0.01	0.02	0.03	0.04
			-0.5	1.25e-04 1.11e-04 9.91e-05 1.21e-04 1.27e-04
			-0.6	8.67e-05 8.76e-05 1.07e-04 1.21e-04 1.18e-04
			-0.7	1.02e-04 1.06e-04 1.16e-04 1.16e-04 1.64e-04
			-0.8	8.12e-05 9.58e-04 1.49e-04 1.56e-04 1.29e-04
			-0.9	1.37e-05 1.25e-04 1.29e-04 1.59e-04 1.79e-04
		CVA WWR(B III) 1.43e-04 1.45e-04 1.39e-04 1.12e-04 1.44e-04
		CVA	1.02e-04 1.03e-04 9.95e-05 8.62e-05 8.85e-04

TABLEAU 2.6 -LIBOR SWAP CVA WWR/free for different values of ρ and λ 0 , α = 0.8. PMC approach.

  3, T = 1/2, and di v = 0 ; δ = 0.0001, #(Leap frog step) : 35. True price : 10.9064, n S = 50000, n t = 750 DOC Barrier option estimates statistics. B = 65,X 0 = 100, K = 100, r = 0.1,σ = 0.3, T = 1/2, and di v = 0 ; δ = 0.0009, #(Leap frog step) : 40. True price : 10.9064, n S = 75000, n t = 750

	Stat	MC	PMC	HFMC
	St. dev.	0.088518965 0.08562686 0.065318495
	RMSE	0.007011127 0.008004332	0.0143
	RRMSE	0.001298078 0.000292621 1.87148E-05
	CPU time	3.7251	4.8432	5.90675
	FOM	4097.9	3387.2	4737.6
		[H]
	TABLEAU 3.2 -Stat	MC	PMC	HFMC
	St. dev.	0.062385996 0.044259477 0.038039517
	RMSE	0.037561882 0.051285344 0.037561882
	RRMSE	0.000355199 0.000240548 0.000129293
	CPU time	2.2626	6.0322	7.6832
	FOM	13475.2	10117.7	10711.0

  .[START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF]) REWEIGHING IN FINANCE Initialization : i = 1,...,n X -#(simulations), t = 1,...,N -#(time steps)

	Algorithm 11 : Homotopy Transport Algorithm
	2 Draw {X (i ) 0 } 3 Set {ω (i ) 0 } n X i =1 = 1 n X i =1 from the prior p 0 (x). n X 4 for t = 1,...,N do
	5	for i = 1,...,n X do

1 6

Propagate particles using state equation X

  |n-1 ;

	10 11 12 13	while λ < 1 do Compute SCM b S M ; Calculate an estimate : X t ,λ = 1 n X Compute the matrix b H = ∂h(X (i ) t ,λ ) ∂X t ,λ	∑ i X (i ) t ,λ ;
			dX (i ) t ,λ dλ = -	[ ∂ 2 Ψ(X (i ) t ,λ ) ∂X 2 t ,λ	]-1 [ ∂L(X (i ) t ,λ ) ∂X t ,λ	] T	;
	19	Xt = 1 n X	∑ n X i =1 X (i ) t ,λ
	20	end
	21 end

14

Update the time : λ = λ + ∆λ ; 15

Calculate the flow 16

Transport particles according to its flow : X

(i )

t ,λ = X (i )

t ,λ + ∆λ d X (i ) t ,

λ d λ ; 17 end 18 Update state estimate :

  Homotopy Transport with Particle Reweighing Algorithm 1 Initialization : i = 1,...,n X -#(simulations), t = 1,...,T -#(time steps) Follow steps 6-8 of the Algorithm 2 in the section 4.3.1. Follow steos 7-12 of the Algorithm 1 in the section 4.2.2

		1 ( X(i) t -1 )	(4.53)
	Algorithm 12 : 2 Draw {X (i ) 0 } n X i =1 from the prior p 0 (x). 3 Set {ω (i ) 0 } n X i =1 = 1 n X 4 for t = 1,...,N do
	5	for i = 1,...,n X do
	6
	8	end
	9 end

7

  16.05, t = 40000,

	M = 64
	Stat	MC	PF	Homotopy RW-Homotopy
	St. dev.	0.070351719 0.060799052 0.048943672	0.045246118
	RMSE	0.130446299 0.079273246 0.04921257	0.045762201
	Bias	0.109849318 0.050869665 0.005137504	0.006853309
	RRMSE	0.001067583 0.000392831 0.00015101	0.000130578
	CPU time	0.278895	0.54737	0.26618	0.581495
	FOM	184049.069 126479.8136 403391.758	216062.7397

  to investigate the homotopy transport on non-gaussian examples. ,N (Xt -1 )∇ X m t ,N (X t -1 ) + (Y tm t ,N )∇ X σ m t ,N )(-2X 3 t -1 ∆t 2 + 2(Y tm t ,N )X t -1 ∆t) (4.76) u ′ = 2X 2 t -1 ∆t 2 ( X 2 t -1 ∆t -(Y tm t ,N ) ) + (Y tm t ,N ) ( -6X2t -1 ∆t 2 + 2∆t((Y tm t ,N ))

	BIBLIOGRAPHIE
	∂ψ ∂x	(X) =	1 2	N ∑ t =1 -(2∆t -(Y t -m t ,N )(2σ	p t p t ,N (X t -1 )) σ p t ,N (X t -1 ) 2	)	(4.75)
				u = (Y ) (4.77)
				v = X 4 t -1 ∆t 2 , v ′ = 4X 3 t -1 ∆t 2		(4.78)

t

2.7 Numerical Analysis and Discussion . . . . . . . . . . . . . . . . . . . . . . 44

Algorithm

The algorithm can be described by the following scheme, ∀m = 1,...,M :

(2.100)

At each time step n = 1,...,N we draw M independent random variables from the density q n (•, X (i ) n-1 ) to construct M particles b X (m) n = ( b X (1) n , ..., b X (M) n ). Given generated particles, we select independently M particles b X n = ( b X (1) n , ..., b

n ) with replacement of rejected particles according to their probability weights :

to sample

n .

14 end 15 end

IPS Simulation of the Stochastic Intensity

In this section we formulate the problem of a rare event estimation and show on the example of a stochastic intensity the application of IPS to estimate survival probabilities.

When we deal with credit risks most of the problems are formulated as occurrences of rare events, in the sense that P(τ > t ) is very small, when t is large. If one can transform from a measure P to the measure Q , using Radon-Nikodym derivative transformation, so that in the new measure Q, t is not too large, then one can use central limit theorem approximation.

Interacting particle system, described in previous sections allows us to approximate a rare event probability by sequentially regenerating particles according to their probability weights.

Abstract

In this article we present a sampling method, Hamiltonian flow Monte Carlo, which is based on the simulation of Hamiltonian system to estimate rare events. The classical way to deal with this problem is the change of measure, so that our samples can reach rare events area, however, in practise it is hard to find an optimal importance density and the Monte Carlo simulation performs terribly, when the rare-event probability goes deep into the tail. The application of Hamiltonian dynamics allows to estimate rare events and sample from target distributions defined as the change of measures. Our approach uses a physical simulation of moving particles with momentum under the impact of the energy function to propose Markov chain transitions, that allows rapidly explore a state space. Its fast exploration can be explained by the fact that it extends the state space by an auxilliary momentum variables, P, and then runs a physical simulation to move long distances along probability contours in the extended state space.

The estimates demonstrated a variance reduction of the presented algorithm and its efficiency with respect to a standard Monte Carlo and interacting particle based system(IPS). We tested the algorithm on the case of the barrier option pricing.

Algorithm

The algorithm can be described by the following scheme, ∀m = 1,...,n S :

At each time step n = 1,...,n t we draw n S independent random variables from the den-

n-1 ). Given generated particles, we select, or draw independently n S particles X (m) n = (X (m) 0 , ..., X (m) n ) with replacement of rejected particles according to their probability weights :

And at time step n t , we get the following IPS estimator :

1 Initialization : n S -#(simulations), n t -#(time steps), X 0 -initial value

to sample X (m) n .

14 end 15 end

Hamiltonian flow Monte Carlo

Markov Chains on a Phase Space

From section 2 we know that one of the ways to deal with rare event probabilities is to change a measure :

We will use Hamiltonian system to generate Markov Chain and approximate a Radon-Nikodym derivative dP d Q . Hamiltonian flow Monte Carlo uses a physical simulation of moving particles with momentum under the impact of the energy function to propose Markov chain transitions, that allows rapidly explore state space. Its fast exploration can be explained by the fact that it extends the state space by an auxilliary momentum variables, P, and then runs a physical simulation to move long distances along probability contours in the extended state space.

We remind that, given Markov Chain {X l } l ≥0 , Birkhoff theorem says that

where ϱ is the expectation of f (X) with respect to the unique invariant distribution π of Markov Chain.

Hamiltonian Flow's Integrator and Properties

We will use a configuration space M with periodic boundary conditions. Each point on M will be a set of n S particles : X (1) , ..., X (n S ) and a generic momentum space R n S , in this case the cotangent space is given by T

Ξ u (X 0 , P 0 ) is the solution to the Hamilton's equation :

Hamiltonian system has three main properties : reversibility, conservation of energy and volume preservation.

Sympletic Integration of Hamiltonian Equations

In most cases we can not compute Hamiltonian flow in closed form and that is why we need to discretize the system (3.24). To make sure that we can preserve symplecticness and reversibility, we will discretize using leap-frog integrator, which is a symplectic integrator of Hamiltonian system.

Split Hamiltonian (3.22) into 3 parts :

Taking each of these terms separately to be Hamiltonian function of Hamiltonian system gives rise to equations of motion with trivial dynamics.

where ∆u is the discretization size of Hamiltonian.

Consider a concatenation of three maps :

where Ξ ∆u,H 1 : (X(0), P(0)) → (X(∆u), P(∆u)). Similarily, Ξ ∆u,H 1 = Ξ ∆u,H 3 , and Ξ ∆u,H 2 is calculated to be position update. Since the energy is preserved by the flow, the trajectories evolve on the submanifold of constant energy :

where E 0 = H (X 0 , P 0 ) is the energy of the initilized data.

Hamiltonian Flow Monte Carlo on Rare Events Sets

Let H (X, P) be a Hamiltonian function on R 2n S , where X is a potential, and P is a momentum variable of Hamiltonian system. The algorithm consists of two steps, first sampling from prior distribution values for potential and momentum and then a physical simulation of Hamiltonian dynamics. To make sure that at the end of each physical simulation of time step n + 1 we will have a probability measure, i.e. values will not exceed 1, we will use a Metropolis-Hastings test α n+1 , by choosing the minimum between 1 and the ratio of generated values of potential at time steps n + 1 and n, which is an acceptance probability of potential simulated by Hamiltonian dynamics. If we extend the state space X = {X 1 , ..., X n } and denote the extended space as e X = {X 1 , ..., X n , P 1 , ..., P n }, we can denote the acceptance probability as :

If we assume that the importance measure Q admits the following importance distribution with a kernel K :

Then, the associated Radon-Nikodym derivative will have the following form :

Assume that at each times step n we have n S sample of r.v. {X (m) n } n S m=1 . Now we can define a transition kernel K as follows.

Definition Consider a mapping Ξ n : e

n+1 ) is given by :

This kernel can be interpreted as the probability to move from the point e X If the proposed step is not accepted, then next step is the same as the current step, i.e. e X (m) n+1 = e X (m) n . This procedure allows as to leave the joint distribution of X (m) n and P (m) n invariant. Volume preservation means that the determinant of the Jacobian matrix of a transformation Ξ n is equal to one.

We will need basic property of symplectic integrators, i.e. reversibility. We refer to [13] for the proof of this result.

Proof Rewrite the kernel K as :

where Definition Assume that the assumption 3.3.2 holds and consider Markov Chain (X (m) n ) n≥1 with an initial prior p 1 (X 1 ) and define the following transition kernel p( e

It means that the point e X (m) n moves to a new point e X (m) n+1 only if it is inside a rare event set A n , otherwise we stay at point e X (m) n .

Invariance of e X (m) n says that for any bounded and measurable function f , the distribution of f (M ( e

Under the kernel M of e X n , the final HFMC estimate is given by :

Algorithm

Hamiltonian function is defined by H (X, P) = Ψ(X) + 1 2 P T M -1 P, where Ψ(X) -is a potential energy function, and the second term is a kinetic energy function with a momentum variable P and mass matrix M. Usually one sets a mass matrix M to be an identity matrix I. The proposed samples are obtained by a physical simulation of Hamiltonian dynamics :

We start by simulating M random variables from a prior X 1 = p 0 (•, X (m) 0 ), which is the density of the underlying SDE and generating M random variables from gaussian distribution for momentum {P (m) 0 } M m=1 .

For each step n = 1,...,N we set x (m) H = X (m) n , P (m) H = P (m) n . The proposed new candidates are obtained after L-leapfrog steps of the simulation of Hamiltonian dynamics and they are defined by x * = x (m) H (L) and P * = P (m) H (L). These new set of proposed candidates are then accepted according to the following Metropolis-Hastings test. First generate uniformly distributed random variable u ∼ U ni f (0, 1), then compute α : α = 1 ∧ e (-H (x * ,P * )+H (x (m) H ,P (m) H ))∆t ;

(3.44)

If proposed candidates (x * , P * ) are accepted, i.e. α > u we set X (m) n+1 = x * , and if they are rejected, i.e. α ≤ u, we set X (m) n+1 = x (m) H . At the end, calculate estimator in (3.42). The main steps of the algorithm are summarized in Algorithm 2.

The Metropolis-Hastings test insures a volume preservation. That explains the fact that we don't need to compute a normalizing constant in our algorithm. Volume preservation means that the absolute value of the Jacobian matrix of the leapfrog integrator is

Abstract

In this paper we show the application of an optimal transportation problem to estimate stochastic volatility process by using a flow that optimally transports the set of particles from the prior to a posterior distribution. Importantly, it is well-known that a naive particle filter has curse of dimensionality, when the quality of weights deteriorates as the number of iteration increases. Homotopy transportation approach allows efficiently estimate stochastic volatility in high-dimensional problems. We also show how to direct the flow to a rarely visited areas of the state space by using a particle method (a mutation and a reweighing mechanism). The efficiency of our approach is demonstrated on the simple example of the European option price under the Stein-Stein stochastic volatility model for which a closed form formula is available. Both homotopy and reweighted homotopy methods show a lower variance, root-mean squared errors and a bias compared to other filtering schemes recently developed in the signal-processing literature, including particle filter techniques.

Appendix Flow related computations

In the classical particle filtering approach the desired estimate is approximate by M pratciles : The transition density is given by : p(X t |X t -1 ) = N (X t ; µ X , σ X ) (4.64)

where µ X = X t -1 + κ(µ X -X t -1 )∆t and σ X = σ 2 X ∆t. The unnormalized posterior is given by :