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.

Setup

Functions & packages

source('_data/setup_impacts.R')
source('_scripts/create_df_functions.R')

## set number of bootstrapped samples
boot <- 1000

## turn progress bars on/off
progress <- FALSE
if (progress) pb <- txtProgressBar(min = 0, max = ncell(grid_ca), style = 3)

## set random seed for reproducibility
set.seed(2023)

Load data

load('_scripts/_checkpoints/df_3hr_1209.Rdata')
load('_scripts/_checkpoints/df_24hr_0209.Rdata')

Create AR event catalog

Add inter-event time to df_3hr

df_3hr <- 
  foreach (i = 1:ncell(grid_ca)) %do% {
    if (progress) setTxtProgressBar(pb, i)
    if (i %in% index_ca) {
      temp <- 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)) %>% 
        mutate(inter = !ar, inter.count = add_counter(inter))
      temp <- create_catalog(temp, 'inter', cat = FALSE, interval = 3/24) %>% 
        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)
      df_3hr[[i]] %>% left_join(temp, by = 'ar.count')
    } else NULL
  }

Create df_catalog

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
  }

Combine df_catalog into catalog

## keep WY 1981-2021 (for category analysis)
catalog.full <- df_catalog %>% 
  reduce(rbind) %>% 
  mutate(wy = wateryear(ar.start)) %>%
  filter(wateryear(ar.start) %in% 1981:2021)
  
## keep WY 1997-2021 (for impact analysis)
catalog <- catalog.full %>%
  filter(wateryear(ar.start) %in% 1997:2021) %>%
  mutate(
    ncei_damage = setNA(ncei_damage,0), 
    claims_value = setNA(claims_value,0)) 

Report catalog size

cat('Number of AR events in the catalog: ')
cat(comma(nrow(catalog)))
## Number of AR events in the catalog: 23,267

Report class imbalance

cat('Percent of events with nonzero NCEI loss: ')
cat(
  catalog %>% count(loss = ncei_damage > 0) %>% mutate(pct = n/sum(n)) %>% 
    filter(loss) %>% pull(pct) %>% percent(accuracy = 0.1))
cat('\n')

cat('Percent of events with nonzero NFIP loss: ')
cat(
  catalog %>% count(loss = claims_value > 0) %>% mutate(pct = n/sum(n)) %>% 
    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%

Result #1: Temporal compounding and AR category

Figure 2

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'))
g.or <- map(
  .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))))
g.or[[1]] <- g.or[[1]] + 
  facet_grid('Adjacent ARs' ~ 'AR1', switch = 'y')
g.or[[2]] <- g.or[[2]] +
  facet_grid(.~paste0('AR', ar.cat)) + 
  theme(axis.text.y = element_blank())
g.or[[3]] <- g.or[[3]] + 
  facet_grid(.~paste0('AR', ar.cat)) + 
  theme(axis.text.y = element_blank())
g.or[[4]] <- g.or[[4]] + 
  facet_grid(.~paste0('AR', ar.cat)) + 
  theme(axis.text.y = element_blank())
g.or[[5]] <- g.or[[5]] + 
  facet_grid(.~paste0('AR', ar.cat)) + 
  theme(axis.text.y = element_blank())
g.and <- map(
  .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))))
