Calcul de la CVA par système particulaire

P. Cohort (Zeliade) and B. Lange (LPMA Diderot and Zeliade)

Ce notebook décrit la méthode proposé par Henry-Labordère , et en propose une implémentation, qui reproduit les résultats obtenus dans l'article.

L'equation aux dérivées partielles

Le modèle

Une partie A vend à B un produit dérivé sur \(d\) actifs sous-jacent \(X \in \mathbb{R}_{+}^{d}\). On note \(u\) le prix de ce produit avec le risque de contrepartie, engendré par le risque de défaut de B, de payoff \(\psi\) et de maturité \(T\). Dans le cas d'accords de netting sur plusieurs produits dérivés, \(u\) représente la valeur agrégée de tous ces produits. Afin de se couvrir contre le risque de contrepartie engendré par B, A peut trader dynamiquement une obligation risquée émise par B, notée \(O^B\) (toutes les expressions avec un indicatif \(^B\) font référence à la partie risquée).

Les processus \(X_t\) et \(O^B_t\) satisfont sous la probabilité risque-neutre : \[\frac{dX_t}{X_t}=rdt+\sigma(t,X_t).dW_t,\] \[\frac{dO^{B}_t}{O^{B}_t}=(r+\lambda^B)dt-dJ^{B}_t,\]\(W_t\) est un mouvement Brownien sous la probabilité risque-neutre, \(d\)-dimensionnel, \(J^{B}_t\) est un processus de Poisson d'intensité \(\lambda^B\) et \(r\) le taux sans-risque. On suppose que le modèle est complet et sous la probabilité risque neutre \(\exp(-rt)u(t,X_t)\) est une martingale telle que :

\(\mathcal{L}\) représente le générateur d'Itô de \(X\) (pour \(X\) unidimensionnel; \(\mathcal{L}=\partial_x+\frac{1}{2}\sigma(t,x)^2\partial^{2}_{xx}\)) et \(\widetilde{u}\) la valeur du produit dérivé après le défaut de B. Soit dans le cas général \(\widetilde{u}\) est donné par :

\(R\) est le taux de recouvrement et \(M\) est la valeur marché du produit utilisée dans la procédure de clôture du contrat (\(Y=Y^+-Y^-\)). \

Clôture du contrat au moment du défaut: le closeout

La définition de M dépend donc du type de procédure mise en place lors de la clôture du contrat au moment du défaut. Dans le cas du replacement closeout on prend en compte le risque de contrepartie (\(M=u\)) et dans le cas du risk-free closeout on prend uniquement en compte le prix sans risque du produit :

Avec un certain changement de variable on se ramène à deux systèmes d'EDPs plus simples encore : \[\partial_tu+\mathcal{L}u+S^B(u^+-u) = 0,\ \ u(T,x) = \psi(x)\ :\ EDP1\] \[\partial_tu+\mathcal{L}u+\frac{S^B}{1-R}\left((1-R)\mathbb{E}_{t,x}[\psi]^++R\mathbb{E}_{t,x}[\psi]-u\right) = 0,\ \ u(T,x) = \psi(x)\ :\ EDP2\] avec \(S^B=\lambda^B(1-R)\) qui s'identifie au spread de CDS et \(\mathbb{E}_{t,x}[\psi]\) représente le prix sans risque \(M\) dans le cas du risk-free closeout.

Résolution par système particulaire

La méthode que propose P. Henry-Labordère est une résolution des EDPs dans lesquelles on approche la semi-linéarité que provoque la fonction partie positive par un polynôme. On cherche alors à résoudre l'EDP :

avec \(F(u) = \sum\nolimits_{k=0}^{M} \left(\frac{a_k}{p_k}\right)p_ku^k\) un polynôme de degré \(M\). Nous verrons plus bas le choix d'écrire le polynôme de cette façon.

L'auteur montre d'abord qu'on peut trouver une solution dans le cas général où \(S^B\) est une fonction déterministe dans \(\mathbb{R_+}\) et où \(M=\infty\) (avec des contraintes supplémentaires telles que \(|\psi|<1\)). C'est H. P. McKean qui a été le premier à donner une représentation probabiliste de ces EDPs. On peut interpréter cette EDP de façon probabiliste de la manière suivante :

Processus de branchement

