This markdown script reproduces the figures and numerical results associated with the paper “Temporal Compounding Increases Economic Impacts of Atmospheric Rivers in California” (doi:XX), published in Science Advances. Please note that the figures in the markdown file may not exactly match the published version because of formatting and scale differences.
source('_data/setup_impacts.R')
source('_scripts/create_df_functions.R')
## set number of bootstrapped samples
<- 1000
boot
## turn progress bars on/off
<- FALSE
progress if (progress) pb <- txtProgressBar(min = 0, max = ncell(grid_ca), style = 3)
## set random seed for reproducibility
set.seed(2023)
load('_scripts/_checkpoints/df_3hr_1209.Rdata')
load('_scripts/_checkpoints/df_24hr_0209.Rdata')
<-
df_3hr foreach (i = 1:ncell(grid_ca)) %do% {
if (progress) setTxtProgressBar(pb, i)
if (i %in% index_ca) {
<- df_3hr[[i]] %>%
temp mutate(
ar = case_when(ar.cat==0 ~ FALSE, TRUE ~ ar),
ar.cat = case_when(ar ~ ar.cat),
ar.count = case_when(ar ~ ar.count)) %>%
mutate(inter = !ar, inter.count = add_counter(inter))
<- create_catalog(temp, 'inter', cat = FALSE, interval = 3/24) %>%
temp select(count, duration) %>%
setNames(paste('inter', names(.), sep = '.')) %>%
left_join(temp, ., by = 'inter.count')
<- temp %>%
temp select(ar.count, ar.cat, inter.duration) %>%
mutate(
inter.duration = setNA(inter.duration,0),
prev.inter = c(NA, inter.duration[-nrow(.)]),
next.inter = c(inter.duration[-1], NA)) %>%
filter(!is.na(ar.count)) %>%
group_by(ar.count = toNumber(ar.count)) %>%
summarize(
prev.inter = prev.inter[1],
next.inter = next.inter[length(ar.count)],
ar.cat = ar.cat[1]) %>%
mutate(
prev.cat = c(NA, ar.cat[-nrow(.)]),
next.cat = c(ar.cat[-1], NA),
prev.inter = c(NA, prev.inter[-1]),
next.inter = c(next.inter[-nrow(.)], NA)) %>%
select(-ar.cat)
%>% left_join(temp, by = 'ar.count')
df_3hr[[i]] else NULL
} }
<-
df_catalog foreach (i = 1:ncell(grid_ca)) %do% {
if (progress) setTxtProgressBar(pb, i)
if (i %in% index_ca) {
left_join(
%>%
df_3hr[[i]] ###
mutate(
ar = case_when(ar.cat==0 ~ FALSE, TRUE ~ ar),
ar.cat = case_when(ar ~ ar.cat),
ar.count = case_when(ar ~ ar.count)) %>%
###
filter(ar) %>% select(-ar) %>%
group_by(ar.count) %>%
summarize(
across(starts_with(c('ar', 'seq', 'prev', 'next')), ~.[1]),
precip = sum(precip),
sm = mean(sm)) %>%
select(-ends_with('end')),
%>%
df_24hr[[i]] filter(ar) %>%
group_by(ar.count = toNumber(ar.count)) %>%
summarize(
ar.maxrunoff = max(runoff),
ar.totalrunoff = sum(runoff),
ar.maxnorm_rp2 = max(norm_rp2),
across(c(disaster, pa_amt, ncei_event, ncei_damage, starts_with(c('claims', 'policies'))), Sum),
pop = mean(pop),
wwa_advisory = Sum(wwa_vtec=='advisory'),
wwa_watch = Sum(wwa_vtec=='watch'),
wwa_warning = Sum(wwa_vtec=='warning')),
by = 'ar.count') %>%
mutate(id = i)
else NULL
} }
## keep WY 1981-2021 (for category analysis)
<- df_catalog %>%
catalog.full reduce(rbind) %>%
mutate(wy = wateryear(ar.start)) %>%
filter(wateryear(ar.start) %in% 1981:2021)
## keep WY 1997-2021 (for impact analysis)
<- catalog.full %>%
catalog filter(wateryear(ar.start) %in% 1997:2021) %>%
mutate(
ncei_damage = setNA(ncei_damage,0),
claims_value = setNA(claims_value,0))
cat('Number of AR events in the catalog: ')
cat(comma(nrow(catalog)))
## Number of AR events in the catalog: 23,267
cat('Percent of events with nonzero NCEI loss: ')
cat(
%>% count(loss = ncei_damage > 0) %>% mutate(pct = n/sum(n)) %>%
catalog filter(loss) %>% pull(pct) %>% percent(accuracy = 0.1))
cat('\n')
cat('Percent of events with nonzero NFIP loss: ')
cat(
%>% count(loss = claims_value > 0) %>% mutate(pct = n/sum(n)) %>%
catalog filter(loss) %>% pull(pct) %>% percent(accuracy = 0.1))
cat('\n')
## Percent of events with nonzero NCEI loss: 7.9%
## Percent of events with nonzero NFIP loss: 2.1%
Probability of another AR within ±5 days
<-
prox foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb, i)
if (i %in% index_ca) {
%>%
df_3hr[[i]] filter(wateryear(ts) %in% 1981:2021) %>%
filter(ar.cat > 0) %>%
group_by(ar.count) %>%
summarize(ar.cat = ar.cat[1], before = prev.inter[1]<=5, after = next.inter[1]<=5) %>%
group_by(ar.cat) %>%
summarize(
n = length(ar.cat),
and = Sum(before & after)/n, or = Sum(before | after)/n) %>%
mutate(id = i)
}%>%
} left_join(raster.df(grid_ca), by = c('id' = 'value'))
<- map(
g.or .x = 1:5, .f = ~prox %>%
filter(ar.cat == .x) %>%
ggplot() +
geom_point(aes(x=x, y=y, color = or, size = n)) +
geom_sf(data = cal, fill = NA, color = 'black', size = 0.25) +
scale_color_distiller(
'Percentage of\nAR Events', labels = percent,
palette = 'Purples', limits = c(0,1), direction = 1, na.value = NA) +
scale_size_area(max_size = 2.5) +
guides(size = guide_none()) + #, color = guide_colorbar(direction = 'horizontal')) +
theme(
text = element_text(family = 'Segoe UI', size = 9),
strip.background = element_rect(color = NA, fill = 'grey95'),
strip.placement = 'outside',
strip.text = element_text(family = 'Segoe UI', size = 9),
axis.line = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
plot.margin = margin(0, 0, 0, 0, "pt"),
axis.title = element_blank(), axis.text.x = element_blank(),
legend.title = element_text(margin = margin(0,0,5,0))))
1]] <- g.or[[1]] +
g.or[[facet_grid('Adjacent ARs' ~ 'AR1', switch = 'y')
2]] <- g.or[[2]] +
g.or[[facet_grid(.~paste0('AR', ar.cat)) +
theme(axis.text.y = element_blank())
3]] <- g.or[[3]] +
g.or[[facet_grid(.~paste0('AR', ar.cat)) +
theme(axis.text.y = element_blank())
4]] <- g.or[[4]] +
g.or[[facet_grid(.~paste0('AR', ar.cat)) +
theme(axis.text.y = element_blank())
5]] <- g.or[[5]] +
g.or[[facet_grid(.~paste0('AR', ar.cat)) +
theme(axis.text.y = element_blank())
<- map(
g.and .x = 1:5, .f = ~prox %>%
filter(ar.cat == .x) %>%
ggplot() +
geom_point(aes(x=x, y=y, color = and, size = n)) +
geom_sf(data = cal, fill = NA, color = 'black', size = 0.25) +
scale_color_distiller(
'Percentage of\nAR Events', labels = percent,
palette = 'Blues', limits = c(0,0.5), direction = 1, na.value = 'darkblue') +
scale_size_area(max_size = 2.5) +
guides(size = guide_none()) + #, color = guide_colorbar(direction = 'horizontal')) +
theme(
text = element_text(family = 'Segoe UI', size = 9),
strip.background = element_rect(color = NA, fill = 'grey95'),
strip.placement = 'outside',
strip.text = element_text(family = 'Segoe UI', size = 9),
axis.line = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
plot.margin = margin(0, 0, 0, 0, "pt"),
axis.title = element_blank(),
legend.title = element_text(margin = margin(0,0,5,0))))
1]] <- g.and[[1]] + facet_grid('Sandwiched ARs'~., switch = 'y')
g.and[[2]] <- g.and[[2]] + theme(axis.text.y = element_blank())
g.and[[3]] <- g.and[[3]] + theme(axis.text.y = element_blank())
g.and[[4]] <- g.and[[4]] + theme(axis.text.y = element_blank())
g.and[[5]] <- g.and[[5]] + theme(axis.text.y = element_blank()) g.and[[
<- wrap_plots(g.or) + guide_area() + plot_layout(design = 'abcdef', guides = 'collect')
plot.or <- wrap_plots(g.and) + guide_area() + plot_layout(design = 'abcdef', guides = 'collect')
plot.and
/ plot.and
plot.or ggsave('_figures-impacts/fig2_coincidence.png', width=7.25, height=3.25, units = 'in', dpi = 500)
%>% count(ar.cat) catalog.full
## # A tibble: 5 x 2
## ar.cat n
## <dbl> <int>
## 1 1 21078
## 2 2 10358
## 3 3 4812
## 4 4 1230
## 5 5 323
cat('How much more likely are Category 5 events to be adjacent, compared to Category 1?\n')
<- catalog.full %>%
adjacent count(ar.cat, five = prev.inter <= 5 | next.inter <= 5) %>%
group_by(ar.cat) %>%
mutate(pct = n/sum(n)) %>%
filter(five)
$pct[5] / adjacent$pct[1]) %>% comma(accuracy = 0.1, suffix = 'x')
(adjacentcat('\n')
cat('How much more likely are Category 5 events to be sandwiched, compared to Category 1?\n')
<- catalog.full %>%
sandwich count(ar.cat, five = prev.inter <= 5 & next.inter <= 5) %>%
group_by(ar.cat) %>%
mutate(pct = n/sum(n)) %>%
filter(five)
$pct[5] / sandwich$pct[1]) %>% comma(accuracy = 0.1, suffix = 'x') (sandwich
## How much more likely are Category 5 events to be adjacent, compared to Category 1?
## [1] "1.7x"
##
## How much more likely are Category 5 events to be sandwiched, compared to Category 1?
## [1] "2.6x"
Relationship between primary event category and preceding event category
<- catalog.full %>% filter(prev.inter <= 5)
catalog5
## Report number of events
cat('Number of ARs within five days of another AR: ')
cat(comma(nrow(catalog5)))
<-
ratios map_dfr(
.x = 1:boot,
.f = ~sample(
$ar.cat, size = nrow(catalog5), replace = TRUE) %>%
catalog5%>% prop.table %>% cumsum) %>%
table apply(2, function(x) quantile(x, c(0.05,0.5,0.95))) %>%
%>% as.data.frame %>%
t setNames(c('q05', 'med', 'q95')) %>%
rownames_to_column('ar.cat') %>%
mutate(ar.cat = factor(ar.cat, levels = 5:1)) %>%
mutate(med.diff = diff(c(0,med)))
## Number of ARs within five days of another AR: 10,616
<- ggplot(ratios) +
g.left geom_col(
aes(x = 'All AR Events', y = med.diff, group = ar.cat),
fill = 'white', color = NA, width = 0.65) +
geom_col(
aes(x = 'All AR Events', y = med.diff, group = ar.cat, color = ar.cat, fill = ar.cat),
width = 0.65, alpha = 0.65, size = 0.35) +
geom_errorbar(
data = ratios[-5,], aes(x = 'All AR Events', ymin = q05, ymax = q95),
size = 0.35, width = 0.15) +
geom_hline(
data = ratios[-5,], aes(yintercept = med),
size = 0.25, linetype = 'dotted', show.legend = FALSE) +
geom_ribbon(
data = ratios, aes(x = seq(0,7,length.out=5), ymin = q05[1], ymax = q95[1]),
fill = 'grey25', alpha = 0.15) +
geom_ribbon(
data = ratios, aes(x = seq(0,7,length.out=5), ymin = q05[2], ymax = q95[2]),
fill = 'grey25', alpha = 0.15) +
geom_ribbon(
data = ratios, aes(x = seq(0,7,length.out=5), ymin = q05[3], ymax = q95[3]),
fill = 'grey25', alpha = 0.15) +
geom_ribbon(
data = ratios, aes(x = seq(0,7,length.out=5), ymin = q05[4], ymax = q95[4]),
fill = 'grey25', alpha = 0.15) +
scale_color_manual(values = roma.colors) +
scale_fill_manual(values = roma.colors) +
guides(color = guide_none(), fill = guide_none()) +
theme(
text = element_text(family = 'Segoe UI', size = 9),
axis.title = element_blank(),
axis.ticks.y = element_line(size = 0.25),
plot.margin = margin(t = 2, b = 10, r = -50),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major.y = element_line(size = 0.25))
<- g.left +
g1 scale_y_continuous(labels = percent, expand = expansion(mult = 0)) +
coord_cartesian(xlim = c(0.95,1.35), ylim = c(0,0.85))
<- g.left +
g2 scale_y_continuous(labels = percent) +
coord_cartesian(xlim = c(0.95,1.35), ylim = c(0.85,1)) +
theme(
axis.title.x = element_blank(), axis.text.x = element_blank(),
axis.line.x = element_blank(), axis.ticks.x = element_blank())
<- catalog5 %>%
g.right count(ar.cat, prev.cat = factor(prev.cat, levels = 5:1)) %>%
ggplot() +
geom_col(
aes(x = ar.cat, y = n, group = prev.cat),
position = 'fill', fill = 'white', color = NA, width = 0.65) +
geom_col(
aes(x = ar.cat, y = n, group = prev.cat, color = prev.cat, fill = prev.cat),
position = 'fill', width = 0.65, alpha = 0.65, size = 0.35) +
geom_hline(
data = ratios[-5,], aes(yintercept = med),
size = 0.25, linetype = 'dotted', show.legend = FALSE) +
geom_ribbon(
data = ratios, aes(x = seq(0,7,length.out=5), ymin = q05[1], ymax = q95[1]),
fill = 'grey25', alpha = 0.15) +
geom_ribbon(
data = ratios, aes(x = seq(0,7,length.out=5), ymin = q05[2], ymax = q95[2]),
fill = 'grey25', alpha = 0.15) +
geom_ribbon(
data = ratios, aes(x = seq(0,7,length.out=5), ymin = q05[3], ymax = q95[3]),
fill = 'grey25', alpha = 0.15) +
geom_ribbon(
data = ratios, aes(x = seq(0,7,length.out=5), ymin = q05[4], ymax = q95[4]),
fill = 'grey25', alpha = 0.15) +
scale_color_manual('Previous\nEvent Rank', labels = paste0('AR', 5:1), values = roma.colors) +
scale_fill_manual('Previous\nEvent Rank', labels = paste0('AR', 5:1), values = roma.colors) +
scale_x_continuous(breaks = 1:5, labels = paste0('AR', 1:5)) +
guides(color = guide_legend(reverse = TRUE), fill = guide_legend(reverse = TRUE)) +
theme(
text = element_text(family = 'Segoe UI', size = 9),
legend.position = 'bottom',
legend.margin = margin(10,5,5,5),
legend.background = element_rect(fill = NA),
axis.title.y = element_blank(), axis.text.y = element_blank(),
axis.ticks.y = element_blank(), axis.line.y = element_line(color = NA),
# axis.ticks.y = element_line(size = 0.25),
plot.margin = margin(t = 2, b = 10, l = -50, r = 0),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major.y = element_line(size = 0.25))
<- g.right +
g3 scale_y_continuous(labels = percent, expand = expansion(mult = 0)) +
labs(x = 'Primary Event Rank') +
coord_cartesian(xlim = c(0.65,5.15), ylim = c(0,0.85))
<- g.right +
g4 coord_cartesian(xlim = c(0.65,5.15), ylim = c(0.85,1)) +
scale_y_continuous(labels = percent) +
theme(
axis.title.x = element_blank(), axis.text.x = element_blank(),
axis.line.x = element_blank(), axis.ticks.x = element_blank())
<-
gg_axis ::get_plot_component(ggplot() + labs(y = ' Percent of Occurrences'), 'ylab-l')
cowplot
+ g1 + g4 + g3 + plot_spacer() + guide_area() + plot_spacer() + gg_axis + plot_spacer() +
g2 plot_layout(design = 'hiaec\nhibed\ngifff', widths = c(1,-10,26,-8,80), heights = c(5,8,1), guides = 'collect')
ggsave('_figures-impacts/fig3_prevcat.png', width = 3.55, height = 4.25, units = 'in', dpi = 500)
<- catalog5 %>%
temp count(ar.cat, prev.cat = factor(prev.cat, levels = 5:1)) %>%
group_by(ar.cat) %>% mutate(obs = n/sum(n)) %>%
left_join(ratios %>% transmute(prev.cat = factor(ar.cat), exp = med.diff), by = 'prev.cat') %>%
mutate(diff = obs/exp) %>%
filter(ar.cat == 5)
cat('How many more Category 5 events are preceded by Category 5 events than expected?\n')
%>%
temp filter(prev.cat == 5) %>%
pull(diff) %>% comma(accuracy = 0.1, suffix = 'x')
cat('\n')
cat('How many fewer Category 5 events are preceded by Category 1 events than expected?\n')
%>%
temp filter(prev.cat == 1) %>%
pull(diff) %>% comma(accuracy = 0.01, suffix = 'x')
## How many more Category 5 events are preceded by Category 5 events than expected?
## [1] "2.2x"
##
## How many fewer Category 5 events are preceded by Category 1 events than expected?
## [1] "0.31x"
<-
daily map_dfr(
.x = index_ca,
.f = ~df_24hr[[.x]] %>% select(-ends_with('rp2')) %>% mutate(id = .x)) %>%
mutate(
impactday = claims_num>0 | ncei_event,
lossday = claims_value>0 | ncei_damage>0) %>%
mutate(wy = wateryear(date)) %>%
filter(wy %in% 1997:2021)
cat('Number of daily records: ')
cat(comma(nrow(daily)))
## Number of daily records: 1,506,615
cat('Statewide average percent of days that are labeled impact days, 1996-2021:\n')
%>%
daily count(impactday) %>%
mutate(pct = n/sum(n)) %>%
filter(impactday) %>%
pull(pct) %>%
percent(accuracy = 0.1)
cat('\n')
cat('Minimum and maximum individual individual grid cells:\n')
%>%
daily count(id, impactday) %>%
group_by(id) %>%
mutate(pct = n/sum(n)) %>%
filter(impactday) %>%
pull(pct) %>% range %>%
percent(accuracy = 0.1)
## Statewide average percent of days that are labeled impact days, 1996-2021:
## [1] "4.8%"
##
## Minimum and maximum individual individual grid cells:
## [1] "0.7%" "13.1%"
<- felm(
f2 ~ ar*seq - seq | factor(id) + factor(wy),
impactday data = daily, exactDOF = TRUE)
summary(f2, robust = TRUE)
##
## Call:
## felm(formula = impactday ~ ar * seq - seq | factor(id) + factor(wy), data = daily, exactDOF = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35020 -0.06714 -0.03596 -0.00723 1.04170
##
## Coefficients:
## Estimate Robust s.e t value Pr(>|t|)
## arTRUE 0.109770 0.001522 72.10 <2e-16 ***
## arFALSE:seqTRUE 0.054049 0.000933 57.93 <2e-16 ***
## arTRUE:seqTRUE 0.082585 0.002544 32.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2082 on 1506423 degrees of freedom
## Multiple R-squared(full model): 0.05648 Adjusted R-squared: 0.05636
## Multiple R-squared(proj model): 0.02809 Adjusted R-squared: 0.02797
## F-statistic(full model, *iid*):472.1 on 191 and 1506423 DF, p-value: < 2.2e-16
## F-statistic(proj model): 5323 on 3 and 1506423 DF, p-value: < 2.2e-16
Impact day likelihood regression coefficients
summary(f2, robust = TRUE) %>%
%>% data.frame %>%
coef rownames_to_column() %>%
transmute(
est = Estimate,
se = Robust.s.e,
ar = c(TRUE,FALSE,TRUE),
seq = c(NA,TRUE,TRUE)) %>%
group_by(ar) %>%
mutate(est.cumsum = cumsum(est)) %>%
%>%
ungroup mutate(lower = est.cumsum - 2*se, upper = est.cumsum + 2*se) %>%
ggplot() +
geom_col(aes(x = ar, y = est, group = seq, fill = seq), width = 0.6) +
geom_errorbar(aes(x = ar, ymin = lower, ymax = upper), width = 0.1, size = 0.5) +
scale_y_origin(
'Increase in Probability of Impact', labels = percent) +
scale_x_discrete('', labels = c('Non-AR Day', 'AR Day')) +
scale_fill_manual(
'In a Sequence?', labels = c('Yes', ''),
values = 'gray30', na.value = 'grey90') +
guides(fill = guide_legend(override.aes = list(fill = c('gray30','white')))) +
theme(
text = element_text(family = 'Segoe UI', size = 9),
legend.position = c(0.75, 0.25),
legend.background = element_rect(fill = NA),
panel.grid.major.x = element_line(size = 0.25)) +
coord_flip()
ggsave('_figures-impacts/fig4_likelihood.png', width = 3.55, height = 2.5, dpi = 500)
# colorRampPalette(c('navy','white'))(9)
## NCEI
<- which(catalog$ncei_damage == 0)
index.ncei.zero <- which(catalog$ncei_damage > 0)
index.ncei.nonzero <-
coef.under foreach (b = 1:boot) %do% {
felm(
log10(ncei_damage+1) ~ factor(ar.cat) + seq | factor(id) + factor(wy),
data = catalog %>%
slice(c(index.ncei.nonzero, sample(
size = 2*length(index.ncei.nonzero), replace = FALSE))),
index.ncei.zero, exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),mean)
<- coef.under %>% data.frame %>%
df.ncei transmute(
cat = c(2:5,NA), seq = c(rep(NA,4), TRUE),
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
rbind(data.frame(cat = 1, seq = NA, est = 1, lower = 1, upper = 1))
## NFIP
<- which(catalog$claims_value == 0)
index.nfip.zero <- which(catalog$claims_value > 0)
index.nfip.nonzero <-
coef.under foreach (b = 1:boot) %do% {
felm(
log10(claims_value+1) ~ factor(ar.cat) + seq | factor(id) + factor(wy),
data = catalog %>%
slice(c(index.nfip.nonzero, sample(
size = 2*length(index.nfip.nonzero), replace = FALSE))),
index.nfip.zero, exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),mean)
<- coef.under %>% data.frame %>%
df.nfip transmute(
cat = c(2:5,NA), seq = c(rep(NA,4), TRUE),
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
rbind(data.frame(cat = 1, seq = NA, est = 1, lower = 1, upper = 1))
%>% filter(is.na(seq)) %>%
df.nfip select(cat, lower, est, upper) %>%
arrange(cat) %>%
%>%
gt tab_header('NFIP Category Loss Multipliers') %>%
fmt_markdown(cat) %>%
fmt_number(c(lower, est, upper), n_sigfig = 3) %>%
cols_label(
cat = 'AR Category',
lower = 'Lower Bound',
est = 'Median Estimate',
upper = 'Upper Bound') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
rbind(
%>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NFIP')) %>%
df.nfip select(cat, est, data) %>%
pivot_wider(names_from = data, values_from = est) %>%
arrange(cat) %>%
mutate(diff = (NCEI / NFIP) - 1) %>%
%>%
gt tab_header('Category Loss Multipliers') %>%
fmt_markdown(cat) %>%
fmt_number(c(NCEI, NFIP), , n_sigfig = 3) %>%
fmt_percent(diff, decimals = 1) %>%
cols_label(
cat = 'AR Category',
diff = 'Difference') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
NFIP Category Loss Multipliers | |||
---|---|---|---|
AR Category | Lower Bound | Median Estimate | Upper Bound |
1 |
1.00 | 1.00 | 1.00 |
2 |
2.00 | 3.46 | 5.98 |
3 |
10.8 | 20.8 | 40.1 |
4 |
27.7 | 78.9 | 225 |
5 |
347 | 1,570 | 7,140 |
Category Loss Multipliers | |||
---|---|---|---|
AR Category | NCEI | NFIP | Difference |
1 |
1.00 | 1.00 | 0.0% |
2 |
2.28 | 3.46 | −33.9% |
3 |
14.7 | 20.8 | −29.2% |
4 |
45.4 | 78.9 | −42.5% |
5 |
177 | 1,570 | −88.8% |
cat('Sequence loss multipliers\n')
cat(' NCEI: ')
cat(df.ncei %>% filter(!is.na(seq)) %>% pull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat(' (90% CI: ')
cat(df.ncei %>% filter(!is.na(seq)) %>% pull(lower) %>% comma(accuracy = 0.1, suffix = 'x'))
cat(', ')
cat(df.ncei %>% filter(!is.na(seq)) %>% pull(upper) %>% comma(accuracy = 0.1, suffix = 'x'))
cat(')\n')
cat(' NFIP: ')
cat(df.nfip %>% filter(!is.na(seq)) %>% pull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat(' (90% CI: ')
cat(df.nfip %>% filter(!is.na(seq)) %>% pull(lower) %>% comma(accuracy = 0.1, suffix = 'x'))
cat(', ')
cat(df.nfip %>% filter(!is.na(seq)) %>% pull(upper) %>% comma(accuracy = 0.1, suffix = 'x'))
cat(')\n')
## Sequence loss multipliers
## NCEI: 4.2x (90% CI: 3.2x, 5.5x)
## NFIP: 3.4x (90% CI: 1.9x, 5.8x)
## multiplicative difference from a Cat 1-2
rbind(
%>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NFIP')) %>%
df.nfip %>%
rownames_to_column filter(cat == 2) %>%
select(data, cat, est)
## multiplicative difference from a Cat 2-3
rbind(
%>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NFIP')) %>%
df.nfip group_by(data) %>%
arrange(data, cat) %>%
mutate(est = est/est[cat==2]) %>%
filter(cat == 3) %>%
%>%
ungroup select(data, cat, est)
## multiplicative difference from a Cat 3-4
rbind(
%>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NFIP')) %>%
df.nfip group_by(data) %>%
arrange(data, cat) %>%
mutate(est = est/est[cat==3]) %>%
filter(cat == 4) %>%
%>%
ungroup select(data, cat, est)
## multiplicative difference from a Cat 4-5
rbind(
%>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NFIP')) %>%
df.nfip group_by(data) %>%
arrange(data, cat) %>%
mutate(est = est/est[cat==4]) %>%
filter(cat == 5) %>%
%>%
ungroup select(data, cat, est)
## data cat est
## 1 NCEI 2 2.284765
## 2 NFIP 2 3.455397
## # A tibble: 2 x 3
## data cat est
## <chr> <dbl> <dbl>
## 1 NCEI 3 6.44
## 2 NFIP 3 6.01
## # A tibble: 2 x 3
## data cat est
## <chr> <dbl> <dbl>
## 1 NCEI 4 3.08
## 2 NFIP 4 3.80
## # A tibble: 2 x 3
## data cat est
## <chr> <dbl> <dbl>
## 1 NCEI 5 3.90
## 2 NFIP 5 20.0
Effect of AR intensity category and sequences on loss
<-
g1 rbind(
%>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NFIP')) %>%
df.nfip mutate(cat = toNumber(cat)) %>%
ggplot() +
geom_hline(yintercept = 1, color = 'grey50') +
geom_point(aes(x = cat, y = est, group = data, color = data, shape = data), size = 1.5) +
geom_line(aes(x = cat, y = est, group = data, color = data), size = 0.75) +
geom_ribbon(
aes(x = cat, ymin = lower, ymax = upper, group = data, fill = data, alpha = data)) +
geom_line(
data = data.frame(x = 1:5, y = c(0.1,0.4,3,20,260)),
aes(x=x, y=y/y[1], linetype = 'Corringham et al. (14)')) +
scale_x_continuous(labels = paste0('AR', 1:5)) +
scale_y_log10('Loss Multiplier', labels = comma_format(suffix = 'x')) + logticks('l') +
coord_cartesian(ylim = c(0.3,5e3)) +
scale_color_manual(values = roma.colors[c(2,5)]) +
scale_fill_manual(values = roma.colors[c(2,5)]) +
scale_shape_manual(values = c(17,16)) +
scale_alpha_manual(values = c(0.35, 0.25)) +
scale_linetype_manual(values = 'dashed') +
labs(x = 'Rank Coefficient', color = 'Data', fill = 'Data', shape = 'Data', alpha = 'Data') +
theme(
panel.grid.major.y = element_line(size = 0.25),
legend.position = c(0.25,0.72),
legend.title = element_blank(),
legend.spacing = unit(-15, 'pt'),
legend.key.width = unit(0.25, 'in'),
legend.background = element_rect(fill = NA, color = NA))
<-
g2 rbind(
%>% filter(!is.na(seq)) %>% mutate(data = 'NCEI'),
df.ncei %>% filter(!is.na(seq)) %>% mutate(data = 'NFIP')) %>%
df.nfip ggplot() +
geom_hline(yintercept = 1, color = 'grey50') +
geom_col(aes(x = data, y = est, fill = data), width = 0.65, show.legend = FALSE) +
geom_errorbar(aes(x = data, ymin = lower, ymax = upper), size = 0.5, width = 0.1) +
scale_fill_manual(values = roma.colors[c(2,5)]) +
scale_y_log10(limits = c(0.3,5e3)) + logticks('l') +
labs(x = 'Sequence Coefficient') +
theme(
panel.grid.major.y = element_line(size = 0.25),
axis.title.y = element_blank(), axis.text.y = element_blank(),
axis.ticks.y = element_blank(), axis.line.y = element_blank())
+ g2 + plot_layout(nrow = 1, widths = c(2,1)) +
g1 plot_annotation(tag_levels = 'A') &
theme(
text = element_text(family = 'Segoe UI', size = 9),
plot.tag = element_text(family = 'Segoe UI', face = 'bold', size = 9),
plot.tag.position = c(0.15,0.95))
ggsave('_figures-impacts/fig5_magnitude.png', width = 7.25, height = 3.5, dpi = 500)
Relationship between sequence, category, and loss
%>%
catalog filter(ncei_damage > 0 | claims_value > 0) %>%
count(ar.cat, seq) %>%
pivot_wider(names_from = seq, values_from = n) %>%
mutate(across(-ar.cat, ~setNA(.x,0))) %>%
%>%
gt tab_header('Relationship between sequence, category, and loss') %>%
fmt_markdown(ar.cat) %>%
fmt_number(c(`FALSE`,`TRUE`), decimals = 0, use_seps = TRUE, sep_mark = ',') %>%
cols_label(
ar.cat = 'AR Category',
`FALSE` = html('Outside<br>Sequences'),
`TRUE` = html('Within<br>Sequences')) %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Relationship between sequence, category, and loss | ||
---|---|---|
AR Category | Outside Sequences |
Within Sequences |
1 |
446 | 341 |
2 |
168 | 451 |
3 |
36 | 493 |
4 |
1 | 129 |
5 |
0 | 48 |
cat('Total number of AR events with nonzero NCEI/NFIP loss: ')
cat(catalog %>% filter(ncei_damage > 0 | claims_value > 0) %>% nrow %>% comma)
cat('\n')
## Total number of AR events with nonzero NCEI/NFIP loss: 2,113
<-
coef.ncei foreach (b = 1:boot) %do% {
<- c(index.ncei.nonzero, sample(
index.sample size = 2*length(index.ncei.nonzero), replace = FALSE))
index.ncei.zero, felm(
log10(ncei_damage+1) ~ factor(ar.cat) + factor(prev.cat) | factor(id) + factor(wy),
data = catalog %>%
mutate(prev.cat = ifelse(prev.inter <= 5, prev.cat, 0)) %>%
slice(index.sample),
exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),mean)
<-
coef.nfip foreach (b = 1:boot) %do% {
<- c(index.nfip.nonzero, sample(
index.sample size = 2*length(index.nfip.nonzero), replace = FALSE))
index.nfip.zero, felm(
log10(claims_value+1) ~ factor(ar.cat) + factor(prev.cat) | factor(id) + factor(wy),
data = catalog %>%
mutate(prev.cat = ifelse(prev.inter <= 5, prev.cat, 0)) %>%
slice(index.sample),
exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),mean)
<- rbind(
coef.prevcat %>% data.frame %>% rownames_to_column %>%
coef.ncei slice(-(1:4)) %>%
transmute(
prev.cat = 1:5, source = 'NCEI',
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
rbind(data.frame(prev.cat = 0, source = 'NCEI', est = 1, lower = 1, upper = 1)),
%>% data.frame %>% rownames_to_column %>%
coef.nfip slice(-(1:4)) %>%
transmute(
prev.cat = 1:5, source = 'NFIP',
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
rbind(data.frame(prev.cat = 0, source = 'NFIP', est = 1, lower = 1, upper = 1)))
cat('Relationship between previous event coefficient and primary event coefficient\n')
cat(' NCEI: ')
cat(coef.prevcat %>% filter(prev.cat == 5 & source == 'NCEI') %>% pull(est) %>% round(1))
cat('x (')
cat(
%>% filter(prev.cat == 5 & source == 'NCEI') %>% pull(est)) /
((coef.prevcat %>% filter(cat == 5) %>% pull(est))) %>%
(df.ncei percent(., accuracy = 0.1))
cat(')\n')
cat(' NFIP: ')
cat(coef.prevcat %>% filter(prev.cat == 5 & source == 'NFIP') %>% pull(est) %>% round(1))
cat('x (')
cat(
%>% filter(prev.cat == 5 & source == 'NFIP') %>% pull(est)) /
((coef.prevcat %>% filter(cat == 5) %>% pull(est))) %>%
(df.nfip percent(., accuracy = 0.1))
cat(')\n')
## Relationship between previous event coefficient and primary event coefficient
## NCEI: 18.1x (10.2%)
## NFIP: 302.9x (19.2%)
<- c(0,1,3,5,7,14,30,Inf)
prev.bucket <-
coef.ncei foreach (b = 1:boot) %do% {
<- c(index.ncei.nonzero, sample(
index.sample size = 2*length(index.ncei.nonzero), replace = FALSE))
index.ncei.zero, felm(
log10(ncei_damage+1) ~ factor(ar.cat) + prev.bucket | factor(id) + factor(wy),
data = catalog %>%
mutate(prev.bucket = cut(
%>% factor(levels = rev(levels(.)))) %>%
prev.inter, prev.bucket) slice(index.sample),
exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),mean)
<-
coef.nfip foreach (b = 1:boot) %do% {
<- c(index.nfip.nonzero, sample(
index.sample size = 2*length(index.nfip.nonzero), replace = FALSE))
index.nfip.zero, felm(
log10(claims_value+1) ~ factor(ar.cat) + prev.bucket | factor(id) + factor(wy),
data = catalog %>%
mutate(prev.bucket = cut(
%>% factor(levels = rev(levels(.)))) %>%
prev.inter, prev.bucket) slice(index.sample),
exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),mean)
<-
coef.previnter rbind(
%>% data.frame %>% rownames_to_column %>%
coef.ncei slice(-(1:4)) %>%
transmute(
prev.bucket =
gsub('prev\\.bucket\\(','',rowname) %>% gsub('\\]','', .) %>% gsub(',','-',.),
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
rbind(data.frame(prev.bucket = '30+', est = 1, lower = 1, upper = 1)) %>%
mutate(source = 'NCEI'),
%>% data.frame %>% rownames_to_column %>%
coef.nfip slice(-(1:4)) %>%
transmute(
prev.bucket =
gsub('prev\\.bucket\\(','',rowname) %>% gsub('\\]','', .) %>% gsub(',','-',.),
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
rbind(data.frame(prev.bucket = '30+', est = 1, lower = 1, upper = 1)) %>%
mutate(source = 'NFIP')) %>%
mutate(prev.bucket = factor(
levels = c('0-1','1-3','3-5','5-7','7-14','14-30','30+'))) prev.bucket,
## calculate coefficients relative to no previous AR within 180 days
<- c(0,1,3,5,7,14,30,180,Inf)
prev.bucket <-
coef.ncei foreach (b = 1:boot) %do% {
<- c(index.ncei.nonzero, sample(
index.sample size = 2*length(index.ncei.nonzero), replace = FALSE))
index.ncei.zero, felm(
log10(ncei_damage+1) ~ factor(ar.cat) + prev.bucket | factor(id) + factor(wy),
data = catalog %>%
mutate(prev.bucket = cut(
%>% factor(levels = rev(levels(.)))) %>%
prev.inter, prev.bucket) slice(index.sample),
exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),mean)
<-
coef.nfip foreach (b = 1:boot) %do% {
<- c(index.nfip.nonzero, sample(
index.sample size = 2*length(index.nfip.nonzero), replace = FALSE))
index.nfip.zero, felm(
log10(claims_value+1) ~ factor(ar.cat) + prev.bucket | factor(id) + factor(wy),
data = catalog %>%
mutate(prev.bucket = cut(
%>% factor(levels = rev(levels(.)))) %>%
prev.inter, prev.bucket) slice(index.sample),
exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),mean)
.180 <-
coefrbind(
%>% data.frame %>% rownames_to_column %>%
coef.ncei slice(-(1:4)) %>%
transmute(
prev.bucket =
gsub('prev\\.bucket\\(','',rowname) %>% gsub('\\]','', .) %>% gsub(',','-',.),
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
rbind(data.frame(prev.bucket = '30+', est = 1, lower = 1, upper = 1)) %>%
mutate(source = 'NCEI'),
%>% data.frame %>% rownames_to_column %>%
coef.nfip slice(-(1:4)) %>%
transmute(
prev.bucket =
gsub('prev\\.bucket\\(','',rowname) %>% gsub('\\]','', .) %>% gsub(',','-',.),
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
rbind(data.frame(prev.bucket = '30+', est = 1, lower = 1, upper = 1)) %>%
mutate(source = 'NFIP')) %>%
mutate(prev.bucket = factor(
levels = c('0-1','1-3','3-5','5-7','7-14','14-30','30-180', '180+'))) prev.bucket,
## compare no AR within 30 days vs. 180 days
cat('Effect of an previous AR within 0-1 days, relative to no AR within 30 days\n')
cat(' NCEI: ')
cat(
%>% filter(prev.bucket == '0-1' & source == 'NCEI') %>%
coef.previnter pull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat('\n')
cat(' NCEI: ')
cat(
%>% filter(prev.bucket == '0-1' & source == 'NFIP') %>%
coef.previnter pull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat('\n\n')
cat('Effect of an previous AR within 0-1 days, relative to no AR within 180 days\n')
cat(' NCEI: ')
cat(
.180 %>% filter(prev.bucket == '0-1' & source == 'NCEI') %>%
coefpull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat('\n')
cat(' NCEI: ')
cat(
.180 %>% filter(prev.bucket == '0-1' & source == 'NFIP') %>%
coefpull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat('\n\n')
## Effect of an previous AR within 0-1 days, relative to no AR within 30 days
## NCEI: 13.1x
## NCEI: 12.6x
##
## Effect of an previous AR within 0-1 days, relative to no AR within 180 days
## NCEI: 31.5x
## NCEI: 27.1x
<- c(0,1,3,7,30,Inf)
prev.bucket
<-
coef.ncei foreach (b = 1:boot) %do% {
<- c(index.ncei.nonzero, sample(
index.sample size = 2*length(index.ncei.nonzero), replace = FALSE))
index.ncei.zero, felm(
log10(ncei_damage+1) ~ factor(ar.cat) + prev.combine | factor(id) + factor(wy),
data = catalog %>%
mutate(
prev.cat = ifelse(prev.cat==5, 4, prev.cat),
prev.bucket = cut(prev.inter, prev.bucket) %>% factor(levels = rev(levels(.))),
prev.combine =
ifelse(prev.inter > 30, 'over30', paste(prev.cat, prev.bucket, sep = '_')) %>%
factor(levels = rev(levels(factor(.))))) %>%
slice(index.sample),
exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),Mean)
<-
coef.nfip foreach (b = 1:boot) %do% {
<- c(index.nfip.nonzero, sample(
index.sample size = 2*length(index.nfip.nonzero), replace = FALSE))
index.nfip.zero, felm(
log10(claims_value+1) ~ factor(ar.cat) + prev.combine | factor(id) + factor(wy),
data = catalog %>%
mutate(
prev.cat = ifelse(prev.cat==5, 4, prev.cat),
prev.bucket = cut(prev.inter, prev.bucket) %>% factor(levels = rev(levels(.))),
prev.combine =
ifelse(prev.inter > 30, 'over30', paste(prev.cat, prev.bucket, sep = '_')) %>%
factor(levels = rev(levels(factor(.))))) %>%
slice(index.sample),
exactDOF = TRUE) %>%
%>% coef
summary %>%
} abind(along = 3) %>% apply(c(1,2),Mean)
<-
coef.prevboth rbind(
%>% data.frame %>% slice(-(1:4)) %>% rownames_to_column %>%
coef.ncei mutate(source = 'NCEI'),
%>% data.frame %>% slice(-(1:4)) %>% rownames_to_column %>%
coef.nfip mutate(source = 'NFIP')) %>%
separate(rowname, into = c('prev.cat', 'prev.bucket'), sep = '_') %>%
mutate(
prev.cat = gsub('prev\\.combine','', prev.cat) %>% toNumber,
prev.bucket = gsub('\\(','',prev.bucket) %>% gsub('\\]','',.) %>% gsub(',','-',.)) %>%
transmute(
prev.cat, prev.bucket, source,est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error))
cat('Regression coefficient for Cat 4-5 events within the last 0-1 days\n')
cat(' NCEI: ')
cat(
%>%
coef.prevboth group_by(source) %>%
filter(est == max(est) & source == 'NCEI') %>%
pull(est) %>% comma(accuracy = 1, suffix = 'x'))
cat('\n')
cat(' NFIP: ')
cat(
%>%
coef.prevboth group_by(source) %>%
filter(est == max(est) & source == 'NFIP') %>%
pull(est) %>% comma(accuracy = 1, suffix = 'x'))
cat('\n')
## Regression coefficient for Cat 4-5 events within the last 0-1 days
## NCEI: 252x
## NFIP: 171x
Effect of previous AR characteristics on loss
<- ggplot(coef.prevcat) +
g1 geom_hline(yintercept = 1, color = 'grey50') +
geom_point(aes(x = prev.cat, y = est, group = source, color = source, shape = source), size = 1.5) +
geom_line(aes(x = prev.cat, y = est, group = source, color = source), size = 0.75) +
geom_ribbon(aes(x = prev.cat, ymin = lower, ymax = upper, group = source, fill = source, alpha = source)) +
scale_x_continuous(breaks = 0:5, labels = c('No Previous\nEvent', paste0('AR',1:5))) +
scale_y_log10(labels = comma_format(suffix = 'x')) + logticks('l') +
scale_color_manual('', values = roma.colors[c(2,5)]) +
scale_fill_manual('', values = roma.colors[c(2,5)]) +
scale_shape_manual('', values = c(17,16)) +
scale_alpha_manual('', values = c(0.35, 0.25)) +
ggtitle(waiver(), subtitle = '') +
labs(x = 'Previous Event Rank', y = 'Primary Event Loss Multiplier') +
theme(
text = element_text(family = 'Segoe UI', size = 9),
panel.grid.major.y = element_line(size = 0.25),
legend.key.width = unit(0.25, 'in'),
legend.position = c(0.2,0.9))
<- ggplot(coef.previnter) +
g2 geom_hline(yintercept = 1, color = 'grey50') +
geom_point(aes(x = prev.bucket, y = est, group = source, color = source, shape = source), size = 1.5) +
geom_line(aes(x = prev.bucket, y = est, group = source, color = source), size = 0.75) +
geom_ribbon(aes(x = prev.bucket, ymin = lower, ymax = upper, group = source, fill = source, alpha = source)) +
scale_y_log10(labels = comma_format(suffix = 'x')) + logticks('l') +
scale_color_manual('', values = roma.colors[c(2,5)]) +
scale_fill_manual('', values = roma.colors[c(2,5)]) +
scale_shape_manual('', values = c(17,16)) +
scale_alpha_manual('', values = c(0.35, 0.25)) +
ggtitle(waiver(), subtitle = '') +
labs(x = 'Days Since Previous AR', y = 'Primary Event Loss Multiplier') +
theme(
text = element_text(family = 'Segoe UI', size = 9),
panel.grid.major.y = element_line(size = 0.25),
legend.key.width = unit(0.25, 'in'),
legend.position = c(0.8,0.9))
<- ggplot(coef.prevboth %>% mutate(tag = case_match(source, 'NCEI'~'C', 'NFIP'~'D'))) +
g3 geom_hline(yintercept = 1, color = 'grey50') +
geom_point(
aes(x = prev.cat, y = est, group = prev.bucket, color = prev.bucket),
position = position_dodge(width = 0.5), size = 1.5) +
geom_linerange(
aes(x = prev.cat, ymin = lower, ymax = upper, group = prev.bucket, color = prev.bucket),
position = position_dodge(width = 0.5), size = 0.75) +
facet_wrap(~source, nrow = 1) +
geom_text(
data = data.frame(
x = c(0.6,0.6), y = c(1e4,1e4),
lab = c('C', 'D'), source = c('NCEI','NFIP')),
aes(x=x, y=y, label = lab),
family = 'Segoe UI', size = 9/.pt, fontface = 'bold') +
labs(x = 'Previous Event Rank', y = 'Primary Event Loss Multiplier', color = 'Days Since Previous AR') +
scale_x_continuous(breaks = 1:4, labels = paste0('AR',c(1:3,'4-5'), '\n')) +
scale_y_log10(breaks = 10^(-1:3), labels = paste0(c(0.1, 1, 10, 100, '1,000'), 'x')) +
logticks('l') +
scale_color_scico_d(palette = 'davos', end = 0.8, guide = guide_legend(nrow = 1)) +
coord_cartesian(ylim = c(NA,1e3), clip = 'off') +
theme(
text = element_text(family = 'Segoe UI', size = 9),
plot.margin = margin(20,0,0,0),
# plot.background = element_rect(color = 'pink'),
legend.position = 'bottom',
legend.key.width = unit(0.25, 'in'),
panel.grid.major.y = element_line(size = 0.25),
strip.text = element_text(margin = margin(2,2,2,2)),
legend.margin = margin(0,0,0,0),
strip.background = element_rect(color = NA, fill = 'grey95'))
g3
+ g2 + g3 +
g1 plot_layout(design = 'ab\ncc', heights = c(6,5)) +
plot_annotation(
tag_levels = list(c('A', 'B', ''))) &
theme(
plot.tag = element_text(family = 'Segoe UI', face = 'bold', size = 9),
plot.tag.position = c(0.15, 1))
ggsave('_figures-impacts/fig6_previous.png', width = 7.25, height = 5.75, dpi = 500)
## add sequence position variable
<- catalog %>%
catalog filter(seq) %>%
select(id, ar.count, seq.count) %>%
group_by(id, seq.count) %>%
mutate(seq.position = 1:length(unique(ar.count))) %>%
%>%
ungroup left_join(catalog, ., by = c('id', 'ar.count', 'seq.count'))
<- catalog %>% mutate(seq.position = ifelse(seq.position>10, 10, seq.position)) catalog
cat('Number of AR events within sequences: ')
cat(comma(sum(catalog$seq)))
cat('\n')
cat('Percent of AR events within sequences: ')
cat((sum(catalog$seq)/nrow(catalog)) %>% percent(accuracy = 0.1))
cat('\n')
## Number of AR events within sequences: 12,745
## Percent of AR events within sequences: 54.8%
cat('Percentage of AR events in the 10th sequence position or later: ')
cat(
%>%
catalog filter(seq) %>%
count(seq.position) %>%
mutate(pct = n/sum(n)) %>%
filter(seq.position >= 10) %>%
pull(pct) %>% sum %>% percent(accuracy = 0.1))
cat('\n')
## Percentage of AR events in the 10th sequence position or later: 1.5%
cat('Average number of days from sequence start to the 5th AR: ')
cat(
%>%
catalog filter(seq) %>%
filter(seq.position == 5) %>%
mutate(delay = as.numeric(as.Date(ar.start) - as.Date(seq.start))) %>%
pull(delay) %>% median)
cat('\n')
## Average number of days from sequence start to the 5th AR: 25
<- rbind(
coef.position felm(
log10(ncei_damage+1) ~ factor(seq.position) | factor(id)*factor(wy),
data = catalog %>% filter(seq)) %>%
%>% coef %>% data.frame %>%
summary rownames_to_column('var') %>%
rbind(setNames(data.frame('factor(seq.position)1', 0,0,NA,NA), names(.)), .) %>%
mutate(
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
transmute(var, est, lower, upper, level = 1:nrow(.), source = 'NCEI'),
felm(
log10(claims_value+1) ~ factor(seq.position) | factor(id)*factor(wy),
data = catalog %>% filter(seq)) %>%
%>% coef %>% data.frame %>%
summary rownames_to_column('var') %>%
rbind(setNames(data.frame('factor(seq.position)1', 0,0,NA,NA), names(.)), .) %>%
mutate(
est = 10^Estimate,
lower = 10^(Estimate-2*Std..Error),
upper = 10^(Estimate+2*Std..Error)) %>%
transmute(var, est, lower, upper, level = 1:nrow(.), source = 'NFIP'))
cat('Loss multiplier for the 5th sequence position\n')
cat(' NCEI: ')
cat(
%>%
coef.position group_by(source) %>%
filter(est == max(est) & source == 'NCEI') %>%
pull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat('\n')
cat(' NFIP: ')
cat(
%>%
coef.position group_by(source) %>%
filter(est == max(est) & source == 'NFIP') %>%
pull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat('\n')
## Loss multiplier for the 5th sequence position
## NCEI: 2.8x
## NFIP: 1.6x
Effect of AR sequence position on loss
<- coef.position %>%
g1 ggplot() +
geom_hline(yintercept = 1, color = 'grey50') +
geom_point(aes(x = level, y = est, color = source, shape = source), size = 1.5) +
geom_line(aes(x = level, y = est, color = source), size = 0.75) +
geom_ribbon(aes(x = level, ymin = lower, ymax = upper, fill = source, alpha = source)) +
labs(x = 'Sequence Position', y = 'Loss Multiplier', color = 'Data', fill = 'Data') +
scale_x_continuous(breaks = 1:10, labels = c(1:9, '10+')) +
scale_y_log10(labels = comma_format(suffix = 'x')) + logticks('l') +
scale_color_manual('', values = roma.colors[c(2,5)]) +
scale_fill_manual('', values = roma.colors[c(2,5)]) +
scale_shape_manual('', values = c(17,16)) +
scale_alpha_manual('', values = c(0.35, 0.25)) +
theme(
text = element_text(family = 'Segoe UI', size = 9),
legend.key.width = unit(0.25, 'in'),
legend.background = element_rect(fill = NA),
legend.position = c(0.1,0.8),
panel.grid.major.y = element_line(size = 0.25),
plot.margin = margin(10,5,5,5))
<- catalog %>%
g2 filter(seq) %>%
mutate(delay = as.numeric(as.Date(ar.start) - as.Date(seq.start))) %>%
ggplot() +
geom_hline(yintercept = 0) +
geom_histogram(aes(x = delay), color = NA, fill = 'grey40', size = 0.25, boundary = 0) +
facet_wrap(
~seq.position,
labeller = labeller(
seq.position = c(paste('N =', 1:9), 'N = 10+') %>% setNames(1:10)),
scales = 'free_x', nrow = 1) +
scale_x_origin() + scale_y_origin() +
coord_flip() +
labs(
x = expression(paste('Days to Reach ', N^th, ' Position')),
y = 'Number of Occurrences') +
theme(
text = element_text(family = 'Segoe UI', size = 9),
strip.background = element_rect(fill = 'grey95', color = NA),
strip.text = element_text(margin = margin(2,5,2,5)),
panel.grid.major = element_line(size = 0.25),
panel.spacing.x = unit(10, 'pt'),
plot.margin = margin(10,5,5,5),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
+ g2 + plot_layout(nrow = 2, heights = c(3,2)) +
g1 plot_annotation(tag_levels = 'A') &
theme(
plot.tag = element_text(family = 'Segoe UI', size = 9, face = 'bold'),
plot.tag.position = c(0.075,1.075))
ggsave('_figures-impacts/fig7_position.png', width = 7.25, height = 4, dpi = 500)
## how many are there of each position?
%>% filter(seq) %>% count(seq.position) catalog
## # A tibble: 10 x 2
## seq.position n
## <dbl> <int>
## 1 1 5056
## 2 2 3336
## 3 3 1850
## 4 4 1028
## 5 5 546
## 6 6 310
## 7 7 195
## 8 8 140
## 9 9 98
## 10 10 186