wiki-example-stockassessment/model_forecast.R at main · ices-taf-dev/wiki-example-stockassessment · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

library(stockassessment)
load("model/fit.RData")

# scenarios:

###################
## USER edits here
###################

Ry <- 2012:2022

Fpa <- 0.69
Flim <- 1.47
Fmsy <- 0.37
Fmsyadvice <- Fmsy * 15519 / 24739

#######################
## USER edits end here
#######################

forecast <- list()

# run forecast scenarios
set.seed(12345)
forecast[[length(forecast) + 1]] <-
  forecast(
    fit,
    fscale = c(1, 1, 1, 1),
    rec.years = Ry,
    processNoiseF = FALSE,
    label = "SQ all years",
    addTSB = TRUE
  )

set.seed(12345)
zeroF <- 0.000001
forecast[[length(forecast) + 1]] <-
  forecast(
    fit,
    fscale = c(1, NA, NA, NA),
    fval = c(NA, zeroF, zeroF, zeroF),
    rec.years = Ry,
    processNoiseF = FALSE,
    label = "SQ then zero F",
    addTSB = TRUE
  )

set.seed(12345)
forecast[[length(forecast) + 1]] <-
  forecast(
    fit,
    fscale = c(1, NA, NA, NA),
    fval = c(NA, Fpa, Fpa, Fpa),
    rec.years = Ry,
    processNoiseF = FALSE,
    label = "SQ then Fpa",
    addTSB = TRUE
  )

set.seed(12345)
forecast[[length(forecast) + 1]] <-
  forecast(
    fit,
    fscale = c(1, NA, NA, NA),
    fval = c(NA, Flim, Flim, Flim),
    rec.years = Ry,
    processNoiseF = FALSE,
    label = "SQ then Flim",
    addTSB = TRUE
  )

set.seed(12345)
forecast[[length(forecast) + 1]] <-
  forecast(
    fit,
    fscale = c(1, NA, NA, NA),
    fval = c(NA, Fmsy, Fmsy, Fmsy),
    rec.years = Ry,
    processNoiseF = FALSE,
    label = "SQ then Fmsy",
    addTSB = TRUE
  )

set.seed(12345)
forecast[[length(forecast) + 1]] <-
  forecast(
    fit,
    fscale = c(1, NA, NA, NA),
    fval = c(NA, Fmsyadvice, Fmsyadvice, Fmsyadvice),
    rec.years = Ry,
    processNoiseF = FALSE,
    label = "SQ then Fmsyadvice",
    addTSB = TRUE
  )

save(forecast, file = "model/forecast.RData")