Short_Term_Forecasting_for_advice_using_FLash.Rmd

---
title: 'Short Term Forecasting for advice using FLash'
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
github_document:
  mathjax: TRUE
pdf_document:
  fig_width: 6 
  fig_height: 4 
tags: [FLR]
license: Creative Commons Attribution-ShareAlike 4.0 International Public License
---

```{r, ini, echo=FALSE, results='hide', message=FALSE, warnings=FALSE, cache=FALSE}
library(knitr)
source("R/ini.R")
```
 
Short term forecasts help projecting the stock to forecast the implications of different management choices over a limited number of years. There are three necessary ingredients when projecting in FLR:

- a stock object that you will forecast and about which you have made some assumptions regarding what will happen in the future,
- a stock-recruitment relationship,
- a projection control specifying what the targets are and when to hit them.

In FLR, the method for short term forecasts is called `fwd()`, which takes an `FLStock`, `FLSR`, and `fwdControl`. 

## Required packages

To follow this tutorial you should have installed the following packages:

- CRAN: [ggplot2](https://cran.r-project.org/web/packages/ggplot2/index.html)
- FLR: [FLCore](http://www.flr-project.org/FLCore/),  [FLash](http://www.flr-project.org/FLash/), [FLAssess](http://www.flr-project.org/FLAssess/), [FLBRP](http://www.flr-project.org/FLBRP/),  [ggplotFL](http://www.flr-project.org/ggplotFL/)

if you are using Windows, please use the 32-bit R version

You can do so as follows,

```{r, eval=FALSE}
install.packages(c("FLAssess","FLBRP","ggplotFL"), repos="http://flr-project.org/R")
```

The example used here is for the ple4 data that is available in FLR. This needs to be loaded first.

```{r, eval=TRUE}
# Load the required packages
library(FLAssess)
library(FLash)
library(ggplotFL)
library(FLBRP)
# Load the data
data(ple4)
summary(ple4)
```

# Extending the stock object for the projections

Our ple4 stock goes up to `r range(ple4)[["maxyear"]]`, and in this example we want to make a 3-year projection, so we need to extend the stock by 3 years (`nyears`).

The projection will predict abundances in the future, but what  will the future stock weights, maturity, natural mortality, etc. be like? The assumptions about these will affect the outcome of the projection. Rather than using `window()` or `trim()` that make all future data *NA*, we use the `stf()` function (short term forecast) that has several options that allow one to control the assumptions about the future stock parameters.

These assumptions specify how many years you want to average over to set future values. For example, `wts.nyears` is the number of years over which to calculate the `*.wt`, `*.spwn`, `mat` and `m` slots. By default this is set to 3 years. This is a fairly standard assumption for short term forecasts, i.e. the future weights at age will be the same as the average of the last 3 years. This assumes that what is going to happen in the next few years will be similar to what happened in the last few years.

A simple 3-year forecast for the weights, natural mortality, etc., assuming these are equal to their averages over the last 3 years, is done by:

```{r, eval=TRUE}
maxyr_stk <- range(ple4)[["maxyear"]]
ple4_stf <- stf(ple4,nyears=3,wts.nyears=3, na.rm=TRUE)
maxyr_stf <- range(ple4_stf)[["maxyear"]]
```

Previously, the stock went up to `r maxyr_stk`. Now the stock goes up to `r maxyr_stf` and the future weights are equal to the average of the last three observed years, as can be observed from the last six years of the stock weights of the ple4_stf object.

```{r, eval=TRUE}
range(ple4_stf)
stock.wt(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]
```

For maturity the same assumption (of the future being the same as the average of the last 3 years) holds, but those maturities were assumed constant already.

```{r, eval=TRUE}
mat(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]
```

Notice that the  future fishing mortality has also been set (average of the last 3 years, by default).

```{r figA}
ggplot(harvest(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]) + geom_line(aes(x=age, y=data)) + facet_wrap(~year)
```

The stock numbers at age and the catch numbers at age are not forecast yet - this is what the `fwd()` function will perform later.
  
```{r, eval=TRUE}
stock.n(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]
```

Meanwhile, the landings and discards are average forecast ratios (proportion of total catch) of what is discarded and what is landed over the last 3 years of data.  

```{r, eval=TRUE}
discards.n(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]
landings.n(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]
# Compare above landings.n for the forecast years with
yearMeans((landings.n(ple4)/(landings.n(ple4)+discards.n(ple4)))[,ac((maxyr_stk-2):maxyr_stk)])
# Furthermore, landings and discards proportions sum to 1
landings.n(ple4_stf)[,ac((maxyr_stf-2):maxyr_stf)] + discards.n(ple4_stf)[,ac((maxyr_stf-2):maxyr_stf)]
```

# The stock-recruitment relationship (SRR)

A short term forecast does not use an SRR (in the traditional sense). Instead, it generally assumes that recruitment in the future is some mean (e.g. geometric mean) of the historic recruitments. However, we still need to have an SRR that contains this mean value, which is what we mimic for this example. First, we estimate the geometic mean recruitment, that we then add to an SRR object  

```{r, eval=TRUE}
mean_rec <- exp(mean(log(rec(ple4))))
ple4_sr <- as.FLSR(ple4, model="geomean")
params(ple4_sr)['a',] <- mean_rec
params(ple4_sr)
```

# The control object

The final thing we need to set up is the control object. This tells the projection what to do, i.e. what level of fishing mortality to use.

A standard scenario for a short term forecast is the 'status quo' scenario, assuming that the future mean fishing mortality will be the same as the mean of the last X years (X depending on the stock). We will set up this scenario as a simple example using the Fbar in the last year (`r range(ple4)[["maxyear"]]`), being `r round(mean(fbar(ple4)[,as.character(maxyr_stk)]),3)` per year.

```{r figB}
round(fbar(ple4),3)
ggplot(fbar(ple4), aes(x=year,y=data)) + geom_line()
```

Below, we define the last year of the stock and the status quo Fbar `fbar_SQ`.

```{r, eval=TRUE}
fbar_SQ <- mean(fbar(ple4)[,as.character(maxyr_stk)])
```

Now we introduce the control object: `fwdControl()`. This takes 1 argument - a data.frame that sets:

* The quantity (or type) of the target (we are using Fbar but can also set catch etc. - see medium term forecast tutorial)
* The value of the target, here status quo F 
* The year the target is to be hit 
* Some other things that we will ignore for now
 
Let's make the data.frame

```{r, eval=TRUE}
# Set the control object - year, quantity and value for the moment
ctrl_target <- data.frame(year = 2018:2020, quantity = "f", val = fbar_SQ)
ctrl_f <- fwdControl(ctrl_target)
ctrl_f
```

We see that we have what looks like our `ctrl_target`, but now it has two more columns (min and max). There is another table underneath which is for uncertainty, but that we ignore here for the short term forecast examples.

# Running the STF

Below we run a simple short term forecast (STF) with 'status quo' future fishing mortality. Remember we had to make assumptions about the future (weights, fishing mortality pattern, discard ratio, etc.).

This is done using `fwd()`, which takes three objects: the stock, the control object, the SRR. It returns an updated FLStock object. This update has forecast stock numbers, based on the recruitment, fishing mortality, and natural mortality assumptions. Because we now have stock numbers, we can calculate forecast ssb.

```{r, eval=TRUE}
ple4_sq <- fwd(ple4_stf, ctrl = ctrl_f, sr = ple4_sr)
```

Below, we check if the recruitment in the forecast stock indeed corresponds to the mean recruitment

```{r, eval=TRUE}
mean_rec
rec(ple4_sq)[,ac((maxyr_stf-5):maxyr_stf)]
```

Similarly, we check if the fbar in the forecast stock indeed corresponds to status quo fishing mortality.
 
```{r, eval=TRUE}
round(fbar_SQ,3)
round(fbar(ple4_sq)[,ac((maxyr_stf-5):maxyr_stf)],3)
```

The stock numbers are calculated using the recruitment and future mortality assumptions.

```{r, eval=TRUE}
stock.n(ple4_sq)[,ac((maxyr_stf-5):maxyr_stf)]
```

Note that the future harvest slot is different from the harvest from the one we set up, but that the selection pattern is the same: F at age has been multiplied by a constant value to give us the target Fbar value.
 
```{r, eval=TRUE}
round(harvest(ple4_stf)[,ac((maxyr_stf-2):maxyr_stf)],3)
round(harvest(ple4_sq)[,ac((maxyr_stf-2):maxyr_stf)],3)
harvest(ple4_stf)[,ac((maxyr_stf-2):maxyr_stf)] / harvest(ple4_sq)[,ac((maxyr_stf-2):maxyr_stf)]
```

The catch numbers come from the predicted abundance and harvest rates using the Baranov equation. The catches are then split into landings and discards using the average ratios computed by `stf()`.

```{r, eval=TRUE}
landings.n(ple4_sq)[,ac((maxyr_stf-5):maxyr_stf)]
discards.n(ple4_sq)[,ac((maxyr_stf-5):maxyr_stf)]
```

We can see the projection using a plot here

```{r, eval=TRUE}
plot(ple4_sq) 
```

# Short term forecast with many F scenarios

Typically when running STF you explore several different future F scenarios. The scenarios are based on 'F status quo', which we calculated above as the mean F of the last X years.

For a 3 year STF the F pattern is:

* year 1: fbar_status_quo
* year 2: fbar_status_quo * fbar_multiplier
* year 3: fbar_status_quo * fbar_multiplier

Note that year 1 is typically called the 'intermediate year' (in ICES, this would be the year the Expert Group meets to run the assessment, and status quo F for the intermediate year is a common assumption), and the fbar_multiplier is the same for years 2 and 3

We are going to run several STFs with different values for the fbar_multiplier

```{r, eval=TRUE}
fbar_multiplier <- seq(from = 0, to = 2, by = 0.2)
```

Next we are going to build a data.frame that creates these scenarios. Each column in the dataframe is a year, each row is a scenario. Note that if you project for more than 3 years you will need to add more columns / years to the matrix

```{r, eval=TRUE}
fbar_scenarios <- cbind(rep(fbar_SQ,length(fbar_multiplier)),
                        fbar_multiplier*fbar_SQ,
                        fbar_multiplier*fbar_SQ)
```

Another scenario we are interested in is $F_{0.1}$. We can calculate this using [FLBRP](http://flr-project.org/doc/Reference_points_for_fisheries_management_with_FLBRP.html) (or maybe you already have a value).

```{r, eval=TRUE}
f01 <- c(refpts(brp(FLBRP(ple4)))["f0.1","harvest"])
# Add the F0.1 scenario as a final scenario 
fbar_scenarios <- rbind(fbar_scenarios, c(fbar_SQ,f01,f01))
```

```{r, eval=TRUE}
# Add some names
colnames(fbar_scenarios) <- c("2018","2019","2020")
rownames(fbar_scenarios) <- c(fbar_multiplier, "f01")
fbar_scenarios
```

There are various results we want to extract from the STF, like predicted Catch, SSB and the relative change in these. First, make an empty matrix in which to store the results.

```{r, eval=TRUE}
stf_results <- matrix(NA,nrow = nrow(fbar_scenarios),ncol = 8)
# Set some column names
colnames(stf_results) <- c('Fbar',
    paste0('Catch',maxyr_stk+1), 
    paste0('Catch',maxyr_stk+2),
    paste0('Catch',maxyr_stk+3),
    paste0('SSB',maxyr_stk+2),
    paste0('SSB',maxyr_stk+3),
    paste0('SSB_change_',maxyr_stk+2,'-',maxyr_stk+3,'(%)'),
    paste0('Catch_change_',maxyr_stk,'-',maxyr_stk+2,'(%)'))
```

```{r, eval=TRUE, results="hide"}
# Set up an FLStocks object to store the resulting FLStock each time
stk_stf <- FLStocks()
# Loop over the scenarios (each row in the fbar_scenarios table)
for (scenario in 1:nrow(fbar_scenarios)) {
    cat("Scenario: ", scenario, "\n")
    flush.console()
    # Make a target object with F values for that scenario
    # Set the control object - year, quantity and value for the moment
    ctrl_target <- data.frame(year = (maxyr_stf-2):maxyr_stf,
                              quantity = "f",
                              val = fbar_scenarios[scenario,])
    # ctrl_target
    ctrl_f <- fwdControl(ctrl_target)
    # Run the forward projection. We could include an additional argument, maxF.
    # By default the value of maxF is 2.0. It could be increased to 10.0, say,
    # so that F is less limited, and the bound is not hit (not a problem here).
    ple4_fwd <- fwd(ple4_stf, ctrl = ctrl_f, sr = ple4_sr)#, maxF = 10.0)
    ## Check it has worked - uncomment out to check scenario by scenario
    # plot(ple4_fwd[,ac(2010:2020)])
    # Store the result - if you want to, comment out if unnecessary
    stk_stf[[as.character(scenario)]] <- ple4_fwd

    # Fill results table
    stf_results[scenario,1] <- round(fbar(ple4_fwd)[,ac(2020)],3) # final stf year
    stf_results[scenario,2] <- catch(ple4_fwd)[,ac(maxyr_stk+1)] # 1st stf year
    stf_results[scenario,3] <- catch(ple4_fwd)[,ac(maxyr_stk+2)] # 2nd stf year
    stf_results[scenario,4] <- catch(ple4_fwd)[,ac(maxyr_stk+3)] # final stf year
    stf_results[scenario,5] <- ssb(ple4_fwd)[,ac(maxyr_stk+2)] # 2nd stf year
    stf_results[scenario,6] <- ssb(ple4_fwd)[,ac(maxyr_stk+3)] # final stf year
    
    # change in ssb in last two stf years
    stf_results[scenario,7] <- round((ssb(ple4_fwd)[,ac(maxyr_stk+3)]-ssb(ple4_fwd)[,ac(maxyr_stk+2)])/
                                      ssb(ple4_fwd)[,ac(maxyr_stk+2)]*100,1) 
    
    # change in catch from true year, to 2nd to last stf year
    stf_results[scenario,8] <- round((catch(ple4_fwd)[,ac(maxyr_stk+2)]-catch(ple4_fwd)[,ac(maxyr_stk)])/
                                      catch(ple4_fwd)[,ac(maxyr_stk)]*100,1) 
}
# Give the FLStocks object some names
names(stk_stf) <- rownames(fbar_scenarios)

# Plotting
plot(stk_stf)
```
The different scenarios are plotted above. The table of results is given below.


```{r, eval=TRUE}
stf_results
```

# More information

* You can submit bug reports, questions or suggestions on this tutorial at <https://github.com/flr/doc/issues>.
* Or send a pull request to <https://github.com/flr/doc/>
* For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage, <http://flr-project.org>.

## Software Versions

* `r version$version.string`
* FLCore: `r packageVersion('FLCore')`
* FLAssess: `r packageVersion('FLAssess')`
* FLash: `r packageVersion('FLash')`
* ggplotFL: `r packageVersion('ggplotFL')`
* ggplot2: `r packageVersion('ggplot2')`
* FLBRP: `r packageVersion('FLBRP')`
* **Compiled**: `r date()`

## License

This document is licensed under the [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0) license.

## Author information
**Jan-Jaap Poos** Wageningen UR. Wageningen Marine Research. Haringkade 1, IJmuiden, The Netherlands.

**Katell Hamon** Wageningen UR. Wageningen Economic Research. Alexanderveld 5, The Hague, The Netherlands.