Gaussian Loss Model

Project.toml

name = "GaussianFactorModel"
version = "0.1.0"
uuid = "b4865c2b-8080-47c0-a052-c54cf23146cd"
authors = ["Patrick Haener <contact@haenerconsulting.com>"]

[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

GaussianLossModel.jl

module GaussianFactorModel

using Random
using Distributions
using Statistics
using LinearAlgebra
using Optim

# $x_n = \sum_i a_{n,i} \xi_i + b_n \xi_n\qquad \xi_i$: Gaussian variable
#
# $a_i$: factor loading
#
# $b_n$: idiosyncratic loading
#
# $b_n^2 = 1-\sum_i a_i^2$
#
# $\mathbb{E}[\xi_{i}\xi_j]=\delta_{ij}$

# `factorloadings`: 2 dimensional vector $a_{n,i}$
#
# `idiosyncraticloadings`: vector $b_n$

struct GaussianFactorModelParameters
    factorloadings::Matrix{Float64}
    idiosyncraticloadings::Vector{Float64}
end

function gaussianfactormodelparameters(factorloadings::Matrix{Float64})
    @assert all(x->(x>=-1.0 && x<=1.0), factorloadings) "Illegal value for factor loadings: all loadings must be between -1 and 1"
    i = 1.0 .- mapslices(x->sum(x .* x), factorloadings, dims=[2])
    @assert all(x->(x>=0.0), i) "Illegal value for idiosyncratic loadings: implied imaginary value of idiosyncratic loading"
    GaussianFactorModelParameters(factorloadings, vec(sqrt.(i)))
end

# Single factor model:

function gaussianfactormodelparameters(factorloadings::Vector{Float64})
    gaussianfactormodelparameters(reshape(factorloadings, length(factorloadings), 1))
end

# Single factor model with common driver loading ρ

function gaussianfactormodelparameters(ρ::Float64, numbervariables::Int)
    gaussianfactormodelparameters(fill(ρ, (numbervariables, 1)))
end

# Correlation matrix from `GaussianFactorModelParameter`s

function cor(p::GaussianFactorModelParameters)
    c = p.factorloadings * transpose(p.factorloadings)
    for i in 1:size(p.factorloadings)[1]
        c[i, i] = c[i, i] + p.idiosyncraticloadings[i] * p.idiosyncraticloadings[i]
    end
    c
end

# See https://github.com/JuliaLang/julia/issues/22279
function Base.IteratorSize(p::GaussianFactorModelParameters)
    Base.IsInfinite()
end

# +
function Base.iterate(p::GaussianFactorModelParameters)
    iterate(p, (Random.seed!(), Normal()))
end

function Base.iterate(p::GaussianFactorModelParameters, state::Tuple{Random.AbstractRNG, UnivariateDistribution})
    numberidiosyncratic, numberfactors = size(p.factorloadings)
    factors = rand(state[1], state[2], numberfactors)
    idiosyncratic = rand(state[1], state[2], numberidiosyncratic)
    r = p.factorloadings * factors .+ p.idiosyncraticloadings .* idiosyncratic
   (r, state)
end

end

LossModel

Project.toml

name = "LossModel"
uuid = "11c3f125-d377-4068-8121-5488f83f8688"
authors = ["Patrick Haener <contact@haenerconsulting.com>"]
version = "0.1.0"

[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
GaussianFactorModel = "b4865c2b-8080-47c0-a052-c54cf23146cd"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

LossModel.jl

module LossModel

using Random
using Statistics
using Distributions
using LinearAlgebra

using GaussianFactorModel

# ## Model Parameters

abstract type LossModelParameters end

# ### Simple Model

# For uncorrelated losses with deterministic PD and loss at default

struct SimpleLossModelParameters <: LossModelParameters
    # probability of default
    pds::Vector{Float64}
    # loss when default
    losses::Vector{Float64}
end

function simplelossmodelparameters(pdsandls::Vector{Tuple{Float64, Float64}})
    pds = [x[1] for x in pdsandls]
    ls = [x[2] for x in pdsandls]
    SimpleLossModelParameters(pds, ls) 
end

function simplelossmodelparameters(pds::Vector{Float64}, losses::Vector{Float64})
     @assert length(pds)==length(losses)
    SimpleLossModelParameters(pds, losses)
end

# `SimpleLlosParameters` for `n` entities, each having the same  PD and loss

function simplelossmodelparameters(pd::Float64, loss::Float64, n)
    SimpleLossModelParameters(fill(pd, n), fill(loss, n))
end

# ### Parameters for Gaussian Factor Loss Model

# Deterministic PDs and losses; losses are dependent given by a Gaussian factor model

struct GaussianFactorLossModelParameters <: LossModelParameters
    pds::Vector{Float64}
    losses::Vector{Float64}
    factormodelparameters::GaussianFactorModel.GaussianFactorModelParameters
end

function gaussianfactorlossmodelparameters(pdslosses::Vector{Tuple{Float64, Float64}}, factormodelparameters::GaussianFactorModel.GaussianFactorModelParameters)
    @assert length(pdslosses)==length(factormodelparameters.idiosyncraticloadings) "Parameters imply inconsistent number of variables"
    GaussianFactorLossModelParameters([x[1] for x in pdslosses], [x[2] for x in pdslosses], factormodelparameters)
end

function gaussianfactorlossmodelparameters(pds::Vector{Float64}, losses::Vector{Float64}, factormodelparameters::GaussianFactorModel.GaussianFactorModelParameters)
    @assert length(pds)==length(losses) "Vectors for PDs and for losses must have same length"
    gaussianfactorlossmodelparameters(collect(zip(pds, losses)), factormodelparameters)
end

function gaussianfactorlossmodelparameters(lossmodelparameters::SimpleLossModelParameters, factormodelparameters::GaussianFactorModel.GaussianFactorModelParameters)
    gaussianfactorlossmodelparameters(lossmodelparameters.pds, lossmodelparameters.losses, factormodelparameters)
end

# ## Loss Models

# Returns a vector where the elements of the vector are the losses for a given counterparty for sample 1 to sample `numbersamples`:
#
# $v = [\ell_1,\ldots \ell_n,\ldots,\ell_N]$, $N$ number of counteparties and $\ell_n$ vector containing losses for counterparty $n$ for all samples

function lossespercounterparty(p::M, numbersamples::Int64) where { M<: LossModelParameters }
    numberitems = length(p.pds)
    result = [zeros(Float64, numbersamples) for n in 1:numberitems]
    r = Iterators.take(p, numbersamples)
    for (n, l) in Iterators.enumerate(r)
        for i in 1:numberitems
            result[i][n] = l[i]
        end
    end
    result
end

# Returns a vector where the elements of the vector are the losses for a given sample for each counterparty:
#
# $v = [\ell_1,\ldots \ell_i,\ldots,\ell_{\text{numbersamples}}]$,  $\ell_i$ vector containing losses for sample $i$ for all counterparties

function lossespersample(p::M, numbersamples::Int64) where { M<: LossModelParameters }
    collect(Vector{Float64}, Iterators.take(p, numbersamples))
end

# ### Simple

# +
# See https://github.com/JuliaLang/julia/issues/22279
function Base.IteratorSize(p::SimpleLossModelParameters)
    Base.IsInfinite()
end

function Base.iterate(p::SimpleLossModelParameters)
    iterate(p, Random.seed!())
end

function Base.iterate(p::SimpleLossModelParameters, state::Random.AbstractRNG)
    q = rand(length(p.pds))
    r = map(x -> x[1] <= x[2] ? x[3] : 0.0, zip(q, p.pds, p.losses))
   (r, state)
end

# Mean of losses

function Statistics.mean(p::SimpleLossModelParameters)
 p.pds .* p.losses
end

# Variance of losses

function var(p::SimpleLossModelParameters)
 p.pds .* p.losses .* p.losses - mean(p) .* mean(p)
end

# Covariance of losses

function cov(p::SimpleLossModelParameters)
    diagm(var(p))
end

# ### Gaussian

# +
# See https://github.com/JuliaLang/julia/issues/22279
function Base.IteratorSize(p::GaussianFactorLossModelParameters)
    Base.IsInfinite()
end

function Base.iterate(p::GaussianFactorLossModelParameters)
    iterate(p, (Random.seed!(), Normal()))
end

function Base.iterate(p::GaussianFactorLossModelParameters, state::Tuple{Random.AbstractRNG, UnivariateDistribution})
    q = Base.iterate(p.factormodelparameters, state)
    z = [quantile(Distributions.Normal(), x) for x in p.pds]
    r = map(x -> x[1] <= x[2] ? x[3] : 0.0, zip(q[1], z, p.losses))
   (r, state)
end
# -

SimpleLossModelParameters

function Statistics.mean(p::GaussianFactorLossModelParameters)
    p.pds .* p.losses
end

function var(p::GaussianFactorLossModelParameters)
    p.pds .* p.losses .* p.losses - mean(p) .* mean(p)
end

function cov(p::GaussianFactorLossModelParameters)
    throw("Not implemented yet. Need to calcualte CDF(z1, z2)")
end

# ## Loss Aggregation Functions

# $f([L_1,\ldots,L_i,\ldots, L_{\text{numbersamples}}])$ where $L_i$ is the total loss for sample $i$

function eulerdecomposition(p::P, f, numbersamples::Int64, ϵ::Float64=0.0001) where { P<: LossModelParameters }
    x = lossespercounterparty(p, numbersamples)
    value = f(sum(x))
    decomposition = contributions(x, value, f, numbersamples, ϵ)
    (value, decomposition)
end

# Losses is a vector $[\ell_1,\ldots \ell_n,\ldots,\ell_N]$, $N$ number of counteparties and $\ell_n$ vector containing losses for counterparty $n$ for all samples
#
# The output of function `lossespercounterparty` yield such a vector.
#
# `i` is the index of the counterparty: i=1,..., (number of counterparties)
#

function contribution(losses::Vector{Vector{Float64}}, value::Float64, i::Int64, f, numbersamples::Int64, ϵ::Float64=0.0001) where { P<: LossModelParameters }
    x = copy(losses)
    x_old = x[i]
    x_new = x_old .* (1.0 + ϵ)
    x[i] = x_new
    b = f(sum(x))
    (b-value) / ϵ 
end

function contributions(losses::Vector{Vector{Float64}}, value::Float64, f, numbersamples::Int64, ϵ::Float64=0.0001) where { P<: LossModelParameters }
    [contribution(losses, value, i , f, numbersamples, ϵ) for i in 1:length(losses)]
end

end

Vasicek Loss Model

Project.toml

name = "VasicekLossModel"
uuid = "0f8c25db-2d65-45d2-a5d9-c1b381ca065a"
authors = ["Patrick Haener <contact@haenerconsulting.com>"]
version = "0.1.0"

[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

VasicekLossModel.jl

module VasicekLossModel

using Distributions

const normal = Normal()

function vasiceklossquantile(p::Float64, ρ::Float64, q::Float64)
    cdf(normal, (quantile(normal, p) + sqrt(ρ)*quantile(normal, q))/sqrt(1-ρ))
end

# Fraction of defaulted entities for finite number of entities $n>>1$ and $\rho=0$

function fractiondefaulteduncorrelatedapprox(p::Float64, n::Int, q::Float64)
    quantile(Normal(p, sqrt(p*(1-p)/n)), q)
end

end

LossData

Project.toml

name = "LossData"
uuid = "011f8915-1e18-4f05-bb6b-684427991dcf"
authors = ["Patrick Haener <contact@haenerconsulting.com>"]
version = "0.1.0"

[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
GaussianFactorModel = "b4865c2b-8080-47c0-a052-c54cf23146cd"
LossModel = "11c3f125-d377-4068-8121-5488f83f8688"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
VasicekLossModel = "0f8c25db-2d65-45d2-a5d9-c1b381ca065a"

LossData.jl

module LossData 

using Random
using Distributions
using DataFrames
using GaussianFactorModel
using LossModel

abstract type PartitioningIndex end

struct NettingSetIdentifier{T} <: PartitioningIndex
  id::T
end

function string(id::NettingSetIdentifier{T}) where {T}
    id.id
end

@enum Region EW EAP EECA UK NA SALAC MENA SSA SA

const regionnames = Dict([
        (EW, "Europe (west) area"),
        (EAP, "East Asia & Pacific"),
        (EECA, "Eastern Europe & Central Asia"),
        (UK, "United Kingdom"),
        (NA, "North America"),
        (SALAC, "South America, Latin America & Carribean"),
        (MENA, "Middle East & North Africa"), 
        (SSA, "Sub-Saharan Africa"),
        (SA, "South Asia")
    ])

function string(r::Region)
    get(regionnames, r, Nothing)
end

function region(s::String)
    d = Dict([(get(regionnames, r, Nothing), r) for r in instances(Region)])
    get(d, s, Nothing)
end


@enum Sector FI SR M MQ RW BSO TUS RE

const sectornames = Dict([
    (FI, "Finance Industry"),
    (SR, "Sovereigns/Regionals"),
    (M, "Manufacturing"),
    (MQ, "Mining & Quarrying"),
    (RW, "Retail / Wholesale trade"),
    (BSO, "Buisiness Service & Other"),
    (TUS, "Transport, Utilities & Storage"),
    (RE, "Real Estate")
])

function string(r::Sector)
    get(sectornames, r, Nothing)
end

function sector(s::String)
    d = Dict([(get(sectornames, q, Nothing), q) for q in instances(Sector)])
    get(d, s, Nothing)
end

struct RegionSectorPartitioningIndex <: PartitioningIndex
  region::Region
  sector::Sector
end

function instances(::Type{RegionSectorPartitioningIndex})
  [RegionSectorPartitioningIndex(r, s) for r in instances(Region) for s in instances(Sector)]
end

function nettingsetidentifiers(n::I) where {I<:Integer}
   [NettingSetIdentifier(i) for i in 1:n] 
end

function randomnotionals(nettingsetidentifiers::Vector{NettingSetIdentifier{T}}, maxnotional::Float64) where {T}
    n = rand(Uniform(0.0, maxnotional), length(nettingsetidentifiers))
    collect(zip(nettingsetidentifiers, n))
end

function randompds(nettingsetidentifiers::Vector{NettingSetIdentifier{T}}, maxpd::Float64) where {T}
    p = rand(Uniform(0.0, maxpd), length(nettingsetidentifiers))
    collect(zip(nettingsetidentifiers, p))
end

function randomsimplelossmodelparameters(nettingsetidentifiers::Vector{NettingSetIdentifier{T}},  maxpd::Float64, maxloss::Float64) where {T}
    l = rand(Uniform(0.0, maxloss), length(nettingsetidentifiers))
    p = rand(Uniform(0.0, maxpd), length(nettingsetidentifiers))
    (nettingsetidentifiers, LossModel.SimpleLossModelParameters(l,p))
end

function randomsimplelossmodelparameters(n::Int64,  minmaxpd::Tuple{Float64,Float64}, minmaxloss::Tuple{Float64, Float64})
    p = rand(Uniform(minmaxpd[1], minmaxpd[2]), n)
    l = rand(Uniform(minmaxloss[1], minmaxloss[2]), n)
    LossModel.SimpleLossModelParameters(p, l)
end

function randomsimplelossmodelparameters(n::Int64,  maxpd::Float64, maxloss::Float64)
    randomsimplelossmodelparameters(n, (0.0, maxpd), (0.0, maxloss))
end

function randomlossmodelparameters(nettingsetidentifiers::Vector{NettingSetIdentifier{T}}, maxnotional::Float64, maxpd::Float64) where {T}
    n = rand(Uniform(0.0, maxnotional), length(nettingsetidentifiers))
    p = rand(Uniform(0.0, maxpd), length(nettingsetidentifiers))
    collect(zip(nettingsetidentifiers, n, p))
end

function randomgaussianfactorlossmodelparameters(lossmodelparameters::LossModel.SimpleLossModelParameters, numberdrivers::Int64, minmaxloading::Tuple{Float64,Float64})
    a = rand(Uniform(minmaxloading[1], minmaxloading[2]), (length(lossmodelparameters.pds), numberdrivers))
    g = GaussianFactorModel.gaussianfactormodelparameters(a)
    LossModel.GaussianFactorLossModelParameters(lossmodelparameters.pds, lossmodelparameters.losses, g)
end

function convert(::Type{DataFrame}, nettingsetidentifiers::Vector{NettingSetIdentifier{T}}, p::LossModel.GaussianFactorLossModelParameters) where {T}
    @assert length(nettingsetidentifiers) == length(p.losses) "Number of netting set identifiers not identical to the number of risk losses in the loss model"
    ndrivers = size(p.factormodelparameters.factorloadings)[2]
    drivernames = map(x->"mean$(x-1)", 1:ndrivers)
    df = hcat(DataFrame(GID = map(string, nettingsetidentifiers), EAD = p.losses, PD = p.pds), DataFrame(p.factormodelparameters.factorloadings, drivernames))
    df
end

function convert(::Type{LossModel.GaussianFactorLossModelParameters}, df::DataFrame)
    nettingsets = df[:, "GID"]
    eads = df[:, "EAD"]
    pds = df[:, "PD"]
    numberdrivers = 5
    drivernames = map(x->"mean$(x-1)", 1:numberdrivers)
    loadings = Matrix{Float64}(df[:, drivernames])
    g = GaussianFactorModel.gaussianfactormodelparameters(loadings)
    LossModel.GaussianFactorLossModelParameters(pds, eads, g)
end

end