1 Stressing Portfolios Using Wasserstein Distance

We consider the problem where a modeller conducts sensitivity analysis of a model consisting of random input factors, a corresponding random output of interest, and a baseline probability measure. The modeller seeks to understand how the model (the distribution of the input factors as well as the output) changes under a stress on the output’s distribution. Specifically, for a stress on the output random variable, we derive the unique stressed distribution of the output that is closest in the Wasserstein distance to the baseline output’s distribution and satisfies the stress. We further derive the stressed model, including the stressed distribution of the inputs, which can be calculated in a numerically efficient way from a set of baseline Monte Carlo samples.

The proposed reverse sensitivity analysis framework is model-free and allows for stresses on the output such as (a) the mean and variance, (b) any distortion risk measure including the Value-at-Risk and Expected-Shortfall, and (c) expected utility type constraints, thus making the reverse sensitivity analysis framework suitable for risk models.

The Wasserstein SWIM package can be found at: https://github.com/judymao/SWIM/tree/wasserstein.

This code is based on the paper: https://privpapers.ssrn.com/sol3/papers.cfm?abstract_id=3878879.

For an reverse sensitivity analysis in Python, see: https://github.com/spesenti/SWIM-Wasserstein.

For sensitivity analysis using relative entropy, see: https://github.com/spesenti/SWIM.

1.1 Introductory Example Dataset

We will be using the introductory dataset to illustrate the stresses using Wasserstein distances. For details on the dataset, see Section 2.2 in Pesenti et al. (2021) Scenario Weights for Importance Measurement (SWIM) – an R package for sensitivity analysis, Annals of Actuarial Science, 15(2), 458-483. doi:10.1017/S1748499521000130. We apply different stresses to the portfolio loss \(Y\) with distribution function \(G_Y\) and optimise for the distribution function \(G^*_Y\) that satisfies the stress and is closest to \(G_Y\) in the 2-Wasserstein distance.

# dataset
set.seed(0)
# number of simulated scenarios
n.sim <- 10 ^ 5
# correlation between Z1 and Z2
r <- 0.5
# simulation of Z1  and Z2
# constructed as a combination of independent standard normals U1, U2
U1 <- rnorm(n.sim)
U2 <- rnorm(n.sim)
Z1 <- 100 + 40 * U1
Z2 <- 100 + 20 * (r * U1 + sqrt(1 - r ^ 2) * U2)
# simulation of Z3
Z3 <- rnorm(n.sim, 100, 20)
# portfolio loss Y
Y <- Z1 + Z2 + Z3
# data of baseline model
dat <- data.frame(Z1, Z2, Z3, Y)

1.2 Stressing the Distortion Risk Measure

First, we stress the distortion risk measure of \(Y\). The distortion risk measure \(\rho_{\gamma}\) with distortion weight function \(\gamma\) is defined as

\[\begin{equation*} \rho_{\gamma}(G_Y) = \int_0^1 \breve{G}_Y(u) \gamma(u) du, \end{equation*}\]

where \(G_Y\) is a distribution function for \(Y\), \(\gamma \in \mathbb{L}^2([0,1])\) is a square-integrable function with \(\gamma: [0,1] \to [0,1]\) and \(\int_0^1 \gamma(u) du = 1\).

We will be using the Expected Shortfall (ES) risk measure, which is defined by \[\begin{equation*} \gamma(u) = \frac{1}{1 - \alpha}\mathbb{1}_{\{u \geq \alpha\}}. \end{equation*}\]

We apply a stress of \(10\%\) to the ES risk measure at level \(\alpha=0.95\), where the stress is a relative change in the corresponding values under the baseline model.