g.and[[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())
plot.or <- wrap_plots(g.or) + guide_area() + plot_layout(design = 'abcdef', guides = 'collect')
plot.and <- wrap_plots(g.and) + guide_area() + plot_layout(design = 'abcdef', guides = 'collect')

plot.or / plot.and
ggsave('_figures-impacts/fig2_coincidence.png', width=7.25, height=3.25, units = 'in', dpi = 500)

Report bubble size

catalog.full %>% count(ar.cat)
## # 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

Report statistics

cat('How much more likely are Category 5 events to be adjacent, compared to Category 1?\n')
adjacent <- catalog.full %>%  
  count(ar.cat, five = prev.inter <= 5 | next.inter <= 5) %>% 
  group_by(ar.cat) %>% 
  mutate(pct = n/sum(n)) %>% 
  filter(five)
(adjacent$pct[5] / adjacent$pct[1]) %>% comma(accuracy = 0.1, suffix = 'x')
cat('\n')

cat('How much more likely are Category 5 events to be sandwiched, compared to Category 1?\n')
sandwich <- catalog.full %>%  
  count(ar.cat, five = prev.inter <= 5 & next.inter <= 5) %>% 
  group_by(ar.cat) %>% 
  mutate(pct = n/sum(n)) %>% 
  filter(five)
(sandwich$pct[5] / sandwich$pct[1]) %>% comma(accuracy = 0.1, suffix = 'x')
## 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"

Figure 3

Relationship between primary event category and preceding event category

Bootstrap confidence intervals for “All AR Events”

catalog5 <- catalog.full %>% filter(prev.inter <= 5)

## 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(
      catalog5$ar.cat, size = nrow(catalog5), replace = TRUE) %>% 
      table %>% prop.table %>% cumsum) %>% 
  apply(2, function(x) quantile(x, c(0.05,0.5,0.95))) %>% 
  t %>% as.data.frame %>% 
  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

Generate figure

g.left <- ggplot(ratios) + 
  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))

g1 <- g.left + 
  scale_y_continuous(labels = percent, expand = expansion(mult = 0)) + 
  coord_cartesian(xlim = c(0.95,1.35), ylim = c(0,0.85))
g2 <- g.left + 
  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())
g.right <- catalog5 %>% 
  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))

g3 <- g.right + 
  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))
g4 <- g.right + 
  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 <- 
  cowplot::get_plot_component(ggplot() + labs(y = '          Percent of Occurrences'), 'ylab-l')

g2 + g1 + g4 + g3 + plot_spacer() + guide_area() + plot_spacer() + gg_axis + plot_spacer() +
  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)

Report statistics

temp <- catalog5 %>% 
  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"

Result #2: Sequences and likelihood of impact

Create daily dataframe

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

Report statistics

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%"

Fit regression

f2 <- felm(
  impactday ~ ar*seq - seq | factor(id) + factor(wy), 
  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

Figure 4

Impact day likelihood regression coefficients

summary(f2, robust = TRUE) %>% 
  coef %>% data.frame %>%
  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)

Result #3: Sequences and magnitude of impact

Fit regression

## NCEI
index.ncei.zero <- which(catalog$ncei_damage == 0)
index.ncei.nonzero <- which(catalog$ncei_damage > 0)
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(
          index.ncei.zero, size = 2*length(index.ncei.nonzero), replace = FALSE))),
      exactDOF = TRUE) %>%
      summary %>% coef
  } %>%
  abind(along = 3) %>% apply(c(1,2),mean)
df.ncei <- coef.under %>% data.frame %>%
  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
index.nfip.zero <- which(catalog$claims_value == 0)
index.nfip.nonzero <- which(catalog$claims_value > 0)
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(
          index.nfip.zero, size = 2*length(index.nfip.nonzero), replace = FALSE))),
      exactDOF = TRUE) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),mean)
df.nfip <- coef.under %>% data.frame %>% 
  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))

Report coefficients