Soit une unique particule au point de départ \(X_t\) qui suit la diffusion d'Itô de générateur \(\mathcal{L}\) (générateur de \(X\)) sur \(\mathbb{R}_d\), meurt après un temps exponentiel de moyenne \(S^B\) et donne naissance à \(k\) particules avec une probabilité \(p_k\). Les \(k\) descendants suivent alors des diffusions (de même générateur \(\mathcal{L}\)) indépendantes depuis leur lieu de naissance (lieu de la mort de la particule mère), puis meurent après des temps exponentiels \(S^B\) et donnent naissance à d'autres particules de la même façon, et ainsi de suite...

On appelle ce type de processus un processus de branchement \(d\)-dimensionnel avec un taux de branchement \(S^B\) et il peut aussi être identifié à un arbre aléatoire ou arbre de Galton-Watson pondéré. On note \(\left(\nu_t^1, ..., \nu_t^{N_t}\right) \in \mathbb{R}^{d \times N_{t}}\) les positions des particules en \(t\) et \(N_t\) le nombre de ces particules. Enfin, en \(T\), on compte le nombre \(\omega_k \in \mathbb{N}\) de naissances de type \(k\) dans l'arbre (\(k \in \{0,...,M\}\)). Les descendants sont tirés avec une distribution \(p_k\) arbitraire (une distribution uniforme donnerait \(p_k=\frac{1}{M+1}\) l'auteur propose une distribution optimale; voir plus bas). Enfin, l'auteur donne l'expression de la probabilité \(p_k\) optimale en minimisant la limite de la variance de l'algorithme , il trouve l'expression :

Le résultat principal de est le suivant: on peut montrer sous certaines hypothèses que la fonction

est l'unique solution de viscosité de ()

La force de la méthode est que la solution ne nécessite plus qu'une seule simulation de Monte Carlo pour calculer l'espérance, à l'inverse de la résolution directe nécessitant une double simulation "Monte Carlo de Monte Carlo" .

Implémentation du système particulaire

On défini un classe représentant une particule, et une classe simulant une trajectoire du système jusqu'à une maturité donnée. les temps de naissance des particules sont simulés de façon exacte (pas de discrétisation du temps) suivant une intensité constante.

In [1]:
class Particle:
    def __init__(self, value, t):
        self.value = value
        self.t = t        
In [2]:
from math import log
import numpy as np
import random

class ParticleSystem:
    
    def __init__(self, initialValue, diff, intensity, proba):
        self.initialValue = initialValue
        self.state = [Particle(initialValue, 0)]
        self.diff = diff
        self.intensity = intensity
        self.probaSum = np.cumsum(proba)
        self.omega = [0.0]*len(proba)
                
    def newParticle(self, particle):
        tau = self.defaultTime() 
        v = self.diff.next(particle.value, min(tau, self.T-particle.t))
        return Particle(v, min(particle.t+tau, self.T))
            
    def branch(self, particle):
        k = self.particleNumber()
        output = []
        for i in range(k):
            output.append(self.newParticle(particle))                
        return output
    
    def simulate(self, T):
        self.T = T
        self.omega = [0.0]*len(self.omega)
        self.state = [self.newParticle(Particle(self.initialValue, 0))]
        while len(self.state)>0 and min(self.state, key=lambda p:p.t).t<T:            
            nextState = filter(lambda p:p.t==T, self.state)
            for p in filter(lambda p:p.t<T, self.state):
                q = self.branch(p)
                nextState.extend(q)
                self.omega[len(q)] += 1
            self.state = nextState
    
    def particleNumber(self):
        u = random.random()
        for i in range(len(self.probaSum)):
            if self.probaSum[i]>=u:
                return i
            
    def defaultTime(self):
        return -log(1.0-random.random())/self.intensity
   

Le payoff de la CVA est implémenté pour le cas du replacement closeout :

In [3]:
from operator import mul

class ReplacementCloseOutPayoff:
    
    def __init__(self, psi, a, p):
        self.psi = psi
        self.a = a
        self.p = p

    def payoff(self, s):        
        p0 = reduce(mul, [self.psi(p.value) for p in s.state])
        p1 = 1.0
        for i in range(len(self.a)):
            if self.p[i]>0.0:
                p1 *= pow(self.a[i]/self.p[i], s.omega[i])
        return p0*p1      

La classe suivante calcule par MonteCarlo la moyenne du payoff à maturité (équation )

In [4]:
from math import sqrt
 