stress.es <- stress_wass(type = "RM", x = dat, q_ratio=c(1.1), alpha=c(0.95), k = 'Y')
summary(stress.es, base = TRUE)
## $base
##                        Z1            Z2            Z3             Y
## mean         1.000564e+02  99.940396075  99.984331513 299.981113486
## sd           4.004115e+01  19.996852012  19.981845455  56.638588761
## skewness    -6.084411e-04   0.001170933  -0.002470021  -0.002340723
## ex kurtosis -1.056073e-02  -0.008906052  -0.012511967  -0.009319609
## 1st Qu.      7.283343e+01  86.474503862  86.481609311 261.612122776
## Median       1.001990e+02  99.986573249 100.009073888 300.054804363
## 3rd Qu.      1.271558e+02 113.395668559 113.493418432 338.266983662
## 
## $`stress 1`
##                      Z1          Z2            Z3           Y
## mean        101.2144998 100.3976315 100.211102032 301.8232334
## sd           42.4764660  20.7873489  20.151509321  61.0244260
## skewness      0.2944971   0.1670657   0.015275539   0.4131771
## ex kurtosis   0.4320920   0.2137675  -0.003681972   0.7162995
## 1st Qu.      72.7906607  86.4627613  86.601090177 261.5441663
## Median      100.0977946  99.9916661 100.143470157 299.9136335
## 3rd Qu.     127.0971538 113.6797403 113.812475568 337.9436966

We can see the distributional characteristics using the summary function. We can also view how the probability distributions have changed for both the inputs and output by plotting them.

plot_cdf(stress.es, xCol = "Z1", base = TRUE)
plot_cdf(stress.es, xCol = "Y", base = TRUE)
Figure 1: Stressed and baseline distributions of $Z_1$ (left) and stressed and baseline distributions of $Y$ (right).Figure 1: Stressed and baseline distributions of $Z_1$ (left) and stressed and baseline distributions of $Y$ (right).

Figure 1: Stressed and baseline distributions of \(Z_1\) (left) and stressed and baseline distributions of \(Y\) (right).

We see a heavier right tail on both the stressed distributions of \(Z_1\) and \(Y\), which is necessary to meet the stress on the \(ES_{0.95}\).

We also plot the scenario weights and sensitivity measure to analyze how the distributions of the factors have been affected by the stress. We use the reverse sensitivity measure, where the reverse sensitivity measure is defined as:

\[\begin{equation*} S_i^{\mathbb{Q}^*} = \begin{cases} \frac{\mathbb{E}^{\mathbb{Q}^*}[s(Z_i)] - \mathbb{E}[s(Z_i)]}{\max_{\mathbb{Q} \in Q}\mathbb{E}^{\mathbb{Q}}[s(Z_i)] - \mathbb{E}[s(Z_i)]} \quad \quad & \mathbb{E}^{\mathbb{Q}^*}[s(Z_i)] \geq \mathbb{E}[s(Z_i)] \\ -\frac{\mathbb{E}^{\mathbb{Q}^*}[s(Z_i)] - \mathbb{E}[s(Z_i)]}{\min_{\mathbb{Q} \in Q}\mathbb{E}^{\mathbb{Q}}[s(Z_i)] - \mathbb{E}[s(Z_i)]} & \textrm{otherwise,} \end{cases} \end{equation*}\]

where \(s:\mathbb{R} \rightarrow \mathbb{R}\) is some function, \(Q =\) {\(\mathbb{Q}\) | \(\mathbb{Q}\) probability measure with \(\frac{d\mathbb{Q}}{d\mathbb{P}} \stackrel{\mathbb{P}}{=} \frac{d\mathbb{Q}^*}{d\mathbb{P}}\)} is the set of all probability measures whose RN-derivative have the same distribution as \(\frac{d\mathbb{Q}^*}{d\mathbb{P}}\) under \(\mathbb{P}\).

We use \(s(Z_i) = Z_i\).

plot_weights(stress.es, xCol='Y')
plot_sensitivity(stress.es, type="reverse", xCol=c(1, 2, 3))
Figure 2: Scenario weights for $Y$ (left) and reverse sensitivity measure (right).Figure 2: Scenario weights for $Y$ (left) and reverse sensitivity measure (right).

Figure 2: Scenario weights for \(Y\) (left) and reverse sensitivity measure (right).

