iPython Cookbook - Monte Carlo Risk Analysis - Large Pool Capital Model

Looking at the Large Pool Capital Model and idiosyncratic risks on top of it


In [1]:
import numpy as np
import matplotlib.pyplot as plt

Model calculator factory

We define a factory function that then provides us with a calculation function for our specific model. The factory function most importantly takes the number of Monte Carlo samples to be drawn N. For convenience it also takes a default payoff function as often the payoff will not be varied over the course of one analysis.

Importantly, the vector of Standard Gaussian random variables z that will be the basis of the Monte Carlo simulation is drawn immediately, and it is remembered in the closure This allows us to run the calculation function (the one returned by the factory) with different parameters but against the same sampling set, thereby greatly reducing Monte Carlo errors.

TODO operator module / partial

In [2]:
from scipy.stats import norm
from functools import partial
from operator import le

def calc_factory (N=1000, d=10, defaul_payoff=None):
    z = np.random.standard_normal((N,d+1))
    def run (pd=0.05, lgd=1, rho=0):
        m = np.vstack((rho * np.ones(d),sqrt(1-rho*rho) * np.identity(d)))
        x = np.dot(z,m)
        ppf_pd = norm.ppf(pd)
        def lossf (x):
            if x < ppf_pd:
                return 1
                return 0
        lossf = partial(le,ppf_pd)
        lossv = list(map(lambda xx: list(map(lossf, xx)), x))
        nloss = list(map(lambda xx: sum(xx), lossv))
        loss  = list(map(lambda xx: xx*lgd, nloss))
        return (loss, nloss)
    return run

In [3]:
calc = calc_factory(N=1000, d=1000)

In [4]:
loss, lossn = calc (pd=0.01, lgd=0.4*1000, rho=0.1)

In [5]:
hgn=np.histogram(lossn, bins=50)
xvalsn = np.delete(hgn[1], -1)
plt.bar(xvalsn,hgn[0], width=0.9*(xvalsn[1]-xvalsn[0]))
plt.title('default distribution')
plt.xlabel('number of defaults')

In [6]:
hg=np.histogram(loss, bins=50)
xvals = np.delete(hg[1], -1)
plt.bar(xvals,hg[0], width=0.9*(xvals[1]-xvals[0]))
plt.title('loss distribution')
plt.xlabel('loss amount')

Licence and version

(c) Stefan Loesch / oditorium 2014; all rights reserved (license)

In [7]:
import sys

3.4.0 (default, Apr 11 2014, 13:05:11) 
[GCC 4.8.2]