df.nfip %>% filter(is.na(seq)) %>% 
  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(
  df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
  df.nfip %>% filter(is.na(seq)) %>% mutate(data = '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(
  df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
  df.nfip %>% filter(is.na(seq)) %>% mutate(data = 'NFIP')) %>%
  rownames_to_column %>% 
  filter(cat == 2) %>% 
  select(data, cat, est)

## multiplicative difference from a Cat 2-3 
rbind(
  df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
  df.nfip %>% filter(is.na(seq)) %>% mutate(data = '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(
  df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
  df.nfip %>% filter(is.na(seq)) %>% mutate(data = '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(
  df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
  df.nfip %>% filter(is.na(seq)) %>% mutate(data = '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

Figure 5

Effect of AR intensity category and sequences on loss

g1 <- 
  rbind(
    df.ncei %>% filter(is.na(seq)) %>% mutate(data = 'NCEI'),
    df.nfip %>% filter(is.na(seq)) %>% mutate(data = '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(
    df.ncei %>% filter(!is.na(seq)) %>% mutate(data = 'NCEI'),
    df.nfip %>% filter(!is.na(seq)) %>% mutate(data = '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())

g1 + g2 + plot_layout(nrow = 1, widths = c(2,1)) + 
  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)

Table 1

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

Result #4a: Effect of previous event characteristics on loss

Previous event category

Fit regression

coef.ncei <- 
  foreach (b = 1:boot) %do% {
    index.sample <- c(index.ncei.nonzero, sample(
      index.ncei.zero, size = 2*length(index.ncei.nonzero), replace = FALSE))
    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) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),mean)
coef.nfip <- 
  foreach (b = 1:boot) %do% {
    index.sample <- c(index.nfip.nonzero, sample(
      index.nfip.zero, size = 2*length(index.nfip.nonzero), replace = FALSE))
    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) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),mean)

coef.prevcat <- rbind(
  coef.ncei %>% data.frame %>% rownames_to_column %>% 
    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)),
  coef.nfip %>% data.frame %>% rownames_to_column %>% 
    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)))

Report statistics

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(
  ((coef.prevcat %>% filter(prev.cat == 5 & source == 'NCEI') %>% pull(est)) / 
  (df.ncei %>% filter(cat == 5) %>% pull(est))) %>% 
  percent(., accuracy = 0.1))
cat(')\n')

cat('  NFIP: ')
cat(coef.prevcat %>% filter(prev.cat == 5 & source == 'NFIP') %>% pull(est) %>% round(1))
cat('x (')
cat(
  ((coef.prevcat %>% filter(prev.cat == 5 & source == 'NFIP') %>% pull(est)) / 
  (df.nfip %>% filter(cat == 5) %>% pull(est))) %>% 
  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%)

Between-event interval

Fit regression

prev.bucket <- c(0,1,3,5,7,14,30,Inf)
coef.ncei <- 
  foreach (b = 1:boot) %do% {
    index.sample <- c(index.ncei.nonzero, sample(
      index.ncei.zero, size = 2*length(index.ncei.nonzero), replace = FALSE))
    felm(
      log10(ncei_damage+1) ~ factor(ar.cat) + prev.bucket | factor(id) + factor(wy),
      data = catalog %>% 
        mutate(prev.bucket = cut(
          prev.inter, prev.bucket) %>% factor(levels = rev(levels(.)))) %>%
        slice(index.sample), 
      exactDOF = TRUE) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),mean)
coef.nfip <- 
  foreach (b = 1:boot) %do% {
    index.sample <- c(index.nfip.nonzero, sample(
      index.nfip.zero, size = 2*length(index.nfip.nonzero), replace = FALSE))
    felm(
      log10(claims_value+1) ~ factor(ar.cat) + prev.bucket | factor(id) + factor(wy),
      data = catalog %>% 
        mutate(prev.bucket = cut(
          prev.inter, prev.bucket) %>% factor(levels = rev(levels(.)))) %>%
        slice(index.sample), 
      exactDOF = TRUE) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),mean)

coef.previnter <- 
  rbind(
    coef.ncei %>% data.frame %>% rownames_to_column %>% 
      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'),
    coef.nfip %>% data.frame %>% rownames_to_column %>% 
      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(
    prev.bucket, levels = c('0-1','1-3','3-5','5-7','7-14','14-30','30+')))

Report statistics

## calculate coefficients relative to no previous AR within 180 days
prev.bucket <- c(0,1,3,5,7,14,30,180,Inf)
coef.ncei <- 
  foreach (b = 1:boot) %do% {
    index.sample <- c(index.ncei.nonzero, sample(
      index.ncei.zero, size = 2*length(index.ncei.nonzero), replace = FALSE))
    felm(
      log10(ncei_damage+1) ~ factor(ar.cat) + prev.bucket | factor(id) + factor(wy),
      data = catalog %>% 
        mutate(prev.bucket = cut(
          prev.inter, prev.bucket) %>% factor(levels = rev(levels(.)))) %>%
        slice(index.sample), 
      exactDOF = TRUE) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),mean)