From the scenario weights, we see that very high weights are placed on higher values of \(Y\) to meet the stressed risk measure. From the sensitivity plot, we can see that \(Z_3\) is the least affected by the stress, which is induced by the fact that \(Z_3\) is independent of \((Z_1, Z_2)\).

1.3 Stressing the Mean and Standard Deviation

Then, we stress the mean and standard deviation of \(Y\). We apply a \(2\%\) stress to the mean and a \(-2\%\) stress to the standard deviation.

curr_mean <- mean(dat[, 'Y'])
curr_sd <- sd(dat[, 'Y'])

stress.mean_sd <- stress_wass(type = "mean sd", x = dat, new_mean=curr_mean*1.02, new_sd=curr_sd*0.98, k = 'Y')
summary(stress.mean_sd)
## $`stress 1`
##                        Z1            Z2            Z3             Y
## mean        103.806246963  1.014389e+02 100.729468692 305.974613453
## sd           39.173922052  1.972443e+01  19.911941399  55.064642864
## skewness     -0.001024437 -7.402582e-04  -0.002933337  -0.002622802
## ex kurtosis  -0.037263442 -1.431740e-02  -0.010979468  -0.059234681
## 1st Qu.      77.214893501  8.813963e+01  87.289287785 268.696551005
## Median      103.888056399  1.014556e+02 100.732819210 305.991493477
## 3rd Qu.     130.366797732  1.146836e+02 114.162611869 343.310412087

We can see the distributional characteristics using the summary function. We can also view how the probability distributions have changed for both the inputs and output by plotting the histograms.

plot_hist(stress.mean_sd, xCol = "Z1", base = TRUE)
plot_hist(stress.mean_sd, xCol = "Y", base = TRUE)
Figure 3: Stressed and baseline histograms of $Z_1$ (left) and stressed and baseline histograms of $Y$ (right).Figure 3: Stressed and baseline histograms of $Z_1$ (left) and stressed and baseline histograms of $Y$ (right).

Figure 3: Stressed and baseline histograms of \(Z_1\) (left) and stressed and baseline histograms of \(Y\) (right).

We can see that the distribution has shifted to the right to match the new mean, as well as become less spread, reflecting the decreased standard deviation.

We then plot the scenario weights and sensitivity measures using the reverse sensitivity measure, with \(s(Z_i) = Z_i\).

plot_weights(stress.mean_sd, xCol='Y')
plot_sensitivity(stress.mean_sd, type="reverse", xCol=c(1, 2, 3))
Figure 4: Scenario weights for $Y$ (left) and reverse sensitivity measure (right).Figure 4: Scenario weights for $Y$ (left) and reverse sensitivity measure (right).

Figure 4: Scenario weights for \(Y\) (left) and reverse sensitivity measure (right).

From the scenario weights, we see the weights are less than 1 for smaller values of \(Y\), and greater than 1 for larger values of \(Y\), to meet the stressed mean. From the sensitivity plot, we can see that \(Z_3\) is the least affected by the stress, as \(Z_3\) is independent of \((Z_1, Z_2)\).

1.4 Stressing the ES, Mean and Standard Deviation

Next, we stress the ES, mean, and standard deviation of \(Y\). We apply a \(5\%\) stress to the \(ES_{0.95}\), a \(5\%\) stress to the mean and a \(-5\%\) stress to the standard deviation.

curr_mean <- mean(dat[, 'Y'])
curr_sd <- sd(dat[, 'Y'])

stress.ES_mean_sd <- stress_wass(type='RM mean sd', x=dat, alpha=0.95, q_ratio = 1.05,
                                 new_mean=curr_mean*1.05, new_sd=curr_sd*0.95, k='Y')
summary(stress.ES_mean_sd)
## $`stress 1`
##                      Z1           Z2            Z3           Y
## mean        109.3437532 103.63158452 101.822457430 314.7977952
## sd           38.2565125  19.41281516  19.831362519  53.3577119
## skewness      0.1120712   0.05165827   0.002893368   0.1686804
## ex kurtosis   0.1235913   0.05569387  -0.005093122   0.2186537
## 1st Qu.      83.6794186  90.53547494  88.461080120 278.9437542
## Median      108.8887491 103.52829233 101.785570073 314.0700589
## 3rd Qu.     134.0733260 116.45928511 115.246109120 348.9563313

