spict icon indicating copy to clipboard operation
spict copied to clipboard

Problem with catch scenarios for longer management periods

Open alko989 opened this issue 3 years ago • 0 comments

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)

alko989 avatar Mar 08 '22 14:03 alko989