class MonteCarlo:
    
    def __init__(self, system, payoff):
        self.system = system 
        self.payoff = payoff    
            
    def compute(self, N, T):
        mean = 0.0
        square = 0.0
        for i in range(N):
            self.system.simulate(T)
            if (len(self.system.state)>0):
                tmp = self.payoff.payoff(self.system)
                mean += tmp
                square += tmp*tmp                   
        mean /= N
        square /= N                
        return (mean, sqrt(square-mean*mean))

Application : digitales dans le modèle de Black&Scholes

Nous avons fait les tests uniquement pour la résolution de \(EDP1\). Nous procédons aux mêmes tests que ceux faits dans . On prend 2 polynômes \(F_1(u)=\frac{1}{2}\left(u^3-u^2\right)\) et \(F_2(u)=\frac{1}{3}\left(u^3-u^2-u^4\right)\), un payoff \(\psi={1}_{x>1}\) et on regarde la résolution Monte Carlo pour \(\mathcal{L}\) le générateur d'Itô d'un mouvement Brownien géométrique de volatilité \(\sigma = 0.2\) et de valeur de départ \(X_t=1\), le spread de CDS \(S^B=500bp\) (intensité du processus de Poisson) et la maturité \(T\)=10ans. nous donne la distribution optimale \((p_0=0,p_1=0,p_2=p_3=\frac{1}{2})\) pour le premier polynôme et \((p_0=0,p_1=0,p_2=p_3=p_4=\frac{1}{3})\) pour le second (distributions uniformes). On a de plus \(||\psi||_\infty=1\). On peut vérifier que les conditions de convergence de l'algorithme sont satisfaites dans les 2 cas.

In [5]:
from math import exp

class BlackScholes:
    def __init__(self, sigma, r):
        self.sigma = sigma
        self.r = r
        
    def next(self, s, deltaT):
        deltaW = random.normalvariate(0.0, sqrt(deltaT))
        return s*exp((self.r-0.5*self.sigma*self.sigma)*deltaT + self.sigma*deltaW)
In [6]:
from ipy_table import apply_theme, make_table
from scipy.stats import norm

# Data
## BlackScholes start value
s0 = 1
## Black volatility
sigma = 0.2
## Riskfree rate
r = 0.0
## Branching probability
p = [0.0,0.0,0.5,0.5]
## Approximation polynomial
a = [0.0,0.0,-0.5,0.5]
## Default intensity
intensity = 0.05
## Payoff
def Psi(x):
    if x>1.0:
        return 1.0
    else:
        return 0.0
## CVA maturity
T = 10.0
## MonteCarlo sample size
N = 2**18
## Confidence level
level = 0.99

# Objects
bs = BlackScholes(sigma, r)
payoff = ReplacementCloseOutPayoff(Psi, a, p)
system = ParticleSystem(s0, bs, intensity, p)
mc = MonteCarlo(system, payoff)


titles = []
titles.append("Valeur de N")
res = []
res.append("Résultat")
conf = []
conf.append("Intervalle de confiance")

for i in range(12, 20, 2):
    N = 2**i
    v = mc.compute(N,T)
    res.append(v[0])
    titles.append("$2^{%i}$"%i)
    conf.append(v[1] * norm.isf(1-0.5*(1+level)) / sqrt(N))

print "Résultats de l'implémentation"
toDisplay = [titles, res, conf]
make_table(toDisplay)
apply_theme('basic')

# CVA
#for i in range(1):
#    v = mc.compute(N,T)
#    confidence = v[1] * norm.isf(1-0.5*(1+level)) / sqrt(N) 
#    print "CVA = [ ", v[0]-confidence, " , ", v[0]+confidence, " ]" 
Résultats de l'implémentation

Out[6]:
In [7]:
print "Résultats de l'article"
ref_res = ["Résultat", 0.2078, 0.2225, 0.2197, 0.2190]
ref_conf = ["Intervalle de confiance", 0.0078, 0.0039, 0.0019, 0.0010]
toDisplay = [titles, ref_res, ref_conf]
make_table(toDisplay)
apply_theme('basic')

powered = [12, 14, 16, 18]
del res[0]
del conf[0]
del ref_res[0]
Résultats de l'article

In [8]:
errorbar(powered, res, yerr=conf, label="Implementation")
plot(powered, ref_res, 'ks', label="Article")

legend(loc=4)

title("Comparison");
In [8]:
Valeur de N$2^{12}$$2^{14}$$2^{16}$$2^{18}$
Résultat0.21880.22660.21850.2190
Intervalle de confiance0.02000.00990.00490.0025