We can see the distributional characteristics using the summary function. We can also view how the probability distributions have changed for both the inputs and output by plotting them.

plot_cdf(stress.ES_mean_sd, xCol = "Z1", base = TRUE)
plot_cdf(stress.ES_mean_sd, xCol = "Y", base = TRUE)
Figure 5: Stressed and baseline distributions of $Z_1$ (left) and stressed and baseline distributions of $Y$ (right).Figure 5: Stressed and baseline distributions of $Z_1$ (left) and stressed and baseline distributions of $Y$ (right).

Figure 5: Stressed and baseline distributions of \(Z_1\) (left) and stressed and baseline distributions of \(Y\) (right).

We can see from the distributions that it resembles a combination of the separate stresses on the ES and the mean/standard deviation, with a rightward shift of the distribution and heavy right tail.

We then plot the scenario weights and sensitivity measures using the reverse sensitivity measure, with \(s(Z_i) = Z_i\).

plot_weights(stress.ES_mean_sd, xCol = "Y")
plot_sensitivity(stress.ES_mean_sd, type="reverse", xCol=c(1, 2, 3))
Figure 6: Scenario weights for $Y$ (left) and reverse sensitivity measure (right).Figure 6: Scenario weights for $Y$ (left) and reverse sensitivity measure (right).

Figure 6: Scenario weights for \(Y\) (left) and reverse sensitivity measure (right).

From the scenario weights, we see that the weights increase as \(Y\) increases, with greater weights for the extreme values of \(Y\). From the sensitivity plot, we can see that \(Z_3\) is the least affected by the stress, as \(Z_3\) is independent of \((Z_1, Z_2)\).

1.5 Stressing the HARA Utility and ES

Last, we stress the Hyperbolic absolute risk aversion (HARA) utility function and ES of \(Y\). The HARA utility function is defined by \(u(x) = \frac{1 - \eta}{\eta}(\frac{ax}{1-\eta} + b)^{\eta}\). We choose parameters \(\eta = 0.5\), \(\alpha = 1\) and \(b = 5\), which guarantees concavity.

We apply a \(10\%\) stress to the \(ES_{0.95}\) and a \(10\%\) stress to the HARA utility with parameters \(a=1\), \(b=5\) and \(\eta=0.5\).

stress.hara <- stress_wass(type = "HARA RM", x = dat, a=1, b=5, eta=0.5, q_ratio=1.1,
                           hu_ratio=1.1, alpha=0.95, k = 'Y')
summary(stress.hara)
## $`stress 1`
##                       Z1           Z2           Z3           Y
## mean        137.13193361 114.79097972 107.34857135 359.2714847
## sd           36.19219997  18.80122727  19.67318681  49.5626263
## skewness     -0.03860068  -0.01876674  -0.02311969  -0.0667497
## ex kurtosis  -0.08463995  -0.01392352  -0.04373633  -0.1267752
## 1st Qu.     112.61963522 102.18888079  94.16400282 325.4788378
## Median      137.28902824 114.82229230 107.35772965 360.0546934
## 3rd Qu.     161.89549118 127.52229063 120.67892897 394.0697650

We can see the distributional characteristics using the summary function. We can also view how the probability distributions have changed for both the inputs and output by plotting them.

plot_cdf(stress.hara, xCol = "Z1", base = TRUE)
plot_cdf(stress.hara, xCol = "Y", base = TRUE)
Figure 7: Stressed and baseline distributions of $Z_1$ (left) and stressed and baseline distributions of $Y$ (right).Figure 7: Stressed and baseline distributions of $Z_1$ (left) and stressed and baseline distributions of $Y$ (right).

