Example quantrr risk analysis of the “Widget Management System” (WMS).
library(quantrr)
library(readxl)
library(janitor)
library(validate)
library(dplyr)
library(formattable)
library(purrr)
library(ggplot2)
library(plotly)
# TODO: workaround for https://github.com/r-lib/lintr/issues/2790, update when
# https://github.com/data-cleaning/validate/pull/197 is released (1.1.6+)
is.nzchar <- nzchar # nolint: object_name_linter.
# set the relative file path and name to import
report_file <- "data/widgetsys.xlsx"The Widget Management System (WMS) is over 30 years old and its architecture has not changed significantly since the original implementation. Over the years, the widget system has become an integral part of our services in managing widgets for our clients. In reviewing the system, three major risks were identified: First, the age of the technology prevents updating components of the system that no longer meet contemporary cybersecurity standards, which increases the risk of a breach. Second, the system is less reliable and experiences frequent outages, typically about 2 major outages per year, which results in lost revenue, contractual penalties, and overtime pay to recover from the incident. Third, limitations of the widget system have started to affect sales - we have recently lost a customer due to the functional obsolescence of the widget system, and expect to both lose more existing and prospective customers in the future due to increased competition in the widget management market.
Import and validate data from Excel.
The data was collected from 3 Technology SMEs, 3 Business SMEs, and one SME with experience in both. Experts were calibrated, informed by historical and industry data, and only gave estimates for areas in which they were confident in answering.
risks <- read_xlsx(report_file, sheet = "Risks") |>
clean_names()
validate_risks <- local({
validate_rules <- validator(
risk_char = is.character(risk),
risk_not_na = !is.na(risk),
risk_not_blank = is.nzchar(risk, keep_na = TRUE),
desc_char = is.character(description),
desc_not_na = !is.na(description),
desc_not_blank = is.nzchar(risk, keep_na = TRUE)
)
confront(risks, validate_rules)
})
check_validation(validate_risks, sheet = "Risks")
estimates <- read_xlsx(report_file, sheet = "Estimates") |>
clean_names() |>
rename(
lambda = frequency_per_yer, p05 = low_5_percent, p95 = high_95_percent, p50 = most_likely
)
validate_estimates <- local({
validate_rules <- validator(
risk_not_na = !is.na(risk),
risk_match = risk %in% risks$risk,
expert_char = is.character(expert),
expert_not_na = !is.na(expert),
expert_not_blank = is.nzchar(expert, keep_na = TRUE),
lambda_num = is.numeric(lambda),
lambda_pos = lambda > 0,
p05_num = is.numeric(p05),
p05_pos = p05 > 0,
p95_num = is.numeric(p95),
p95_pos = p95 > 0,
p50_num = is.numeric(p50),
p50_pos = p50 > 0
)
confront(estimates, validate_rules)
})
check_validation(validate_estimates, sheet = "Estimates")Risk descriptions:
formattable(risks, align = "l")| risk | description |
|---|---|
| Cybersecurity Breach | Risk of a cybersecurity breach of the inventory system. |
| Technology Outage | Risk of an inventory system outage. |
| Loss of Customer | Risk of losing and existing customer or failure to acquire a new customer due to functional limitations of the inventory system. |
Forecast risk using Monte Carlo simulation. The average events, losses, ‘typical’ losses (geometric mean), and percentage of years with no losses for each risk are summarized below:
consensus <- estimates |>
group_by(risk) |>
summarize(across(lambda:p50, ~ mean(.x, na.rm = TRUE)))
consensus_params <- consensus |>
mutate(as_tibble(lnorm_param(.data$p05, .data$p95, .data$p50)))
forecast <- consensus_params |>
select(c("risk", "lambda", "meanlog", "sdlog")) |>
pmap(calc_risk) |>
list_rbind()
forecast |>
group_by(risk) |>
mutate(no_losses = events == 0) |>
summarize(
avg_events = mean(events), avg_losses = mean(losses), typ_losses = gmean(losses),
no_losses = mean(no_losses)
) |>
mutate(
across(avg_losses:typ_losses, ~ currency(.x, digits = 0L)),
no_losses = percent(no_losses)
) |>
formattable(align = "l")| risk | avg_events | avg_losses | typ_losses | no_losses |
|---|---|---|---|---|
| Cybersecurity Breach | 0.25798 | $6,187,796 | $362,774 | 77.28% |
| Loss of Customer | 1.99890 | $10,374,350 | $4,284,572 | 13.14% |
| Technology Outage | 1.98968 | $96,197 | $35,069 | 13.63% |
Losses by risk separately and in aggregate:
forecast |>
# remove zero losses to plot using log10 scale
filter(losses > 0) |>
ggplot(aes(losses)) +
facet_grid(vars(risk)) +
geom_histogram(color = "black", fill = "white", bins = 100) +
scale_x_log10(labels = scales::label_currency(scale_cut = scales::cut_short_scale())) +
labs(x = NULL, y = NULL, title = "Losses by Risk") +
theme_minimal()
forecast |>
filter(losses > 0) |>
ggplot(aes(losses, fill = risk)) +
geom_histogram(bins = 100) +
scale_x_log10(labels = scales::label_currency(scale_cut = scales::cut_short_scale())) +
labs(x = NULL, y = NULL, fill = "Risk", title = "Total Losses") +
scale_fill_viridis_d() +
theme_minimal()
Plot loss exceedance curves for all risks and combined risk.
Plot loss exceedance curves for each risk:
risk_le <- forecast |>
group_by(risk) |>
mutate(probability = 1 - percent_rank(losses)) |>
filter(losses > 0) |>
ggplot(aes(losses, probability)) +
facet_grid(vars(risk)) +
geom_line() +
scale_y_continuous(labels = scales::label_percent(), limits = c(0, 1)) +
scale_x_log10(labels = scales::label_currency(scale_cut = scales::cut_short_scale())) +
labs(x = NULL, y = NULL, title = "Loss Exceedance by Risk") +
theme_minimal()
risk_le
Interactive plot:
ggplotly(risk_le)Plot loss exceedance curve for combined risk:
combined_le <- forecast |>
group_by(year) |>
summarize(total_losses = sum(losses)) |>
mutate(probability = 1 - percent_rank(total_losses)) |>
filter(total_losses > 0) |>
ggplot(aes(total_losses, probability)) +
geom_line() +
scale_y_continuous(labels = scales::label_percent(), limits = c(0, 1)) +
scale_x_log10(labels = scales::label_currency(scale_cut = scales::cut_short_scale())) +
labs(x = NULL, y = NULL, title = "Combined Loss Exceedance") +
theme_minimal()
combined_le
Interactive plot:
ggplotly(combined_le)Additional details on the risk quantification analysis.
Data validation results for Risks tab:
plot(validate_risks)
Data validation results for Estimates tab:
plot(validate_estimates)
All risk estimates:
estimates |>
mutate(across(p05:p50, ~ currency(.x, digits = 0L))) |>
formattable(align = "l")| risk | expert | lambda | p05 | p95 | p50 |
|---|---|---|---|---|---|
| Cybersecurity Breach | Technology Expert 1 | 0.25 | NA | NA | NA |
| Cybersecurity Breach | Technology Expert 2 | 0.33 | NA | NA | NA |
| Cybersecurity Breach | Technology Expert 3 | 0.20 | NA | NA | NA |
| Cybersecurity Breach | Business Expert 1 | NA | $2,000 | $45,000,000 | $300,000 |
| Cybersecurity Breach | Business Expert 2 | NA | $1,400 | $34,000,000 | $200,000 |
| Cybersecurity Breach | Business Expert 3 | NA | $2,100 | $54,000,000 | $500,000 |
| Cybersecurity Breach | Unicorn Expert | 0.25 | $1,900 | $44,000,000 | $275,000 |
| Technology Outage | Technology Expert 1 | 2.00 | NA | NA | NA |
| Technology Outage | Technology Expert 2 | 3.00 | NA | NA | NA |
| Technology Outage | Technology Expert 3 | 1.00 | NA | NA | NA |
| Technology Outage | Business Expert 1 | NA | $700 | $200,000 | $15,000 |
| Technology Outage | Business Expert 2 | NA | $500 | $150,000 | $10,000 |
| Technology Outage | Business Expert 3 | NA | $675 | $180,000 | $11,000 |
| Technology Outage | Unicorn Expert | 2.00 | $800 | $220,000 | $16,000 |
| Loss of Customer | Technology Expert 1 | NA | NA | NA | NA |
| Loss of Customer | Technology Expert 2 | NA | NA | NA | NA |
| Loss of Customer | Technology Expert 3 | NA | NA | NA | NA |
| Loss of Customer | Business Expert 1 | 2.00 | $100,000 | $20,000,000 | $1,500,000 |
| Loss of Customer | Business Expert 2 | 3.00 | $150,000 | $30,000,000 | $2,000,000 |
| Loss of Customer | Business Expert 3 | 1.00 | $75,000 | $15,000,000 | $1,200,000 |
| Loss of Customer | Unicorn Expert | 2.00 | $75,000 | $15,000,000 | $1,000,000 |
Using a simple average of all experts that provided an estimate (not blank/NA), this gives us a consensus estimate for the three risks of:
consensus |>
mutate(across(p05:p50, ~ currency(.x, digits = 0L))) |>
formattable(align = "l")| risk | lambda | p05 | p95 | p50 |
|---|---|---|---|---|
| Cybersecurity Breach | 0.2575 | $1,850 | $44,250,000 | $318,750 |
| Loss of Customer | 2.0000 | $100,000 | $20,000,000 | $1,425,000 |
| Technology Outage | 2.0000 | $669 | $187,500 | $13,000 |
The consensus estimates for p05 and p95 result in the following parameters for log-normal loss magnitude. The p50 estimate is used to calculate the percentage difference from the actual median (mdiff), a measure of estimate accuracy:
consensus_params |>
mutate(across(p05:p50, ~ currency(.x, digits = 0L)), mdiff = percent(mdiff)) |>
formattable(align = "l")| risk | lambda | p05 | p95 | p50 | meanlog | sdlog | mdiff |
|---|---|---|---|---|---|---|---|
| Cybersecurity Breach | 0.2575 | $1,850 | $44,250,000 | $318,750 | 12.564153 | 3.064840 | 11.41% |
| Loss of Customer | 2.0000 | $100,000 | $20,000,000 | $1,425,000 | 14.162084 | 1.610574 | 0.76% |
| Technology Outage | 2.0000 | $669 | $187,500 | $13,000 | 9.323472 | 1.713260 | 16.09% |
A summary() of the forecast results.
summary(forecast) year risk events losses
Min. : 1 Length:150000 Min. : 0.000 Min. :0.000e+00
1st Qu.:12501 Class :character 1st Qu.: 0.000 1st Qu.:0.000e+00
Median :25000 Mode :character Median : 1.000 Median :3.561e+04
Mean :25000 Mean : 1.416 Mean :5.553e+06
3rd Qu.:37500 3rd Qu.: 2.000 3rd Qu.:1.464e+06
Max. :50000 Max. :10.000 Max. :1.358e+10