sprat-xsa-experimental

1
library("msy")
2

3
# Blim = 0.2 * SSBmax
4
result.brp.blim <- max(ssb(result.xsa.stock)) * 0.2
5
# get fbar values as vector
6
fvals <- as.vector(fbar(result.xsa.stock))
7

8
# found extremes of whiskers by boxplot status to remove unreal data
9
wstats <- boxplot.stats(fvals)$stats
10
# remove values out of box of extremums 
11
fvals <- fvals[fvals > wstats[[1]]]
12
fvals <- fvals[fvals < wstats[[5]]]
13

14
# get variance (standard deviation) of Fbar
15
fsd <- sd(fvals)
16

17
# Bpa = Blim * e^(1.645*stdev)
18
result.brp.bpa <- result.brp.blim * exp(qnorm(0.95)*fsd)
19

20

21
# fit SR models via eqsim
22
result.eqsim.srmodel <- eqsr_fit(result.xsa.stock, 
23
                                 n = 5000, 
24
                                 models = c("Ricker", "Bevholt", "Segreg")
25
                                 #remove.years = c(2013,2014,2015)
26
                                 )
27

28
lastYear <- as.numeric(result.xsa.stock@range["maxyear"])
29

30
# calculate reference point based on model
31
result.eqsim.brp <- eqsim_run(result.eqsim.srmodel, 
32
                              bio.years = c(lastYear-3, lastYear), 
33
                              extreme.trim = c(0.05, 0.95), 
34
                              Bpa = result.brp.bpa,
35
                              Blim = result.brp.blim)
36

37
# create new FLR sr-model, segreg from eqsim params
38
tmp.model <- filter(result.eqsim.srmodel$sr.det, model=="Segreg")
39
result.sr.segreg <- FLSR(model = "segreg")
40
result.sr.segreg@params["a"] <- as.numeric(tmp.model["a"])
41
result.sr.segreg@params["b"] <- as.numeric(tmp.model["b"])
42

43
# get mean recruitment of latest 5 years without terminal year
44
tmp.rec.mean <- exp(mean(log(window(rec(result.xsa.stock), 2014, 2018))))
45
# build geomean SR model
46
result.sr.geomean <- FLSR(model = "geomean", params=FLPar(tmp.rec.mean))
47

48
# build test model
49
test.geomean <- fwd(data.stf, sr=result.sr.geomean, target.forecast.F_SQ)
50
test.segreg <- fwd(data.stf, sr=result.sr.segreg, target.forecast.F_SQ)
51