spict
spict copied to clipboard
Problem with catch scenarios for longer management periods
There are sometimes issues when the management period is longer than one year and Catch fractiles are used. The F/Fmsy at the end of the management period is higher than 1.
It seems to affect assessments with low biomass in the end of the time series.
A possible solution is to split Cp and dtcp to yearly chunks instead of a single long catch duration.
ReprEx
library(spict)
#> Loading required package: TMB
#> Welcome to spict_v1.3.5@8a4e1b
set.seed(1234)
nt <- 50
inp <- list(nseasons=4, splineorder=3)
s1 <- seq(0, 10, by=1)
s2 <- seq(11, nt-3-1/inp$nseasons, by=1/inp$nseasons)
s3 <- seq(nt-3, nt-1, by = 1)
inp$timeC <- c(s1,s2,s3)
inp$dtc <- c(rep(1, length(s1)),
rep(0.25, length(s2)),
rep(1, length(s3)))
inp$timeI <- seq(0.1, nt-1/inp$nseasons, by=0.5)
inp$ini <- list(logK=log(1000), logm=log(123), logq=log(1), logn=log(2),
logbkfrac=log(0.7), logsdf=log(0.123), logF0=log(0.4321),
logphi=log(c(0.3, 0.5, 1.8)))
inp$dteuler <- 1/8
inpsim <- sim.spict(inp)
#> The stabilise option was used for fitting / is specified in the input list (inp$stabilise). This activates six very vague priors. Acknolowdging these priors when simulating is not yet implemented.
#> Additional priors were used for fitting / are specified in the input list (inp$priors): logn, logalpha, logbeta. Acknolowdging these priors when simulating is not yet implemented.
inpsim$dtc <- inp$dtc ## Very important, sim.spict drops it
inpsim$maneval <- nt+5
inpsim$maninterval <- c(nt,nt+5)
inpsim$priors$logbkfrac <- c(log(0.7), 0.001, 1)
inpsim$priors$logn <- c(log(2), 0.001, 1)
inpsim$dtc <- inp$dtc
inpsim <- check.inp(inpsim)
plotspict.data(inpsim)
#> Error in xy.coords(x, y, xlabel, ylabel, log): 'x' and 'y' lengths differ
fit <- fit.spict(inpsim)
fit <- add.man.scenario(fit, "Fmsy")
fit <- add.man.scenario(fit, "C35", fractiles = list(catch = 0.35))
fit <- add.man.scenario(fit, "All35", fractiles = list(catch = 0.35, ffmsy=0.35, bbmsy=0.35))
fit <- add.man.scenario(fit, "FB35_noC", fractiles = list( ffmsy=0.35, bbmsy=0.35))
plot2(fit) ## This does not, indxmax is not defined

sumspict.manage(fit)
#> SPiCT timeline:
#>
#> Observations Management
#> 0.00 - 50.00 50.00 - 55.00
#> |-----------------------| ----------------------|
#>
#> Management evaluation: 55.00
#>
#> Predicted catch for management period and states at management evaluation time:
#>
#> C B/Bmsy F/Fmsy
#> 1. Fish at Fmsy 166.1 0.27 1.00
#> 2. C35 148.1 0.02 2.53
#> 3. All35 145.3 0.02 2.66
#> 4. FB35_noC 163.3 0.30 0.92
Created on 2022-03-08 by the reprex package (v2.0.1)