Figure 7: Stressed and baseline distributions of \(Z_1\) (left) and stressed and baseline distributions of \(Y\) (right).

The stressed distribution has shifted to the right, which reflects the stress on the mean induced by a positive stress on the HARA utility. We also see the heavier right tail associated with stressing the risk measure.

We then plot the scenario weights and sensitivity measures using the reverse sensitivity measure, with \(s(Z_i) = Z_i\).

plot_weights(stress.hara, xCol="Y")
plot_sensitivity(stress.hara, type="reverse", xCol=c(1, 2, 3))
## Warning in plot_sensitivity(stress.hara, type = "reverse", xCol = c(1, 2, : No s
## passed in. Using Gamma sensitivity instead.
Figure 8: Scenario weights for $Y$ (left) and reverse sensitivity measure (right).Figure 8: Scenario weights for $Y$ (left) and reverse sensitivity measure (right).

Figure 8: Scenario weights for \(Y\) (left) and reverse sensitivity measure (right).

From the scenario weights, we see that greater weights are placed on greater values of \(Y\), with large weights placed on the largest values of \(Y\). As expected, from the sensitivity plot we can see that \(Z_3\) is the least affected by the stress, as \(Z_3\) is independent of \((Z_1, Z_2)\).

1.6 Comparing Wasserstein Distance and K-L Divergence

We now compare the stressed distribution derived from the Wasserstein distance with that of the Kullback–Leibler (KL) divergence that is implemented in the R Package SWIM. We compare stresses on the mean and standard deviation, applied to \(Z_1\). We set the stressed mean to 100, and the stressed standard deviation to 50.

stress.mean_sd_wass <- stress_wass(type = "mean sd", x = dat, new_mean=100, new_sd=50, k = 1)
stress.mean_sd_kl <- stress(type = "mean sd", x = dat, k = 1, new_means=100, new_sd = 50)
##   cols required_moment achieved_moment     abs_error     rel_error
## 1    1             100             100 -6.394885e-13 -6.394885e-15
## 2    1           12500           12500 -1.109584e-10 -8.876668e-15

We then compare the resulting distributions.

plot_cdf(stress.mean_sd_wass, xCol = "Z1", base = TRUE)
plot_cdf(stress.mean_sd_kl, xCol = "Z1", base = TRUE)
## Registered S3 method overwritten by 'spatstat.geom':
##   method     from
##   print.boxx cli
Figure 9: Stressed and baseline distributions of $Z_1$ using Wasserstein Distance (left) and stressed and baseline distributions of $Z_1$ using K-L Divergence (right).Figure 9: Stressed and baseline distributions of $Z_1$ using Wasserstein Distance (left) and stressed and baseline distributions of $Z_1$ using K-L Divergence (right).

Figure 9: Stressed and baseline distributions of \(Z_1\) using Wasserstein Distance (left) and stressed and baseline distributions of \(Z_1\) using K-L Divergence (right).

From the distribution plots, the results are very similar, as expected since both methods must meet the same stresses, thus resulting in very similar stressed distributions.

We then compare how the scenario weights are distributed for both methods.

plot_weights(stress.mean_sd_wass, xCol = "Z1")
plot_weights(stress.mean_sd_kl, xCol = "Z1")
Figure 10: Scenario weights for $Z_1$ using Wasserstein Distance (left) and scenario weights for $Z_1$ using K-L Divergence (right).Figure 10: Scenario weights for $Z_1$ using Wasserstein Distance (left) and scenario weights for $Z_1$ using K-L Divergence (right).

Figure 10: Scenario weights for \(Z_1\) using Wasserstein Distance (left) and scenario weights for \(Z_1\) using K-L Divergence (right).

We now see a difference between using the Wasserstein distance and K-L divergence. The weights for the K-L divergence are very symmetrical about the stressed mean 100, whereas the weights when using the Wasserstein distance show greater variety, which might be due to the use of kernel density approximation. Nevertheless, both follow the pattern of smaller weights near \(Z_1=100\) and greater weights towards the extreme values of \(Z_1\).