deparse1 <- function (expr, collapse = " ", width.cutoff = 500L, ...) paste(deparse(expr, width.cutoff, ...), collapse = collapse) ## needed on stockassessment.org library(stockassessment) load("run/model.RData") thisY <- as.numeric(substring(Sys.Date(),1,4)) # intermediate year lastY <- thisY-1 # last year of assessment ## General settings: seed <- 12345 ny <- 5 # number of years in forecast (including last y of assessment) nsim <- 1000 ave.years <- max(fit$data$years) + (-4:0) # years for weights year.base <- max(fit$data$years) # years for selectivity ## Reference points: Fpa <- 0.050 Fmsy <- Fpa Blim <- 197344 Bpa <- 274225 MSYBtrig <- Bpa ## Intermediate year assumption: TAC_F <- 788 TAC_NS3a <- 328566 EUextra <- 560 # from https://eur-lex.europa.eu/legal-content/EN/TXT/PDF/?uri=OJ:L_202600249 load("data/TAC_utilization.RData") catch_WBSS_43a <- sum(table[as.character(lastY),"WBSS_IVaE"], table[as.character(lastY),"WBSS_C.fleet"], table[as.character(lastY),"WBSS_D.fleet"]) catch_NSAS_43a <- sum(table[as.character(lastY),"A.fleet.NSAS."], table[as.character(lastY),"B.Fleet"], table[as.character(lastY),"NSAS_C.fleet"], table[as.character(lastY),"NSAS_D.fleet"]) splitABCD <- catch_WBSS_43a/sum(catch_WBSS_43a, catch_NSAS_43a) TAC_all <- splitABCD*(TAC_NS3a+EUextra)+TAC_F FC<-list() ## F=0 ---- set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, #constraints = c(paste0("C=", TAC_all), rep(paste0("F=", 1e-6), ny-2)), catchval=c(TAC_all,rep(NA,ny-2)), fval=c(NA,rep(1e-6, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("F=0"), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) ## ICES MSY approach: ICES AR with precautionary considerations, i.e., F=0 below Blim ---- SSB_intermY <- attributes(FC[[1]])$tab[as.character(thisY),"ssb:median"] if(SSB_intermY=MSYBtrig) Fm <- Fmsy else Fm <- Fmsy*SSB_intermY/MSYBtrig set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, #constraints = c(paste0("C=", TAC_all), rep(paste0("F=", Fm), ny-2)), catchval=c(TAC_all,rep(NA,ny-2)), fval=c(NA,rep(Fm, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("ICES MSY approach: AR with F=0 below Blim"), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) # what about in year adv+1? should it be the ICES AR? if Yes need to optimize for future years! ## ICES EU Baltic Sea multiannual plan (MAP) ---- SSB_intermY <- attributes(FC[[1]])$tab[as.character(thisY),"ssb:median"] if(SSB_intermY>=MSYBtrig) Fm <- Fmsy else Fm <- Fmsy*SSB_intermY/MSYBtrig set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, #constraints = c(paste0("C=", TAC_all), rep(paste0("HCR=", Fmsy, "~", MSYBtrig, "SSB-1"), ny-2)), catchval=c(TAC_all,rep(NA,ny-2)), fval=c(NA,rep(Fm, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("EU Baltic Sea multiannual plan (MAP): F=Fmsy*SSB", thisY, "/MSYBtrigger"), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) # ## ICES EU Baltic Sea multiannual plan (MAP) with Flower ---- # SSB_intermY <- attributes(FC[[1]])$tab[as.character(thisY),"ssb:median"] # if(SSB_intermY>=MSYBtrig) Fm <- Fmsylower else Fm <- Fmsylower*SSB_intermY/MSYBtrig # set.seed(seed) # FC[[length(FC)+1]] <- modelforecast(fit, # #constraints = c(paste0("C=", TAC_all), rep(paste0("HCR=", Fmsy, "~", MSYBtrig, "SSB-1"), ny-2)), # catchval=c(TAC_all,rep(NA,ny-2)), # fval=c(NA,rep(Fm, ny-2)), # nosim=nsim, # ave.years = ave.years, # year.base = year.base, # label=paste0("MAP: F=Fmsylower*SSB", thisY, "/MSYBtrigger"), # resampleFirst = TRUE, # deterministicF = FALSE, # processNoiseF = FALSE) ## F=Fmsy=Fpa ---- set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, constraints = c(paste0("C=", TAC_all), rep(paste0("F=", Fmsy), ny-2)), #catchval=c(TAC_all,rep(NA,ny-2)), #fval=c(NA,rep(Fmsy, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("F=FMSY=Fpa=",Fmsy), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) ## SSB(adv+1)=Blim ---- ## SSB(adv+1)=Bpa=MSYBtrigger ---- ## F=Fsq i.e., intermediate year F ---- FthisY <- attributes(FC[[1]])$tab[which(rownames(attributes(FC[[1]])$tab)==thisY),1] # F in intermediate year set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, constraints = c(paste0("C=", TAC_all), rep(paste0("F=", FthisY), ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("F=F", thisY, "=",FthisY), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) ## F=0.005 ---- Flow1 <- 0.005 set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, constraints = c(paste0("C=", TAC_all), rep(paste0("F=", Flow1), ny-2)), #catchval=c(TAC_all,rep(NA,ny-2)), #fval=c(NA,rep(Fmsy, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("F=",Flow1), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) ## F=0.01 ---- Flow2 <- 0.01 set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, constraints = c(paste0("C=", TAC_all), rep(paste0("F=", Flow2), ny-2)), #catchval=c(TAC_all,rep(NA,ny-2)), #fval=c(NA,rep(Fmsy, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("F=",Flow2), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) ## F=0.02 ---- Flow3 <- 0.02 set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, constraints = c(paste0("C=", TAC_all), rep(paste0("F=", Flow3), ny-2)), #catchval=c(TAC_all,rep(NA,ny-2)), #fval=c(NA,rep(Fmsy, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("F=",Flow3), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) ## F=0.03 ---- Flow4 <- 0.03 set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, constraints = c(paste0("C=", TAC_all), rep(paste0("F=", Flow4), ny-2)), #catchval=c(TAC_all,rep(NA,ny-2)), #fval=c(NA,rep(Fmsy, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("F=",Flow4), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) ## C=interimC ---- set.seed(seed) FC[[length(FC)+1]] <- modelforecast(fit, constraints = c(rep(paste0("C=", TAC_all), ny-1)), #catchval=c(TAC_all,rep(NA,ny-2)), #fval=c(NA,rep(Fmsy, ny-2)), nosim=nsim, ave.years = ave.years, year.base = year.base, label=paste0("Constant intermediate year catch = ",round(TAC_all), " t"), resampleFirst = TRUE, deterministicF = FALSE, processNoiseF = FALSE) FC <- FC[-1] # remove F=0 because same as scenario 2 extract_prob <- function(FC, Btarget, year){ SSBforecast <- sapply(FC,function(x) (x$ssb)) colnames(SSBforecast) <- rownames(attr(FC, "tab")) probs <- apply(SSBforecast, 2, function(x) mean(x