Simulate Historical Fishery

We’ll use the Hogfish_PR_NOAA OM from the OM library to demonstrate the simulation of the fishery.

# download MSEextra from GitHub
MSEextra()

library(MSEextra)
OM <- MSEextra::Hogfish_PR_NOAA

For demonstration purposes, we will set nsim to 3:

OM@nsim <- 3

The Simulate function simulates the historical fishery and returns an object of class Hist:

Hist <- Simulate(OM, silent=TRUE)

We will step through the key calculations that take place inside the Simulate function.

Ensuring Results are Reproducible

To ensure the simulations are consistent across machines, OM@seed is used to set the seed for the random number generator:

set.seed(OM@seed)

Sampling Custom Parameters

Next nsim values of the OM parameters are sampled from the operating model.

First, the Custom Parameters slot (OM@cpars) is sampled:

SampCpars <- SampleCpars(cpars=OM@cpars, OM@nsim)
## No valid names found in custompars. Ignoring OM@cpars

OM@cpars is not used in this operating model, so SampCpars is a empty list. If OM@cpars was populated with Valid Custom Parameters nsim samples of each custom parameter would be drawn from OM@cpars.

Correlations are maintained when sampling from OM@cpars.

Sampling OM Objects

Next, the Stock, Fleet, Obs, and Imp parameters are sampled from the OM:

StockPars <- SampleStockPars(OM, cpars=SampCpars)
## Note: Maximum age ( 23 ) is lower than assuming 1% of cohort survives
## to maximum age ( 32 )
## Optimizing for user-specified movement
FleetPars <- SampleFleetPars(OM, StockPars, cpars=SampCpars)
ObsPars <- SampleObsPars(OM, cpars=SampCpars, Stock=StockPars)
ImpPars <- SampleImpPars(OM, cpars=SampCpars)

A few things to notice here:

• There is a message alerting us that OM@maxage is lower than the age where about 1% of a cohort would survive to. This is usually not a problem because the model groups all age-classes beyond OM@maxage into a plus-group. But given there is little computational cost in increasing maximum age, users may wish to set OM@maxage <- 32 so that more age-classes are explicitly modeled.

• SampCpars is provided to each of the functions that sample the Stock, Fleet, Obs, and Imp parameters. Parameters that are included in the Custom Parameters are used and these parameters are not sampled from the OM object.

• StockPars contains the biological parameters for the spool-up and projection periods, e.g., StockPars$M_ageArray is an array with M-at-age with dimensions c(OM@nsim, OM@maxage+1, OM@nyears+OM@proyears). • FleetPars contains the exploitation parameters for the spool-up period (e.g, FleetPars$Find) and the selectivity/retention parameters for the spool-up and projection period, e.g., FleetPars$V is an array with vulnerability-at-age with dimensions c(OM@nsim, OM@maxage+1, OM@nyears+OM@proyears). The selectivity and retention parameters may be updated in the projection period if an MP changes the selectivity/retention pattern. • ObsPars contains the observation parameters for all years (e.g., ObsPars$Cobs_y) is a matrix with dimensions c(OM@nsim, OM@nyears+OM@proyears) with catch observation error.

• ImpPars contains the implementation error for the projection years; e.g., ImpPars$TAC_y is a matrix with dimensions c(OM@nsim, OM@proyears) with TAC implementation error. Calculate Unfished Reference Points The unfished reference points are calculated by running the simulation model with no fishing mortality and no recruitment deviations. The equilibrium unfished reference points are calculated for each year. For example, the annual unfished spawning stock biomass (SSB0) is calculated for the 75 historical years and the 40 projection years: matplot(t(Hist@Ref$ByYear$SSB0), type="l", xlab="Year", ylab="SSB0", lwd=2) The small amount of inter-annual variability in SSB0 is caused by inter-annual variability in the life-history parameters (e.g., OM@Msd). The inter-annual variability in SSB0 would be greater if there was more variability in OM@Msd, OM@Linfsd, and OM@Ksd. The SSB0 used to determine depletion is calculated by averaging Hist@Ref$ByYear$SSB0 over several years. See Depletion and Unfished Reference Points for details. Dynamic unfished reference points (i.e., including recruitment deviations) are also calculated and returned in Hist@Ref$Dynamic_Unfished. For example, the dynamic unfished spawning stock biomass is:

matplot(t(Hist@Ref$Dynamic_Unfished$SSB0), type="l", xlab="Year",
ylab="Dynamic SSB0", lwd=2)

Calculate MSY Reference Points

Maximum sustainable yield (MSY) reference points are calculated using the information in StockPars and FleetPars. MSY reference points are calculated in each year, assuming equilbrium conditions (no recruitment deviations) and the life-history and selectivity parameters in each year.

For example,

matplot(t(Hist@Ref$ByYear$MSY), type="l", xlab="Year", ylab="MSY", lwd=2)

See MSY Reference Points for more details.

Optimize for Depletion

openMSE was designed to allow easy access to MSE modeling for data-limited fisheries. By definition, historical fishing mortality and current depletion is not known in a data-limited fishery.

Instead, the general pattern in fishing mortality is specified in the Fleet Object and a range of depletion values in the Stock Object (or in Custom Parameters if distributions other than uniform are required).

For each simulation, the model uses the relevant parameters (relative pattern in fishing mortality, size-at-age in each year, selectivity, recruitment deviations, etc) and optimizes for the catchability coefficient q so that the depletion at the end of the spool-up period matches the value drawn from OM@D.

We demonstrate with an example for a single simulation:

sim <- 1
# Depletion value for this simulation
Hist@SampPars$Stock$Depletion[sim]
## [1] 0.06173859
# SSB time-series (sum over areas)