coef.nfip <- 
  foreach (b = 1:boot) %do% {
    index.sample <- c(index.nfip.nonzero, sample(
      index.nfip.zero, size = 2*length(index.nfip.nonzero), replace = FALSE))
    felm(
      log10(claims_value+1) ~ factor(ar.cat) + prev.bucket | factor(id) + factor(wy),
      data = catalog %>% 
        mutate(prev.bucket = cut(
          prev.inter, prev.bucket) %>% factor(levels = rev(levels(.)))) %>%
        slice(index.sample), 
      exactDOF = TRUE) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),mean)

coef.180 <- 
  rbind(
    coef.ncei %>% data.frame %>% rownames_to_column %>% 
      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'),
    coef.nfip %>% data.frame %>% rownames_to_column %>% 
      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(
    prev.bucket, levels = c('0-1','1-3','3-5','5-7','7-14','14-30','30-180', '180+')))
## 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(
  coef.previnter %>% filter(prev.bucket == '0-1' & source == 'NCEI') %>% 
    pull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat('\n')
cat('  NCEI: ')
cat(
  coef.previnter %>% filter(prev.bucket == '0-1' & source == 'NFIP') %>% 
    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(
  coef.180 %>% filter(prev.bucket == '0-1' & source == 'NCEI') %>% 
    pull(est) %>% comma(accuracy = 0.1, suffix = 'x'))
cat('\n')
cat('  NCEI: ')
cat(
  coef.180 %>% filter(prev.bucket == '0-1' & source == 'NFIP') %>% 
    pull(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

Interactions between previous event characteristics

Fit regression

prev.bucket <- c(0,1,3,7,30,Inf)

coef.ncei <- 
  foreach (b = 1:boot) %do% {
    index.sample <- c(index.ncei.nonzero, sample(
      index.ncei.zero, size = 2*length(index.ncei.nonzero), replace = FALSE))
    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) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),Mean)
coef.nfip <- 
  foreach (b = 1:boot) %do% {
    index.sample <- c(index.nfip.nonzero, sample(
      index.nfip.zero, size = 2*length(index.nfip.nonzero), replace = FALSE))
    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) %>% 
      summary %>% coef
  } %>% 
  abind(along = 3) %>% apply(c(1,2),Mean)

coef.prevboth <- 
  rbind(
    coef.ncei %>% data.frame %>% slice(-(1:4)) %>% rownames_to_column %>% 
      mutate(source = 'NCEI'),
    coef.nfip %>%  data.frame %>% slice(-(1:4)) %>% rownames_to_column %>% 
      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))

Report statistics

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

Figure 6

Effect of previous AR characteristics on loss

g1 <- ggplot(coef.prevcat) +
  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))

g2 <- ggplot(coef.previnter) + 
  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))

g3 <- ggplot(coef.prevboth %>% mutate(tag = case_match(source, 'NCEI'~'C', 'NFIP'~'D'))) +
  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

g1 + g2 + g3 + 
  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)

Result #4b: Effect of sequence position on loss

## 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 <- catalog %>% mutate(seq.position = ifelse(seq.position>10, 10, seq.position))

Report statistics

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

Fit regression

coef.position <- rbind(
  felm(
    log10(ncei_damage+1) ~ factor(seq.position) | factor(id)*factor(wy), 
    data = catalog %>% filter(seq)) %>% 
    summary %>% coef %>% data.frame %>% 
    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)) %>% 
    summary %>% coef %>% data.frame %>% 
    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

Figure 7

Effect of AR sequence position on loss

g1 <- coef.position %>% 
  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))

g2 <- catalog %>% 
  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))

g1 + g2 + plot_layout(nrow = 2, heights = c(3,2)) + 
  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? 

catalog %>% filter(seq) %>% count(seq.position)
## # 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