Simulate the system forward — ms

Annual Schaefer growth on start-of-year biomass; catches removed monthly (Baranov-type, C = B (1 - exp(-F))). Each month every segment allocates its effort among metiers by expected profit computed from the biomass it perceives (info: "current" = true current biomass, "lagged" = last year's realised CPUE by metier).

Usage

ms_simulate(
  sys,
  caps,
  years = 100,
  gamma = 1.5,
  behaviour = c("dynamic", "fixed"),
  info = c("current", "lagged"),
  phi = 0.5,
  open = 1,
  hcr = NULL,
  hcr_lag = 2,
  B_init = NULL,
  E_init = NULL,
  capacity = NULL,
  cost_seg = NULL,
  proc_cv = 0,
  rho = 0,
  proc_r_power = 1,
  proc_rho_sp = 0,
  proc_sd = NULL,
  proc_dep = c("none", "jabba", "power"),
  proc_dep_pars = c(a = 0.128, c = 0.223, d = 1.326),
  proc_dep_power = 0,
  proc_dep_ref = 0.5,
  proc_sd_max = 1,
  proc_center = c("mean", "median"),
  impl_cv = 0,
  impl_bias = 0,
  effort_cv = 0,
  effort_rho = 0,
  seed = NULL,
  crn = NULL,
  keep = c("final", "all"),
  avg_years = 10,
  F_cap = sys$F_cap %||% Inf,
  F_cap_p = NULL
)

Arguments

sys: an ms_system.
caps: n x G matrix of annual effort caps by segment
years: number of years.
gamma: preference exponent (0 = profit-blind)
behaviour: "dynamic" (reallocate every month) or "fixed" (shares frozen at the unfished allocation: the classical shared-q model)
info: expected value from the current biomass or last year's catch rates.
phi: entry/exit rate of each segment: E_t+1 = min(cap, E_t * exp(phi * margin_g)), margin_g = segment profit / segment revenue; 0 = effort always at the cap
open: métier access, 0-1 (the open share of each métier's ground or season): a vector of length M (all months), a 12 x M matrix (by month, e.g. a seasonal closure), or a function of the year returning either.
hcr: optional function(state) -> n x G multiplier on caps
hcr_lag: years between the data a rule sees and the year its advice applies (default 2: data of year y - 2, assessed in y - 1, advice for y). Before the projection the data years are taken to be the initial state.
B_init: initial biomass (n x S, or one value per stock); default sys$dep_init * K (the system's starting depletion, 1 if not set)
E_init: initial effort by segment (vector of length G or n x G; default the caps)
capacity: fleet days by segment (registered vessels x days); default the caps. With vessel heterogeneity (sys$vessel_cv > 0) and phi > 0, the days used move towards capacity x the share of vessels that break even (lognormal vessel efficiency, mean 1), and the active vessels' mean efficiency is the mean of that upper tail. phi is then the annual adjustment rate (0-1).
cost_seg: daily cost by segment (length G or n x G), overrides sys$cost_seg
proc_cv: process error CV at the median r (JABBA form: the whole process equation is multiplied by exp(eta), eta ~ N(-sigma^2/2, sigma^2))
rho: AR(1) autocorrelation of the process error in time
proc_r_power: the process CV of species i is proc_cv * (r_i / median r)^power: more productive species are more variable (0 = same CV for all)
proc_rho_sp: correlation of the process error across species (shared environment)
proc_sd: log-scale process SD by species (length S), overrides proc_cv and proc_r_power; e.g. proc_sd_from_r(sys$r) (0.05 for sharks, 0.4 for anchovy)
proc_dep: depletion-dependent process error: "jabba" = JABBA proc.dyn (jabba2jags.R): sd_t = sd_i * relsig(B/Bmsy), relsig(x) = raw(x) / raw(1), raw(x) = a + (1 - a) / (1 + (x / c)^d), empirical ICES a = 0.128, c = 0.223, d = 1.326 (proc_dep_pars); relsig = 1 at Bmsy, ~1.5 at 0.5 Bmsy, -> 4.3 at 0. "power": see proc_dep_power. Capped at proc_sd_max.
proc_dep_pars: parameters a, c, d of the JABBA proc.dyn function.
proc_dep_power, proc_dep_ref, proc_sd_max: power form (proc_dep = "power"): sd_i,t = sd_i * ((B_i,t / K_i) / proc_dep_ref)^(-proc_dep_power), capped at proc_sd_max; sd_i applies at B/K = proc_dep_ref (Bmsy). proc_dep_power > 0: stocks become more volatile at low abundance; < 0: less volatile
proc_center: "mean": exp(sd eps - sd^2/2) (mean-unbiased, the median drifts down as sd grows: stochastic depensation); "median": exp(sd eps), as in JABBA's dlnorm(Pmean, iPV) process equation
impl_cv, impl_bias: implementation error on the effort caps by segment: realised cap = cap * exp(bias + N(-s^2/2, s^2)), s from impl_cv
effort_cv, effort_rho: effort/catchability noise by segment, in every run (regulated or not): an annual lognormal multiplier on effective effort (weather, fuel, markets, catchability), AR(1) effort_rho, mean 1. Unlike the implementation error it acts when no cap binds.
seed: random seed.
crn: common random numbers: rows i and i + crn share the same process (and rule) error draws, e.g. crn = nsim when the rows are nsim simulations repeated for each of several control settings
keep: "final" (equilibrium summaries) or "all" (trajectories)
avg_years: years kept when keep = "final".
F_cap: cap on each stock's realised F, in units of Fmsy (scalar or one per stock; default sys$F_cap, else Inf). F/Fmsy is measured as the harvest-rate ratio h/hmsy (see ffmsy()); the cap on h is converted to instantaneous F and applied monthly as one twelfth. Extra effort on a capped stock catches nothing more, so its catch per day falls (economic extinction: the stock is left to rare by-catch). In Schaefer equilibrium B/Bmsy = 2 - F/Fmsy, so a cap below 2 sets a floor: F_cap = 1.95 keeps every stock at or above 0.05 Bmsy.
F_cap_p: shape of the cap: Inf = hard cap, min(F, F_cap); a finite p = smooth cap, (F^-p + F_cap^-p)^(-1/p), which approaches F_cap without a kink (p = 8 lowers F by about 1.5% at 1.5 Fmsy when F_cap = 1.95). Default NULL: hard without process error (equilibrium and tuning), p = 8 with process error (stochastic runs).

Value

An object of class ms_sim: arrays B, C, F (years x rows x stocks), E_metier, E_seg_metier, E_seg, matrices revenue, profit, the system, and init (the year-0 biomass and days, rows x stocks / segments).