##### FORECAST SCRIPT TO COPY TO STOCKASSESSMENT.ORG BELOW ##### ## Forecasts for catch option to table 3 # The following catch options are explored below: # thisY (intermediate year) # MSY approach (F=0 below Blim) # MAP 2018: F=Fmsy x SSBy-1/MSY_trig # MAP 2018: F=Fmsy.lower x SSBy-1/MSY_trig # MAP 2018: F=Fmsy.upper x SSBy-1/MSY_trig # F=Fmsy # F=Fpa # F=Flim # NOT APPLIED F= F{SSB (thisY+1) = Blim} # NOT APPLIED F= F{SSB (thisY+1) = Bpa} # NOT APPLIED F= F{SSB (thisY+1) = MSY.Btrigger} # F= F(thisY) # F=0 # F=0.05 # F=0.10 # F=0.15 # Constant thisY TAC # library(stockassessment) source("src/forecast2.R") # forecast function of Vanessa with HCR in forecast (used for scenario when Fy depends on SSBy-1) source("src/myforecast.R") # change in original forecast function to avoid optimization problem when catch=0 for few fleets load("run/model.RData") FC<-list() # Forecast make empty list) thisY <- as.numeric(substring(Sys.Date(),1,4)) # intermediate year lastY <- thisY-1 # last year of assessment TACa_NSWB <- 388542 # 510323 in 2024 # 396556 in 2023 # 427628 in 2022 # 356357 in 2021 TACc_NSWB <- 22793 # 29735 in 2024 # 23250 in 2023 # 25021 in 2022 # 21604 in 2021 TACd_NSWB <- 6659 # same as 2024 TACf <- 788 # same as 2022-2024 # 1575 in 2021 splitC <- 0.303283963698 # 0.440036067093 in 2024 # 0.569706976430 in 2023 # 0.644895257674001 in 2022 # last 3 year split !!!! As many decimals as possible! splitD <- 0.142040973185 # 0.215400861011 in 2024 # 0.302424213474 in 2023 # 0.30000367652937 in 2022 # last 3 year split !!!! As many decimals as possible! utiliz_A <- 1 # 100% utilization # utiliz_C <- 0.50 # 0.50 in 2020 given by DK and NOR fishers # New rule brexit # utiliz_D <- 0.030999366525 # 0.065245347342 in 2023 # 0.0706991180910504 in 2022 # 3 year average utilization of TAC, given with split utiliz_F <- 1 # 100% utilization # ny_mean <- 3 # last number of years for average TAC in intermediate year # TACa_WB <- mean(catchbyfleettable(fit, obs.show = TRUE)[(fit$data$noYears-ny_mean+1):fit$data$noYears,13]) # average of last ny_mean observed catch # changes in 2022 for A fleet for same assumption with NSAS NOR_3a <- 200 # 0.1*3102 in 2023 EU_3a <- 969 # same as 2023 prop_EU_C_lastY <- 0.774873040413 # 0.424504914 in 2024 # 0.474510878 in 2023 prop_EU_D_lastY <- 0.225126959587 # 0.575495086 in 2024 # 0.525489122 in 2023 transferC <- TACc_NSWB-(prop_EU_C_lastY*EU_3a+NOR_3a) # TACc_NSWB-(EU_3a+NOR_3a) in 2023 splitA <- 0.008570511880 # 0.009084268546 in 2024 # 0.012620318649 in 2023 # prop of WBSS in total A fleet catch TACa_WB <- (transferC+TACa_NSWB)*splitA TACc_WB <- splitC*TACc_NSWB TACd_WB <- splitD*TACd_NSWB TAC_A <- utiliz_A*TACa_WB #TAC_C <- utiliz_C*TACc_WB # in 2020 #TAC_C <- TACc_WB-((0.5*2881+3000)*splitC) # new 2021 EU/UK/NO consultations (50% transfer NOR TAC (50% of 2881), and 3000 t for EU) #TAC_D <- utiliz_D*TACd_WB TAC_D <- (prop_EU_D_lastY*EU_3a)*splitD # new 2023 TAC_C <- (prop_EU_C_lastY*EU_3a+NOR_3a)*splitC# (prop_EU_C_lastY*EU_3a+0.1*NOR_3a)*splitC # new 2022 TAC_F <- utiliz_F*TACf TAC_all <- sum(TAC_A, TAC_C, TAC_D, TAC_F) Ry <- (lastY-5):(lastY-1) # years to sample for recruitment (last 5 years except lastY) Fmsy <- 0.31 Flower <- 0.216 Fupper <- 0.379 #Fpa <- 0.35 Fpa <- 0.41 #Fp05 now Flim <- 0.45 Blim=120000 MSYBtrig=150000 Flow=0.1 nsim <- 2 seed <- 12345 estimate <- median determ <- TRUE nf <- sum(fit$data$fleetTypes==0) # number of catching fleets ny <- 5 # number of years in forecast (including last y of assessment) ### From 2022: remove F constraint of 50% catch set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) #cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,NA,NA,NA,NA) #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(1e-10, ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, rec.years=Ry, label="MSY approach (zero catch)", nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below #************************ Here we use forecast2 function created by Vanessa to use HCR in the forecast *************************** FC[[length(FC)+1]] <- forecast2(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), #fval=c(NA,NA,0,1180562, rep(NA,ny-3)), Fscenario=2, # hockey-stick shape for HCR MSYBtrig=c(NA, NA, rep(MSYBtrig,ny-2)), Blim=c(NA, NA, rep(Blim,ny-2)), Fmsy=Fmsy, Flow=Flow, cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label=paste0("MAP 2018: F=FMSY(", Fmsy, ")*SSBy-1/MSYBtrigger"), nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below #************************ Here we use forecast function2 created by Vanessa to use HCR in the forecast *************************** FC[[length(FC)+1]] <- forecast2(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), #fval=c(NA,NA,0.08225852, rep(NA,ny-3)), Fscenario=2, # hockey-stick shape for HCR MSYBtrig=c(NA, NA, rep(MSYBtrig,ny-2)), Blim=c(NA, NA, rep(Blim,ny-2)), Fmsy=Flower, Flow=Flow, cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label=paste0("MAP 2018: F=FMSYlower(", Flower, ")*SSBy-1/MSYBtrigger"), nosim=nsim,estimate=estimate, deterministic=determ) ## Removed in advice 2022 # set.seed(seed) # same seed to repeat output # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch # #catch fraction catch weight keeping total catch w # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we use forecast function2 created by Vanessa to use HCR in the forecast *************************** # FC[[length(FC)+1]] <- forecast2(fit, fscale=c(1,rep(NA,ny-1)), # catchval=c(NA,TAC_all,rep(NA,ny-2)), # #fval=c(NA,NA,0.1443332, rep(NA,ny-3)), # Fscenario=2, # hockey-stick shape for HCR # MSYBtrig=c(NA, NA, rep(MSYBtrig,ny-2)), # Blim=c(NA, NA, rep(Blim,ny-2)), # Fmsy=Fupper, # Flow=Flow, # cf.cv.keep.cv=cf.cv.keep.cv, # #cf.cv.keep.fv=cf.cv.keep.fv, # rec.years=Ry, # label=paste0("MAP 2018: F=FMSYupper(", Fupper, ")*SSBy-1/MSYBtrigger"), # nosim=nsim,estimate=estimate, deterministic=determ) # set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(Fmsy,ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label=paste0("F=FMSY=",Fmsy), nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(Fpa,ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label=paste0("F=Fpa=", Fpa), nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(Flim,ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label=paste0("F=Flim=", Flim), nosim=nsim,estimate=estimate, deterministic=determ) #******** Here 3 options are not projected ********************* # NOT APPLIED F= F{SSB (2021) = Blim} # NOT APPLIED F= F{SSB (2021) = Bpa} # NOT APPLIED F= F{SSB (2021) = MSY.Btrigger} #******** continue with table ********************************** FthisY <- attributes(FC[[1]])$tab[which(rownames(attributes(FC[[1]])$tab)==thisY),1] # F in intermediate year set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(FthisY,ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label=paste0("F=F", thisY, "=", FthisY), nosim=nsim,estimate=estimate,deterministic=determ) # set.seed(seed) # same seed to repeat output # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,NA,NA,NA,NA) # #catch fraction catch weight keeping total catch w # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we iterate manually year by year with a new SSB each year *************************** # FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,NA,NA,NA,NA), # catchval=c(NA,TAC_all,rep(NA,ny-2)), # fval=c(NA,NA,0.203, 0.217, 0.238), # cf.cv.keep.cv=cf.cv.keep.cv, # rec.years=Ry, # label="MAP F-MSY 0.32*B/110kt", # nosim=nsim,estimate=estimate, deterministic=determ) # # # # set.seed(seed) # same seed to repeat output # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,NA,NA,NA,NA) # #catch fraction catch weight keeping total catch w # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we iterate manually year by year with a new SSB each year *************************** # FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,NA,NA,NA,NA), # catchval=c(NA,TAC_all,rep(NA,ny-2)), # fval=c(NA,NA,0.146, 0.157, 0.181), # cf.cv.keep.cv=cf.cv.keep.cv, # rec.years=Ry, # label="MAP F-Low 0.23*B/110kt", # nosim=nsim,estimate=estimate, deterministic=determ) # # # # set.seed(seed) # same seed to repeat output # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,NA,NA,NA,NA) # #catch fraction catch weight keeping total catch w # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we iterate manually year by year with a new SSB each year *************************** # FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,NA,NA,NA,NA), # catchval=c(NA,TAC_all,rep(NA,ny-2)), # fval=c(NA,NA,0.260, 0.277, 0.287), # cf.cv.keep.cv=cf.cv.keep.cv, # rec.years=Ry, # label="MAP F-High 0.41*B/110kt", # nosim=nsim,estimate=estimate, deterministic=determ) # # Extra run asked by ACOM in 2020: TACextra = sum(TAC_A,TAC_D) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) for (i in 3:nrow(cf.cv.keep.cv)) cf.cv.keep.cv[i,] <- c(NA,NA,NA,NA,NA,1E-10,TAC_D,1E-10) #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- myforecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all, rep(TACextra,ny-2)), fval=c(NA,NA,rep(NA,ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, rec.years=Ry, label=paste0("Catch for bycatch fleets only"), nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(1e-10, ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, rec.years=Ry, label="F=0", nosim=nsim,estimate=estimate, deterministic=determ) # Extra run 2022: set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(0.01, ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label="F=0.01", nosim=nsim,estimate=estimate, deterministic=determ) # Extra run 2022: set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(0.025, ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label="F=0.025", nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(0.05, ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label="F=0.05", nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(0.1, ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label="F=0.1", nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # cf.cv.keep.fv<-matrix(NA,ncol=nf*2,nrow=ny) # for (t in 3:nrow(cf.cv.keep.fv)) cf.cv.keep.fv[t,] <- c(NA,NA,NA,0.5,NA,NA,NA,NA) # so F takes 50% of the catch #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,TAC_all,rep(NA,ny-2)), fval=c(NA,NA,rep(0.15, ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, #cf.cv.keep.fv=cf.cv.keep.fv, rec.years=Ry, label="F=0.15", nosim=nsim,estimate=estimate, deterministic=determ) set.seed(seed) # same seed to repeat output cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=ny) cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) for (i in 3:nrow(cf.cv.keep.cv)) cf.cv.keep.cv[i,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) #catch fraction catch weight keeping total catch w #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # only specify 3 because total cv is given in "catchval" below FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,rep(NA,ny-1)), catchval=c(NA,rep(TAC_all,ny-1)), fval=c(NA,NA,rep(NA,ny-2)), cf.cv.keep.cv=cf.cv.keep.cv, rec.years=Ry, label=paste0("Constant ", thisY, " TAC"), nosim=nsim,estimate=estimate, deterministic=determ) # set.seed(seed) # same seed to repeat output # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=6) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # #Acatch <- catchbyfleettable(fit, obs.show=TRUE)[,"Obs(Fleet w.o. effort 1)"] # #Atac <- mean(Acatch[length(Acatch)-2:0]) # #Atac <- 291042*0.5/100 # Atac <- 1558 # #Dcatch <- catchbyfleettable(fit, obs.show=TRUE)[,"Obs(Fleet w.o. effort 3)"] # #DNScatch<-c(4.4003, 1.4525, 0.1989)*1000 # #Dtac <- mean(Dcatch[length(Dcatch)-2:0]/(Dcatch[length(Dcatch)-2:0]+DNScatch))*6659 # Dtac <- 0.4*6659 # Ftac <- TAC_F # cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[4,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[5,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[6,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # #catch fraction catch weight keeping total catch wac # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we iterate manually year by year with a new SSB each year *************************** # # # thisfc <- forecast(fit, fscale=c(1,NA,NA,NA,NA,NA), # catchval=c(NA,TAC_all,2*Ftac,2*Ftac,2*Ftac,2*Ftac), # fval=c(NA,NA,NA,NA,NA,NA), # cf.cv.keep.cv=cf.cv.keep.cv, # rec.years=Ry, # label="Henrik option 29Oct2018", # nosim=nsim,estimate=estimate, deterministic=determ) # catch_2019<-attr(thisfc, "catchby")["2019",c(1,4,7,10)] # FC[[length(FC)+1]] <- thisfc # # # # # # set.seed(seed) # same seed to repeat output # newsum<-catch_2019[1]+.54*catch_2019[2]+.46*catch_2019[3]+catch_2019[4] # #newsum<-1445+.54*4882+.46*2664+TAC_F # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=6) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # #Acatch <- catchbyfleettable(fit, obs.show=TRUE)[,"Obs(Fleet w.o. effort 1)"] # #Atac <- mean(Acatch[length(Acatch)-2:0]) # #Atac <- 291042*0.5/100 # Atac <- 1558 # #Dcatch <- catchbyfleettable(fit, obs.show=TRUE)[,"Obs(Fleet w.o. effort 3)"] # #DNScatch<-c(4.4003, 1.4525, 0.1989)*1000 # Dtac <- .46*0.4*6659 #mean(Dcatch[length(Dcatch)-2:0]/(Dcatch[length(Dcatch)-2:0]+DNScatch))*6659* # Ftac <- TAC_F # cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[4,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[5,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[6,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # #catch fraction catch weight keeping total catch w # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we iterate manually year by year with a new SSB each year *************************** # FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,NA,NA,NA,NA,NA), # catchval=c(NA,TAC_all,newsum,newsum,newsum,newsum), # fval=c(NA,NA,NA,NA,NA,NA), # cf.cv.keep.cv=cf.cv.keep.cv, # rec.years=Ry, # label="Henrik option02 29Oct2018", # nosim=nsim,estimate=estimate, deterministic=determ) # # # # # # set.seed(seed) # same seed to repeat output # #newsum<-1445+.54*4882+.46*2664+TAC_F # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=6) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # #Acatch <- catchbyfleettable(fit, obs.show=TRUE)[,"Obs(Fleet w.o. effort 1)"] # #Atac <- mean(Acatch[length(Acatch)-2:0]) # #Atac <- 291042*0.5/100 # Atac <- 1558 # #Dcatch <- catchbyfleettable(fit, obs.show=TRUE)[,"Obs(Fleet w.o. effort 3)"] # #DNScatch<-c(4.4003, 1.4525, 0.1989)*1000 # Dtac <- 0.4*6659 #mean(Dcatch[length(Dcatch)-2:0]/(Dcatch[length(Dcatch)-2:0]+DNScatch))*6659* # Ftac1 <- TAC_F # p=0.9637218 # found after minimization # Ftac2 <- p*TAC_F # cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,Atac,NA,0.46*Dtac,Ftac1) # cf.cv.keep.cv[4,] <- c(NA,NA,NA,NA,Atac,NA,0.46*Dtac,Ftac2) # cf.cv.keep.cv[5,] <- c(NA,NA,NA,NA,Atac,NA,0.46*Dtac,Ftac2) # cf.cv.keep.cv[6,] <- c(NA,NA,NA,NA,Atac,NA,0.46*Dtac,Ftac2) # Ctac1=(Ftac1-Atac-Dtac)*.54 # Ctac2=(Ftac2-Atac-Dtac)*.54 # Tot1=Atac+0.46*Dtac+Ctac1+Ftac1 # Tot2=Atac+0.46*Dtac+Ctac2+Ftac2 # #catch fraction catch weight keeping total catch w # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we iterate manually year by year with a new SSB each year *************************** # FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,NA,NA,NA,NA,NA), # catchval=c(NA,TAC_all,Tot1,Tot2,Tot2,Tot2), # fval=c(NA,NA,NA,NA,NA,NA), # cf.cv.keep.cv=cf.cv.keep.cv, # rec.years=Ry, # label="Henrik option03 29Oct2018", # nosim=nsim,estimate=estimate, deterministic=determ) # # # # set.seed(seed) # same seed to repeat output # #newsum<-1445+.54*4882+.46*2664+TAC_F # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=6) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # Atac <- 1437 # Dtac <- 0.4*6659 #mean(Dcatch[length(Dcatch)-2:0]/(Dcatch[length(Dcatch)-2:0]+DNScatch))*6659* # Ftac <- TAC_F # Ctac <- 16637 # cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[4,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[5,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[6,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # Tot <- Atac+Ctac+Dtac+Ftac # #catch fraction catch weight keeping total catch w # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we iterate manually year by year with a new SSB each year *************************** # FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,NA,NA,NA,NA,NA), # catchval=c(NA,TAC_all,Tot,Tot,Tot,Tot), # fval=c(NA,NA,NA,NA,NA,NA), # cf.cv.keep.cv=cf.cv.keep.cv, # rec.years=Ry, # label="Henrik option04 22Nov2018", # nosim=nsim,estimate=estimate, deterministic=determ) # # # set.seed(seed) # same seed to repeat output # #newsum<-1445+.54*4882+.46*2664+TAC_F # cf.cv.keep.cv<-matrix(NA,ncol=nf*2,nrow=6) # cf.cv.keep.cv[2,] <- c(NA,NA,NA,NA,NA,TAC_C,TAC_D,TAC_F) # Atac <- 1437 # Dtac <- 0.46*0.4*6659 #mean(Dcatch[length(Dcatch)-2:0]/(Dcatch[length(Dcatch)-2:0]+DNScatch))*6659* # Ftac <- TAC_F # Ctac <- 0.54*16637 # cf.cv.keep.cv[3,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[4,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[5,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # cf.cv.keep.cv[6,] <- c(NA,NA,NA,NA,Atac,NA,Dtac,Ftac) # Tot <- Atac+Ctac+Dtac+Ftac # #catch fraction catch weight keeping total catch w # #There are ncol 4 fleets * 2 ways to define catches 1=fraction; 2=tons # # only specify 3 because total cv is given in "catchval" below # #************************ Here we iterate manually year by year with a new SSB each year *************************** # FC[[length(FC)+1]] <- forecast(fit, fscale=c(1,NA,NA,NA,NA,NA), # catchval=c(NA,TAC_all,Tot,Tot,Tot,Tot), # fval=c(NA,NA,NA,NA,NA,NA), # cf.cv.keep.cv=cf.cv.keep.cv, # rec.years=Ry, # label="Henrik option05 22Nov2018", # nosim=nsim,estimate=estimate, deterministic=determ) # save(FC, file="run/forecast.RData")