These notes are based on the Exam 7 Syllabus reading Enterprise Risk Analysis for Property and Liability Insurance Companies by Brehm, Gluck, Kreps, Major, Mango, Shaw, Venter, White, and Witcraft.
easypackages::packages("data.table", "DT", "ggplot2")
options(scipen = 999)
A crucial component is tail dependency, which is the degree of dependency in extreme events. For example, inflation may be generally uncorrelated across two lines of business, unless inflation overall is extremely high. Later sections show how tail dependency can be modelled using copulas.
A key consideration is that extreme reference points, such as a total loss of capital, are where the models are least reliable; there is little data to support modelling the tail.
An internal risk model (IRM) is a DFA model that produces an aggregate loss distribution for the firm. Once this has been produced, there are three ways in which it is used for decision making.
The ratios are compared to fixed values for all companies, regardless of the risks to which the company is exposed. Furthermore, it only considers underwriting risks.
Examples provided by the Insurance Regulatory Information System (IRIS) include:The IRIS test compares these ratios to a range of reasonable values, and if four or more fall outside the range, then additional regulatory scrutiny is warranted.
A correctly executed IRM will impact risk decision making in planning, pricing, reinsurance purchase, capacity allocation, and interactions with rating agencies. Given the wide-ranging impact of IRM, it requires careful planning. Major steps in the process are listed below.
Although this section appears earlier in the learning objectives, I have placed it after the discussion of risk measures in these notes, since some risk measures are used in the discussion.
number.of.simulations = 10000
mu_1 = 8
sigma_1 = 0.25
mu_2 = 7
sigma_2 = 0.5
tau = 0.25
set.seed(12345) # Fix random seed to allow replicability of results
simulation_number = 1:number.of.simulations
adjustment = rnorm(number.of.simulations, mean = 0, sd = tau)
line_1_loss = rlnorm(number.of.simulations, meanlog = mu_1 + adjustment, sdlog = sigma_1)
line_2_loss = rlnorm(number.of.simulations, meanlog = mu_2 + adjustment, sdlog = sigma_2)
company_loss = line_1_loss + line_2_loss
simulation.results = data.table(simulation_number, adjustment, line_1_loss, line_2_loss, company_loss)
Furthermore, for purposes of assessing loss of capital under each simulation, assume that the company charges premium equal to the expected value in each line, with a 12% profit margin:
line.1.premium = 1.12 * exp(mu_1 + sigma_1^2/2)
line.2.premium = 1.12 * exp(mu_2 + sigma_2^2/2)
premium.collected = line.1.premium + line.2.premium
simulation.results[, impact_on_capital := company_loss - premium.collected]
simulation.results[, line_1_capital := line_1_loss - line.1.premium]
simulation.results[, line_2_capital := line_2_loss - line.2.premium]
datatable(simulation.results)
A key consideration is allocating risk to different business units within the company, which may be needed for various reasons:
In this case, when \(Y = \sum_j X_j\), the co-measure \(r\) is defined as \[ r(X_j) = E[h(X_j)L(Y)|g(Y)] \] By the additivity condition, \[ \rho(Y) = \sum_j r(X_j) \] Note that for any given risk measure, the functions \(L\) and \(h\) are not unique, so this does not uniquely define a decomposition into co-measures.
An alternative approach is to use a marginal method, which quantifies the change in overall company risk due to a small change in a business unit’s volume. If a unit with an above-average ratio of profit to risk increases its volume, then the overall ratio of profit to risk will increase. (This is easiest to achieve when a business unit can uniformly change its volume, e.g. with a quota share treaty.) This approach is most effective when the risk measure is proportional to the market value of the risk. The marginal method can be defined as a directional derivative: \[ r_j(Y) = \lim_{\epsilon \rightarrow 0} \frac{\rho(Y + \epsilon X_j) - \rho(Y)}{\epsilon} \]
If a risk measure is scalable (homogeneous of degree 1), then multiplying a random variable by a scalar multiplies the risk measure by the same amount. This typically occurs for measures that are expressed in units of currency. In this case, a marginal risk attribution method exists, and it is a co-measure as well. One way to calculate this is to apply l’H^{o}pital’s Rule: take the derivative of \(\rho(Y + \epsilon X_j)\) with respect to \(\epsilon\), and set \(\epsilon = 0\).
Tail-based measures assess the risk associated with large losses only. Value at Risk (VaR) is just a percentile of the loss distribution. In other words, if \(F(x)\) is the cumulative distribution function for the loss distribution, then \[ \mathrm{VaR} = F^{-1}(p) \] at probability level \(p\). For example, to calculate VaR at the 99th percentile:
VaR_99 = quantile(simulation.results$impact_on_capital, probs = 0.99)
VaR_99_line_1 = quantile(simulation.results$line_1_capital, probs = 0.99)
VaR_99_line_2 = quantile(simulation.results$line_2_capital, probs = 0.99)
print(paste0("VaR-99 for the company is ", VaR_99))
## [1] "VaR-99 for the company is 4661.1768865033"
print(paste0("VaR-99 for line 1 alone is ", VaR_99_line_1))
## [1] "VaR-99 for line 1 alone is 3433.49978735007"
print(paste0("VaR-99 for line 2 alone is ", VaR_99_line_2))
## [1] "VaR-99 for line 2 alone is 2650.04451047233"
We can express VaR as a sum of co-measures: \[ r(X_j) = E[X_j | F(Y) = p] \] However, this is not a very practical approach since there is exactly one simulation for which \(F(Y) = p\), so this is does not produce a stable result. VaR is scalable, and the directional derivative gives the same result as above.
Tail Value at Risk is the average value of losses above a given probability level: \[ TVaR = E[X | X \geq \mathrm{VaR}] = \frac{\int_p^{\infty} xf(x)dx}{1-p} \] Using the simulated data:
TVaR_99 = simulation.results[impact_on_capital >= VaR_99, mean(impact_on_capital)]
TVaR_99_line_1 = simulation.results[line_1_capital >= VaR_99_line_1, mean(line_1_capital)]
TVaR_99_line_2 = simulation.results[line_2_capital >= VaR_99_line_2, mean(line_2_capital)]
print(paste0("TVaR-99 for the company is ", TVaR_99))
## [1] "TVaR-99 for the company is 5690.81822767097"
print(paste0("TVaR-99 for line 1 alone is ", TVaR_99_line_1))
## [1] "TVaR-99 for line 1 alone is 4332.5867110888"
print(paste0("TVaR-99 for line 2 alone is ", TVaR_99_line_2))
## [1] "TVaR-99 for line 2 alone is 3528.84379408856"
print(paste0("Proportion of risk allocated to line 1 is ", TVaR_99_line_1 / (TVaR_99_line_1 + TVaR_99_line_2)))
## [1] "Proportion of risk allocated to line 1 is 0.55111938065667"
TVaR is generally more useful than VaR, because we don’t see a VaR-99 loss once every 100 years. Instead, once every 100 years, we see an average loss equal to TVaR.
We can decompose TVaR as a sum of the co-measure \[ r(X_j) = E[X_j | Y > \mathrm{VaR}] \] since in this case, \(\mathrm{TVaR} = E[\sum_j X_j | Y > \mathrm{VaR}] = \sum_j r(X_j)\). In order to allocate TVaR among lines of business, one approach is to calculate the average contribution to the loss, conditional on the overall company results exceeding the 99th percentile:
co_TVaR_99_line_1 = simulation.results[impact_on_capital >= VaR_99, mean(line_1_capital)]
print(paste0("co-TVaR-99 for line 1 is ", co_TVaR_99_line_1))
## [1] "co-TVaR-99 for line 1 is 3407.91649867727"
co_TVaR_99_line_2 = simulation.results[impact_on_capital >= VaR_99, mean(line_2_capital)]
print(paste0("co-TVaR-99 for line 2 is ", co_TVaR_99_line_2))
## [1] "co-TVaR-99 for line 2 is 2282.9017289937"
print(paste0("Proportion of risk allocated to line 1 is ", co_TVaR_99_line_1 / TVaR_99))
## [1] "Proportion of risk allocated to line 1 is 0.598844728883917"
Note that the alternate TVaR approach does not apply well in practice to VaR, because there is only one simulation that occurs exactly at the 99th percentile of the overall company performance.
A closely-related metric is Excess Tail Value at Risk (XTVaR), which is equal to TVaR less the mean. This measure reflects the idea that the mean is already financed by other funding, so capital is only needed for losses above the mean:
XTVaR_99 = TVaR_99 - simulation.results[, mean(impact_on_capital)]
XTVaR_99_line_1 = TVaR_99_line_1 - simulation.results[, mean(line_1_capital)]
XTVaR_99_line_2 = TVaR_99_line_2 - simulation.results[, mean(line_2_capital)]
print(paste0("XTVaR-99 for the company is ", XTVaR_99))
## [1] "XTVaR-99 for the company is 6069.59787861799"
print(paste0("XTVaR-99 for line 1 is ", XTVaR_99_line_1))
## [1] "XTVaR-99 for line 1 is 4605.54371673404"
print(paste0("XTVaR-99 for line 2 is ", XTVaR_99_line_2))
## [1] "XTVaR-99 for line 2 is 3634.66643939035"
Similarly to TVaR, we can decompse XTVaR as a sum of the co-measure \[ r(X_j) = E[X_j - E[X_j] | Y > \mathrm{VaR}] \]
mean.line.1 = simulation.results[, mean(line_1_capital)]
mean.line.2 = simulation.results[, mean(line_2_capital)]
co_XTVaR_99_line_1 = simulation.results[impact_on_capital >= VaR_99, mean(line_1_capital - mean.line.1)]
co_XTVaR_99_line_2 = simulation.results[impact_on_capital >= VaR_99, mean(line_2_capital - mean.line.2)]
print(paste0("co-XTVaR-99 for line 1 is ", co_XTVaR_99_line_1))
## [1] "co-XTVaR-99 for line 1 is 3680.87350432251"
print(paste0("co-XTVaR-99 for line 2 is ", co_XTVaR_99_line_2))
## [1] "co-XTVaR-99 for line 2 is 2388.72437429548"
This is also a marginal risk measure, though the directional derivative is not obvious.
Expected Policyholder Deficit (EPD) is the expected value of defaulted losses, if capital available were set at VaR: \[ \mathrm{EPD} = (\mathrm{TVaR} - \mathrm{VaR}) (1-p) \] Note that \(\mathrm{TVaR} - \mathrm{VaR}\) is the conditional expected loss, given the loss exceeds \(\mathrm{VaR}\). Multiplying by the probability of exceeding \(\mathrm{VaR}\) converts this into an unconditional expectation:
EPD_99 = (TVaR_99 - VaR_99) * 0.01
EPD_99_line_1 = (TVaR_99_line_1 - VaR_99_line_1) * 0.01
EPD_99_line_2 = (TVaR_99_line_2 - VaR_99_line_2) * 0.01
print(paste0("The expected policyholder deficit is ", EPD_99))
## [1] "The expected policyholder deficit is 10.2964134116767"
print(paste0("The expected policyholder deficit for line 1 is ", EPD_99_line_1))
## [1] "The expected policyholder deficit for line 1 is 8.99086923738736"
print(paste0("The expected policyholder deficit for line 2 is ", EPD_99_line_2))
## [1] "The expected policyholder deficit for line 2 is 8.78799283616238"
To define a co-measure, let \(B = F_Y^{-1}(p)\) and express EPD as \[ \rho(Y) = (1-p) E[Y-B | Y > B] \] This is only scalable when \(B\) is a fixed quantile of \(Y\). In this case, the marginal decomposition is \[ r(X_j) = \alpha[\mathrm{coTVaR - coVaR}] \] This follows from the definition of EPD, and linearity of the directional derivative.
A simple moment-based measure is the variance of the decrease in capital: Typically, standard deviation is preferred over variance because it is in units of currency rather than currency squared:
sd_capital = simulation.results[, sd(impact_on_capital)]
sd_capital_1 = simulation.results[, sd(line_1_capital)]
sd_capital_2 = simulation.results[, sd(line_2_capital)]
print(paste0("The standard deviation of the change in capital is ", sd_capital, " and the variance is ", sd_capital^2))
## [1] "The standard deviation of the change in capital is 1572.86167913589 and the variance is 2473893.86169418"
print(paste0("The standard deviation of the change in capital for line 1 only is ", sd_capital_1, " and the variance is ", sd_capital_1^2))
## [1] "The standard deviation of the change in capital for line 1 only is 1162.73013844201 and the variance is 1351941.37484137"
print(paste0("The standard deviation of the change in capital for line 2 only is ", sd_capital_2, " and the variance is ", sd_capital_2^2))
## [1] "The standard deviation of the change in capital for line 2 only is 772.313690179503 and the variance is 596468.436038681"
Note that due to the correlation between the lines, we would underestimate the standard deviation if we were to measure each line individually and combine them:
sqrt(sd_capital_1^2 + sd_capital_2^2)
## [1] 1395.855
Since we can express \[ \mathrm(Var)(Y) = E[(Y - E[Y])(Y - E[Y])] \] by setting \(h(Y) = L(Y) = Y - E[Y]\), we obtain the co-variance as a co-measure: \[ \mathrm{Cov}(X_j, Y) = r(X_j) = E[(X_j - E[X_j])(Y - E[Y])] \] In this example:
mean_X_1 = simulation.results[, mean(line_1_capital)]
mean_X_2 = simulation.results[, mean(line_2_capital)]
mean_Y = mean_X_1 + mean_X_2
mean_X_1_Y = simulation.results[, mean(line_1_capital * impact_on_capital)]
mean_X_2_Y = simulation.results[, mean(line_2_capital * impact_on_capital)]
cov_X_1 = mean_X_1_Y - mean_X_1 * mean_Y
cov_X_2 = mean_X_2_Y - mean_X_2 * mean_Y
print(paste0("The covariance of line 1 with the overall company result is ", cov_X_1))
## [1] "The covariance of line 1 with the overall company result is 1614521.93190841"
print(paste0("The covariance of line 2 with the overall company result is ", cov_X_2))
## [1] "The covariance of line 2 with the overall company result is 859124.540399601"
Variance is not scalable, since it is in units of currency squared, so the directional derivative approach does not apply.
For standard deviation, there are two reasonable candidates for a co-measure. The first is \[ r_1(X_j) = E\left[X_j \frac{\mathrm{sd}(Y)}{E[Y]}\right] \] The second is \[ r_2(X_j) = E\left[(X_j - E[X_j]) \frac{Y - E[Y]}{\mathrm{sd}(Y)}\right] = \frac{\mathrm{Cov}(X_j, Y)}{\mathrm{sd}(Y)} \] The marginal risk measure can be determined through differentiation of \[ \rho(Y + \epsilon X_j) = (\mathrm{Var}(Y + \epsilon X_j))^{1/2} \] First, re-express this as \[ \rho(Y + \epsilon X_j) = (\mathrm{Var}(Y) + 2\epsilon\mathrm{Cov}(X_j, Y) + \epsilon^2 \mathrm{Var}(X_j)^2)^{1/2} \] The derivative with respect to \(\epsilon\) is \[ \frac{d\rho(Y + \epsilon X_j)}{d\epsilon} = \frac12 (\mathrm{Var}(Y) + 2\epsilon\mathrm{Cov}(X_j, Y) + \epsilon^2 \mathrm{Var}(X_j)^2)^{-1/2}(2 \mathrm{Cov}(X_j, Y) + 2\epsilon\mathrm{Var}(X_j)^2) \] Setting \(\epsilon = 0\), we obtain \[ r(X_j) = \frac{\mathrm{Cov}(X_j, Y)}{\mathrm{sd}(Y)} \] indicating that the function \(r_2\) defined above is the unique marginal co-measure.
A disadvantage of this approach is that it treats favourable deviations the same as unfavourable ones. The semistandard deviation uses only the unfavourable deviations:
semi_sd_capital = simulation.results[impact_on_capital > 0, sd(impact_on_capital)]
semi_sd_capital_1 = simulation.results[line_1_capital > 0, sd(line_1_capital)]
semi_sd_capital_2 = simulation.results[line_2_capital > 0, sd(line_2_capital)]
print(paste0("The semistandard deviation of the change in capital is ", semi_sd_capital))
## [1] "The semistandard deviation of the change in capital is 1282.5678877533"
print(paste0("The semistandard deviation of the change in capital for line 1 only is ", semi_sd_capital_1))
## [1] "The semistandard deviation of the change in capital for line 1 only is 953.399528035576"
print(paste0("The semistandard deviation of the change in capital for line 2 only is ", semi_sd_capital_2))
## [1] "The semistandard deviation of the change in capital for line 2 only is 762.407926475597"
Quadratic risk measures may not fully reflect market attitudes toward risk. This can be addressed by using higher moments (e.g. skewness), or by using exponential moments: \[ E[Y \exp(cY/E[Y])] \] Setting \(c=0.5\), this value may be calculated as follows:
E_Y = simulation.results[, mean(impact_on_capital)]
simulation.results[, exponential_moment_contribution := impact_on_capital * exp(0.5 * impact_on_capital / E_Y)]
exponential.moment = simulation.results[, mean(exponential_moment_contribution)]
print(paste0("The exponential moment is ", exponential.moment))
## [1] "The exponential moment is -12128.6536034656"
The exponential moment has the advantage that it reflects all losses, but responds more to larger losses, in contrast to tail-based measures that do not consider intermediate-sized losses. A simple co-measure in this case is \[ r(X_j) = E[X_j \exp(cY / E[Y])] \] Since the exponential moment is scalable, then it has a marginal decomposition. A straightforward (but messy) calculation of the derivative of \(\rho(Y+\epsilon X_j)\) gives \[ \frac{\partial \rho(Y + \epsilon X_j)}{\partial \epsilon}|_{\epsilon = 0} = E[X_j \exp(cY/E[Y])] + \frac{c E[X_j]}{E[Y]}E\left[y\exp(cY/E[Y]) \left(\frac{X_j}{E[X_j]} - \frac{Y}{E[Y]}\right) \right] \]
A generalized moment is an expectation that is not expressable as a power of the variable. For example, if we let \((Y>X)\) denote the function that is equal to 1 when \(Y>X\), and 0 otherwise, then we can express TVaR as a generalized moment: \[ TVaR = E[Y(F_Y(Y) > p)] \] A spectral measure is a generalized moment that can be written in the form \[ \rho = E[Y \eta(F(Y))] \] Therefore, TVaR is an example of a spectral measure. A variation on VaR can be defined by weighting based on the distance from the target percentile: \[ \eta(p) = \exp(-\theta(p - 0.99)^2) \]
option.2.premium = 240
option.3.premium = 110
reinsurance.simulation = simulation.results[, .(simulation_number, line_1_loss, line_2_loss, company_loss)]
reinsurance.simulation[, recoverable_2 := 0.8 * pmax(0, pmin(line_1_loss, 8000) - 4000) + pmax(0, pmin(line_2_loss, 7000) - 3000)]
reinsurance.simulation[, recoverable_3 := pmax(0, pmin(company_loss, 12000) - 7000)]
reinsurance.simulation[, net_results_1 := line.1.premium + line.2.premium - company_loss]
reinsurance.simulation[, net_results_2 := line.1.premium + line.2.premium - company_loss + recoverable_2 - option.2.premium]
reinsurance.simulation[, net_results_3 := line.1.premium + line.2.premium - company_loss + recoverable_3 - option.3.premium]
reinsurance.summary = reinsurance.simulation[, .(net_results_1 = mean(net_results_1),
sd_1 = sd(net_results_1),
min_1 = min(net_results_1),
max_1 = max(net_results_1),
safety_1 = quantile(net_results_1, probs = 0.99),
net_results_2 = mean(net_results_2),
sd_2 = sd(net_results_2),
min_2 = min(net_results_2),
max_2 = max(net_results_2),
safety_2 = quantile(net_results_2, probs = 0.99),
net_results_3 = mean(net_results_3),
sd_3 = sd(net_results_3),
min_3 = min(net_results_3),
max_3 = max(net_results_3),
safety_3 = quantile(net_results_3, probs = 0.99))]
reinsurance.summary
## net_results_1 sd_1 min_1 max_1 safety_1 net_results_2
## 1: 378.7797 1572.862 -9934.364 3737.36 2951.696 322.735
## sd_2 min_2 max_2 safety_2 net_results_3 sd_3 min_3
## 1: 1230.044 -5344.04 3497.36 2711.696 354.3525 1352.56 -5044.364
## max_3 safety_3
## 1: 3627.36 2841.696
Observations about this simulation include:
Visualizing probability densities can also be helpful:
ggplot(data = reinsurance.simulation, aes(x = net_results_1)) + geom_density(color = "red") + geom_density(aes(x = net_results_2), color = "green") + geom_density(aes(x = net_results_3), color = "blue")
Notice that the treaties that provide more protection are shifted to the left, indicating a lessened ability to benefit in good years. Cumulative distribution plots may be more insightful:
ggplot(data = reinsurance.simulation, aes(x = net_results_1)) + stat_ecdf(geom = "step", color = "red") + stat_ecdf(aes(x = net_results_2), geom = "step", color = "green") + stat_ecdf(aes(x = net_results_3), geom = "step", color = "blue")
At a given probability level, the curve furthest to the right performs best. Notice that for most of the graph, Option 3 outperforms Option 2, though there is a range of probabilities in the 5% to 20% range where Option 2 performs better. We can also look at these values in tabluar format:
reinsurance.cdfs = data.table(probability = 1:50 * 0.01)
reinsurance.cdfs$option_1 = apply(reinsurance.cdfs, 1, function(r) quantile(reinsurance.simulation$net_results_1, probs = r['probability']))
reinsurance.cdfs$option_2 = apply(reinsurance.cdfs, 1, function(r) quantile(reinsurance.simulation$net_results_2, probs = r['probability']))
reinsurance.cdfs$option_3 = apply(reinsurance.cdfs, 1, function(r) quantile(reinsurance.simulation$net_results_3, probs = r['probability']))
datatable(reinsurance.cdfs)
Option 2 outperforms Option 3 between the 3-in-100 year level and the 28-in-100 year level. Note that these do not necessarily correspond to the same events, since the CDFs were calculated independently for each reinsurance option. This is essential: we are ranking results according to their own ending probability distributions, not ranking them by the difference with another program. Variations on this approach include:
These ideas can be illustrated using the simulated data, setting a capital requirement of XTVaR at the 90th percentile:
option.1.mean.loss = mean(reinsurance.simulation$company_loss)
option.2.mean.loss = mean(reinsurance.simulation$company_loss - reinsurance.simulation$recoverable_2)
option.3.mean.loss = mean(reinsurance.simulation$company_loss - reinsurance.simulation$recoverable_3)
option.1.VaR = quantile(reinsurance.simulation$company_loss, probs = 0.9)
option.2.VaR = quantile(reinsurance.simulation$company_loss - reinsurance.simulation$recoverable_2, probs = 0.9)
option.3.VaR = quantile(reinsurance.simulation$company_loss - reinsurance.simulation$recoverable_3, probs = 0.9)
option.1.TVaR = reinsurance.simulation[company_loss > option.1.VaR, mean(company_loss)]
option.2.TVaR = reinsurance.simulation[company_loss - recoverable_2 > option.2.VaR, mean(company_loss - recoverable_2)]
option.3.TVaR = reinsurance.simulation[company_loss - recoverable_3 > option.3.VaR, mean(company_loss - recoverable_3)]
option.1.XTVaR = option.1.TVaR - option.1.mean.loss
option.2.XTVaR = option.2.TVaR - option.2.mean.loss
option.3.XTVaR = option.3.TVaR - option.3.mean.loss
option.1.capital.req = option.1.XTVaR
option.2.capital.req = option.2.XTVaR
option.3.capital.req = option.3.XTVaR
option.2.ROE = (mean(reinsurance.simulation$recoverable_2) - option.2.premium) / (option.2.capital.req - option.1.capital.req)
print(paste0("The ROE for option 2 is ", option.2.ROE))
## [1] "The ROE for option 2 is 0.0536912923887152"
option.3.ROE = (mean(reinsurance.simulation$recoverable_3) - option.3.premium) / (option.3.capital.req - option.1.capital.req)
print(paste0("The ROE for option 3 is ", option.3.ROE))
## [1] "The ROE for option 3 is 0.0317172414690246"
Since options 2 and 3 result in a recovery of captial relative to option 1, based on the above, option 3 is preferred. However, the ROE should still be compared to the company’s cost of capital: if it is less than 3.1%, then it would be better to not purchase reinsurance.
The above approach only considered capital requirements over a single year. For long-tailed business, the reserves typically absorb capital for more than one year, and this should be considered when analyzing the value of reinsurance. An additional concern is the accumulation of risks that are correlated across accident years.Copulas provide a method for combining individual distributions into a single multivariate distribution. A key consideration is not only how much correlation there is, but where in the distribution it takes place. A copula is a function \(C(u,v)\) that expresses a joint probability distribution \(F(x,y)\) in terms of the marginal distributions \(F_X(x)\) and \(F_Y(y)\): \[ F(x,y) = C(F_X(x), F_Y(y)) \] Note that this requires that \(C\) be a distribution function on the unit square. Any multivariate distribution can be expressed in this way for some copula. Provide the conditional distribution function \(C(v|u)\) can be inverted, copulas can be simulated by first simulating a random draw of \(u\in [0,1]\), then generating a random \(p \in [0,1]\), then calculating \(C^{-1}(p|u)\). Then, the marginal distributions can be inverted to get \(x = F_X^{-1}(u)\) and \(y = F_Y^{-1}(v)\). Examples of copulas include:
Given a parameter \(a\), this is defined by \[ C(u,v) = -a^{-1} \log(1 + g_ug_v / g_1) \] where \(g_z = e^{-az} - 1\). The Kendall \(\tau\) value for this copula is \[ \tau = 1 - \frac4a + \frac{4}{a^2} \int_0^a \frac{t}{e^t-1}dt \] This can be visualized as follows:
g = function(z, a) {
return(exp(-a * z) - 1)
}
frank.copula = function(u, v, a) {
return(- log(1 + g(u, a) * g(v, a) / g(1, a)) /a)
}
u = 0:100
v = 0:100
copula.illustration = data.table(merge(u, v, by = NULL))
colnames(copula.illustration) = c("u_integer", "v_integer")
copula.illustration[, u := u_integer / 100]
copula.illustration[, v := v_integer / 100]
copula.illustration$frank_copula_CDF = apply(copula.illustration, 1, function(r) frank.copula(u = r['u'], v = r['v'], a = 5))
ggplot(copula.illustration, aes(x = u, y = v, fill = frank_copula_CDF)) + geom_raster() + scale_fill_gradientn(colours = terrain.colors(10)) + geom_contour(aes(z = frank_copula_CDF))
For comparison, look at a the cumulative distribution function for a copula in which \(X\) and \(Y\) are independent:
copula.illustration$independent_CDF = apply(copula.illustration, 1, function(r) r['u'] * r['v'])
ggplot(copula.illustration, aes(x = u, y = v, fill = independent_CDF)) + geom_raster() + scale_fill_gradientn(colours = terrain.colors(10)) + geom_contour(aes(z = independent_CDF))
Properties of the Frank copula include:
The Gumbel copula is defined for \(a\geq 1\): \[ C(u,v) = \exp(-((-\log(u))^a + (- \log(v))^a)^{1/a}) \] The Kendall \(\tau\) value for this copula is \[ \tau = 1 - \frac1a \]
gumbel.copula = function(u, v, a) {
return(exp(-((-log(u))^a + (-log(v))^a)^(1/a)))
}
copula.illustration$gumbel_copula_CDF = apply(copula.illustration, 1, function(r) gumbel.copula(u = r['u'], v = r['v'], a = 1.5))
ggplot(copula.illustration, aes(x = u, y = v, fill = gumbel_copula_CDF)) + geom_raster() + scale_fill_gradientn(colours = terrain.colors(10)) + geom_contour(aes(z = gumbel_copula_CDF))
Properties of the Gumbel copula include:
The Heavy Right Tail copula, for \(a>0\), is \[ C(u,v) = u + v - 1 + ((1 - u)^{-1/a} + (1 - v)^{-1/a} - 1)^{-a} \] The Kendall \(\tau\) value for this copula is \[ \tau = \frac{1}{2a + 1} \]
hrt.copula = function(u, v, a) {
return(u + v - 1 + ((1 - u)^(-1/a) + (1 - v)^(-1/a) -1)^(-a))
}
copula.illustration$hrt_copula_CDF = apply(copula.illustration, 1, function(r) hrt.copula(u = r['u'], v = r['v'], a = 1))
ggplot(copula.illustration, aes(x = u, y = v, fill = hrt_copula_CDF)) + geom_raster() + scale_fill_gradientn(colours = terrain.colors(10)) + geom_contour(aes(z = hrt_copula_CDF))
Properties of the heavy right tail copula include:
Let \(N(x)\) denote the CDF for the standard normal distribution, and let \(B(x, y;a)\) denote the CDF for the bivariate normal distribution. The Normal copula is \[ C(u,v) = B(N^{-1}(u), N^{-1}(v); a) \]
Properties of the normal copula include:The general approach to evaluation of copulas involves calculating functions of the copulas, then looking at these functions for the data. Several copulas are compared, and the one that is the best fit for the data is selected. For example, the right and left tail concentrations are defined as \[ R(z) = \mathrm{Pr}(U> z | V>z) \] and \[ L(z) = \mathrm{Pr}(U < z | V < z) \] The right-tail concentration can be re-expressed as \[ R(z) = \mathrm{Pr}(U>z, V>z) / (1-z) = \frac{1 - 2z + C(z,z)}{1-z} \] and the left tail concentration as \[ L(z) = \mathrm{Pr}(U<z, V<z) / z = \frac{C(z,z)}{z}. \] To illustrate these ideas, compare the Gumbel and Heavy Right Tail copulas. (Note that the parameters selected above ensure that \(\tau = 1/3\) for both copulas, so the comparison will be fair.)
#right.tail.concentration = function(z, copula_data) {
# cp.numerator = copula_data[u_integer == z & v_integer >= z, .(sum(1 - copula_value))][[1]]
# cp.denominator = copula_data[v_integer == z, .(sum(1 - copula_value))][[1]]
# return(cp.numerator / cp.denominator)
#}
gumbel.right = function(z, a) {
(1 - 2 * z + gumbel.copula(z, z, a)) / (1 - z)
}
hrt.right = function(z, a) {
(1 - 2 * z + hrt.copula(z, z, a)) / (1 - z)
}
independent.right = function(z) {
(1 - 2 * z + z^2) / (1 - z)
}
#left.tail.concentration = function(z, copula_data) {
# cp.numerator = copula_data[u_integer == z & v_integer <= z, .(sum(copula_value))][[1]]
# cp.denominator = copula_data[v_integer == z, .(sum(copula_value))][[1]]
# return(cp.numerator / cp.denominator)
#}
gumbel.left = function(z, a) {
gumbel.copula(z, z, a) / z
}
hrt.left = function(z, a) {
hrt.copula(z, z, a) / z
}
independent.left = function(z) {
z
}
tail.concentration.data = data.table(z = 0.001 * 1:999)
tail.concentration.data$gumbel_right = apply(tail.concentration.data, 1, function(r) gumbel.right(r['z'], a = 1.5))
tail.concentration.data$hrt_right = apply(tail.concentration.data, 1, function(r) hrt.right(r['z'], a = 1))
tail.concentration.data$gumbel_left = apply(tail.concentration.data, 1, function(r) gumbel.left(r['z'], a = 1.5))
tail.concentration.data$hrt_left = apply(tail.concentration.data, 1, function(r) hrt.left(r['z'], a = 1))
tail.concentration.data$independent_right = apply(tail.concentration.data, 1, function(r) independent.right(r['z']))
tail.concentration.data$independent_left = apply(tail.concentration.data, 1, function(r) independent.left(r['z']))
ggplot(data = tail.concentration.data, aes(x = z, y = gumbel_right)) + geom_line(colour = "blue") + geom_line(aes(y = hrt_right), colour = "red") + geom_line(aes(y = gumbel_left), colour = "blue") + geom_line(aes(y = hrt_left), colour = "red") + geom_line(aes(y = independent_left), colour = "green") + geom_line(aes(y = independent_right), colour = "green")
Although both the left and right tail concentrations are shown above, they are often graphed so that the left tail concentration is shown below 50%, and the right tail concentration is shown above 50%. These graphs can be compared to the graphs of data to select the copula that best captures the correlation. Note that for many (but not all) copulas, \(\lim_{z\rightarrow 1} R(z) = 0\). As a result, these graphs may be misleading since for some copulas the slope is very steep near 1; there can be significant tail dependence even if the limit is zero. Investigating \(R(z)\) for values close to 1 indicates copula influence on large loss relationships. Other metrics that can be used to select copulas include the overall correlation (e.g. Kendall \(\tau\)), and the \(J\) and \(\chi\) functions.
Let \(T = 1+J\). Since \(E[1+J]=1\), then \(E[A] = E[T]\). Assuming \(J\) is independent of \(A\), then \[ \mathrm{Var}((1+J)A) = \mathrm{Var}(1+J)\mathrm{Var}(A) + \mathrm{Var}(1+J) E[A]^2 + \mathrm{Var}(A) E[(1+J)]^2 \] Therefore, \[ \frac{\mathrm{Var}((1+J)A)}{E[(1+J)A]^2} = \frac{\mathrm{Var}(A)}{E[A]^2}(1 + \mathrm{Var}(1+J)) + \mathrm{Var}(1+J) \] As the coefficient of variation of \(J\) increases, we can assess how the overall coefficient of variation for the company changes as follows:
CV.poisson = function(severity_CV, expected_frequency, cv_J = 0) {
original_cv_squared = 1 / expected_frequency * (severity_CV^2 + 1)
revised_cv_squared = original_cv_squared * (1 + cv_J^2) + cv_J^2
return(sqrt(revised_cv_squared))
}
expected.frequencies = data.frame(expected_frequency = c(2000, 20000, 200000))
CV.J.assumptions = data.frame(CV_assumption = c(0, 0.01, 0.03, 0.05))
CV.example = merge(expected.frequencies, CV.J.assumptions, by = NULL)
CV.example$aggregate_cv = apply(CV.example, 1, function(r) CV.poisson(severity_CV = 7, expected_frequency = r['expected_frequency'], cv_J = r['CV_assumption']))
CV.example = as.data.table(CV.example)
datatable(CV.example)
Note that the relative impact on the larger company is greater – although the introduction of \(J\) increases the volatility of all sizes of companies, the small company was already fairly volatile. Assuming a normally-distributed loss ratio of 65%, we can determine percentiles of the loss ratio as follows:
mean.LR = 0.65
CV.example[, LR_percentile_90 := mean.LR + mean.LR * aggregate_cv * qnorm(0.9)]
CV.example[, LR_percentile_95 := mean.LR + mean.LR * aggregate_cv * qnorm(0.95)]
CV.example[, LR_percentile_99 := mean.LR + mean.LR * aggregate_cv * qnorm(0.99)]
datatable(CV.example)
Without considering parameter risk, the large company is unrealistically stable. Smaller companies are impacted less.
These ideas can be illustrated using simulated data.
set.seed(12345)
random.data = data.table(observation = rgamma(1000, shape = 20, scale = 500))
head(random.data)
## observation
## 1: 11085.668
## 2: 11379.376
## 3: 9510.158
## 4: 9801.123
## 5: 15535.713
## 6: 11660.322
gamma.loglikelihood = function(observations, shape, scale) {
logdensity = sapply(observations, function(x) dgamma(x, shape = shape, scale = scale, log = TRUE))
return(sum(logdensity))
}
shape.parameter.range = data.frame(shape = 5 + 0:10 * 3)
scale.parameter.range = data.frame(scale = 100 + 0:8 * 100)
parameter.ranges = merge(shape.parameter.range, scale.parameter.range, by = NULL)
parameter.ranges$loglikelihood = apply(parameter.ranges, 1, function(r) gamma.loglikelihood(observations = random.data$observation, shape = r['shape'], scale = r['scale']))
ggplot(data = parameter.ranges, aes(x = shape, y = scale, fill = loglikelihood)) + geom_raster() + scale_fill_gradientn(colours = terrain.colors(10)) + geom_contour(aes(z = loglikelihood))
The wide range of high values for the loglikelihood function suggests that there is uncertainty in the parameter estimates. In particular, it is difficult to distinguish from the data whether this is from a high-shape, low-scale Gamma distribution, or a low-shape, high-scale Gamma distribution.
datatable(as.data.table(parameter.ranges)[order(-loglikelihood)])
As an example, adopting a “maintain market share” strategy during a soft market will result in a situation in which premium volume remains the same, but the company takes on more exposure. Once the company can no longer avoid recognition of its growing claims exposure, it may face a rating downgrade or insolvency. A ratings downgrade may force policyholders to switch carriers (failure of stability and availability). Insolvency means policyholders will only see partial recoveries on claims (failure of reliability and affordability, since customers now have to pay a portion of the claim).
A cycle management strategy generally consists of four components:Agency theory attempts to determine ways to align management and owner interests, and to understand the implications of possible divergence. For example, a management incentive plan tied to market capitalization may lead management to be more risk-seeking than the owners, essentially gambling with someone else’s money. On the other hand, stock options may make management more risk-averse than diversified shareholders, resulting in a company run more like a family enterprise. As mentioned above, in insurance this is not just an issue for management incentives, but also underwriter incentives.
Note that at a high level this is essentially the same process as standard actuarial analysis, except that the exposures are typically not covered by any insurance program.
The key distinction between this approach and a traditional business plan is that it is a set of less detailed plans, rather than one highly-detailed (and almost certainly wrong) plan.
In an insurance context, a key application of scenario planning is in determining the planned mix of business, i.e. the planned combination of written premium, price changes, and loss ratios by line of business. The approach assumes that the company can monitor price changes on renewal business, and that loss ratio is a function of price. An example provided in the textbook assumes the following:
base.loss.ratio = 0.8
cost.trend = 1.06
plan.loss.ratio = base.loss.ratio * cost.trend
print(paste0("The plan loss ratio is ", plan.loss.ratio*100, "%"))
## [1] "The plan loss ratio is 84.8%"
A scenario analysis might have produced the following three scenarios:
scenario.table = data.table(Scenario = c("Optimistic", "Realistic", "Pessimistic"), price_change = c(0.05, 0, -0.1), likelihood = c(0.1, 0.5, 0.4))
scenario.table
## Scenario price_change likelihood
## 1: Optimistic 0.05 0.1
## 2: Realistic 0.00 0.5
## 3: Pessimistic -0.10 0.4
Under each scenario, we can calculate the actual loss ratio under each price change to see how it deviates from the plan.
scenario.table[, actual_loss_ratio := plan.loss.ratio / (1 + price_change)]
scenario.table
## Scenario price_change likelihood actual_loss_ratio
## 1: Optimistic 0.05 0.1 0.8076190
## 2: Realistic 0.00 0.5 0.8480000
## 3: Pessimistic -0.10 0.4 0.9422222
The company can now develop response plans for each scenario, e.g. writing more premium in the optimistic scenario, and less premium under the pessimistic scenario. This contrasts with an approach in which the company “sticks to plan” regardless of the price change, which could result in the company writing signifciant premium volume at a high loss ratio. A key benefit of this approach include:
The underwriting cycle is the recurring pattern of increases (hardening) and decreases (softening) in insurance prices and profits. It is illustrated in the following dataset of industry combined ratios by year. (Original data source is A.M. Best; however, since exact values were not provided in the reading I ``eyeballed’’ it for the purpose of having sample data to illustrate concepts.)
industry.COR <- fread("./Brehm_UW_cycle.csv")
ggplot(data = industry.COR, aes(x = Year, y = COR)) + geom_point() + geom_line() + labs(title = "U.S. Property and Casualty Insurance Industry Combined Ratios")
The ratio between premiums and expected losses is viewed as the “price” of insurance for economic purposes, since it accounts for variations in premium that result from coverage changes (limits, retentions, terms and conditions). Key features of the cycle include:
Note that underwriting cycle dynamics vary by stage, which makes it difficult to quantify overall.
Drivers of the underwriting cycle include:This approach involves viewing the profitability measure as a time series, and fitting an autoregressive model to the data. Let \(X_t\) denote the response variable at time \(t\). An \(AR(n)\) model of the form \[ X_t = a + \sum_{1\leq i \leq n} b_x X_{t-i} + \sigma \epsilon_i \] is fit to the data, where \(\epsilon_i\) are independent and identically distributed unit normal random variables. Typically \(n = 2\) or \(n = 3\) is sufficient to model the underwriting cycle. We can fit an \(AR(2)\) model to the data as follows:
ar.2.model <- ar(industry.COR$COR, order.max = 2)
The autoregressive coefficients of the model are
ar.2.model$ar
## [1] 0.8847058 -0.2855619
Note that R fits a model of the form \[ X_t - m = a_1 (X_{t-1} - m) + a_2 (X_{t-2} - m) + \epsilon_t \] where \(m\) is the overall mean. To relate this to the form given in the reading, we’ll collect all the intercept terms on the right side:
TS.intercept = ar.2.model$x.mean * (1 - ar.2.model$ar[1] - ar.2.model$ar[2])
TS.intercept
## [1] 0.4228504
Model accuracy on in-time data can be assessed by calculating one-step-ahead forecasts:
industry.COR[, cor_lag_1 := shift(COR, type = "lag")]
industry.COR[, cor_lag_2 := shift(COR, n = 2, type = "lag")]
industry.COR[, OSE_forecast := TS.intercept + cor_lag_1 * ar.2.model$ar[1] + cor_lag_2 * ar.2.model$ar[2]]
The predictions for the next three points are:
ar.2.predict = predict(ar.2.model, n.ahead = 3)
ar.2.predict$pred
## Time Series:
## Start = 39
## End = 41
## Frequency = 1
## [1] 1.002872 1.030247 1.047934
To visualize these predictions along with the actual data, convert this to a data table, then append it to the original dataset.
prediction.table = data.table("Year" = 2004:2006, "OSE_forecast" = ar.2.predict$pred, "SE" = ar.2.predict$se)
prediction.table[, lower_ci := OSE_forecast - 2 * SE]
prediction.table[, upper_ci := OSE_forecast + 2 * SE]
data.with.predictions = rbind(industry.COR, prediction.table, fill = TRUE)
We can plot the one-step-ahead forecasts along with the predictions to visualize the results. A 95% confidence interval around the predictions is shown in purple.
ggplot(data = data.with.predictions, aes(x = Year, y = COR)) + geom_point(color = "blue") + geom_line(color = "blue") + geom_point(aes(x = Year, y = OSE_forecast), color = "red") + geom_line(aes(x = Year, y = lower_ci), color = "purple") + geom_line(aes(x = Year, y = upper_ci), color = "purple") + labs(title = "U.S. Property and Casualty Insurance Industry Combined Ratios With Time Series Predictions")
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 38 rows containing missing values (geom_path).
## Warning: Removed 38 rows containing missing values (geom_path).
An important use case for these models is to simulate possible future values of the target variable, which can be used as part of a larger ERM or strategic planning model of the firm.
Aggregate approaches based on economic considerations are described below.