In ericward-noaa/pinkescape: Simulate Pink Salmon Escapement

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.asp = 0.618
)

library(ggplot2)
library(dplyr)
library(pinkescape)

Comparing variability of recruitment on net benefits

We will simulate a low and high variance scenario, calculate net benefits for each simulation and plot results.

df_lo = sim(rec_std = 0.05)
df_hi = sim(rec_std = 0.2)

# summarize total benefits
df_lo_summary = mutate(df_lo,
                              net_ben = discount*net_benefits) %>%
  group_by(sim) %>%
  summarise(tot_ben = sum(net_ben))

df_hi_summary = mutate(df_hi,
                              net_ben = discount*net_benefits) %>%
  group_by(sim) %>%
  summarise(tot_ben = sum(net_ben))

df_lo_summary$sd = 0.05
df_hi_summary$sd = 0.2
ggplot(rbind(df_lo_summary,df_hi_summary), aes(sd, tot_ben,group=sd)) +
  geom_boxplot()

Calculating means

Each simulation in a particular run uses the same set of random number seeds to model environmental variation, etc. We can control for this random variability by comparing the pairwise simulations (e.g. simulation # 1 from the low recruitment run to simulation # 1 from the high recruitment run) and calculate the distribution of differences across simulations.

Using the example above, df_lo_summary and df_hi_summary contain the total benefits for each simulation.

df_diff <- data.frame(y = df_hi_summary$tot_ben - df_lo_summary$tot_ben)

ggplot(df_diff, aes(y)) + 
  geom_histogram()

m <- dplyr::summarise(df_diff, 
                      mean = mean(y),
                      median = median(y),
                      lo25 = quantile(y, 0.25),
                      hi25 = quantile(y, 0.75),
                      lo_95 = quantile(y, 0.025),
                      hi_95 = quantile(y, 0.975))
knitr::kable(m, caption = "Summary statistics on differences")

Comparing effect of recruitment autocorrelation on net benefits

Similar to above, we can compare a white noise to autocorrelated scenario,

df_lo = sim(rec_acf = 0)
df_hi = sim(rec_acf = 0.7)
# summarize
df_lo_summary = mutate(df_lo,
                              net_ben = discount*net_benefits) %>%
  group_by(sim) %>%
  summarise(tot_ben = sum(net_ben))
df_hi_summary = mutate(df_hi,
                              net_ben = discount*net_benefits) %>%
  group_by(sim) %>%
  summarise(tot_ben = sum(net_ben))
df_lo_summary$acf = 0
df_hi_summary$acf = 0.7
ggplot(rbind(df_lo_summary,df_hi_summary), aes(acf, tot_ben,group=acf)) +
  geom_boxplot()

Adding harvest variability

We can introduce variability in realized harvest by specifying the harvest_CV parameter, which is a lognormal variable multiplied by expected harvest. Realistic values are probably in the 0.5 - 1.0 range.

sim(harvest_CV = 0.5)

Using custom Ricker parameters

Instead of the defaults, we can create a dataframe with custom parameters and pass these in. The default parameters can be found using the get_ricker() function,

pars <- get_ricker()

But let's change the Ricker $a$ and $b$ parameters. We have a custom

pars$a = c(1,1,2,2)
pars$b = runif(4)

These can be passed into the simulation now

custom_sims <- sim(ricker_pars = pars)

Adding MSY rules

We also allow for the decision rules to be altered. By default, the simulations use a constant escapement rule ("equilibrium"). The full arguments for the simulation function (with defaults) are below:

sim(sims = 1000, # number of simulations
                time_steps = 100,
                ricker_pars = NULL,
                pr_12 = 0.1,# probability transitioning 1:2
                pr_21 = 0.1,# probability transitioning 2:1
                run = "odd",
                deterministic_model = TRUE,
                rec_std = 0.1, # recruitment standard deviation in log space
                rec_acf = 0.7, # recruitment autocorrelation in time
                escapement_rule = "both",
                harvest_CV = 0,
                discount_rate = 0.1,
                price_linear_change = 0,
                price_acf = 0.7, # default based on PWS pink salmon 1984 - 2021
                price_cv = 0, # default based on PWS pink salmon 1984 - 2021
                time_lag = 0,
                msy_scenario = "equilibrium")

When the msy_scenario is equilibrium, then harvest at each time step is the difference between the returning abundance and optimal spawning stock size (can be adjusted by manipulating ricker_pars$S_star).

Alternatively, we can change the msy_scenario to msy. For this scenario, in each time step, we find the optimal spawning size S that maximizes the Ricker function. For both the deterministic and stochastic versions of the simulation -- without any price feedbacks, this function is

(1 + ricker_pars$b[x[t - time_lag]]*S) * exp(ricker_pars$a[x[t - time_lag]] + ricker_pars$b[x[t - time_lag]]*S)

where x represents the discrete regime assigned to each time step, t is the year of interest, time_lag represents and optional time lag, and a and b are Ricker parameters.

ericward-noaa/pinkescape documentation built on Sept. 21, 2023, 7:19 a.m.