---
title: "Nonstationary FFA"
output: 
  rmarkdown::html_document:
    theme: readable
vignette: >
  %\VignetteIndexEntry{Nonstationary FFA}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

This vignette demonstrates how to use the `ffaframework` package to perform flood frequency analysis (FFA) under the assumption of nonstationarity.
Readers unfamiliar with stationary FFA workflows should first consult the *Stationary FFA* vignette.

The framework supports three forms of nonstationarity, modelled as linear trends in the distribution parameters:

1. A linear trend in the mean (location parameter).
2. A linear trend in the variance (scale parameter).
3. A linear trend in both the mean and variance.

The shape parameter is treated as a constant considering its difficult estimation due to sample size limitations.

## Case Study

This vignette will explore data from the Bow River at Banff (05BB001) hydrological monitoring station.
The remoteness of this station means that trends annual maxima are caused by changes in climate as opposed to changes in land use or cover.
Data for this station is provided as `Application_3.1.csv` in the `ffaframework` package.

```{r, fig.width = 10, fig.height = 8, fig.align = "center", out.width = "100%"}
library(ffaframework)

csv_path <- system.file("extdata", "Application_3.1.csv", package = "ffaframework")
df <- read.csv(csv_path, comment.char = "#")
df <- subset(df, !is.na(max)) # Remove missing values

head(df)

plot_ams_data(df$max, df$year, title = "Bow River at Banff (05BB001)")
```

### The `trend` List

This vignette assumes a scenario where a trend in the mean has been identified.
To specify a trend in the time series during FFA, create a `trend` list as shown below:

```{r}
trend <- list(location = TRUE, scale = FALSE)
```

**Note**: For guidance on trend detection, refer to the *Trend Identification* vignette.

## Distribution Selection

L-moment-based distribution selection remains applicable under nonstationarity, but requires dataset decomposition (detrending) prior to analysis.
This is accomplished using the `ams_decomposition` function, which takes an annual maximum series (AMS), the corresponding vector of years, and the `trend` object.
The decomposed vector of data is then passed to the selection function.

```{r, fig.width = 10, fig.height = 8, fig.align = "center", out.width = "100%"}
data_decomposed <- ams_decomposition(df$max, df$year, trend)

selection <- select_ldistance(data_decomposed)

print(selection$recommendation)

plot_lmom_diagram(selection)
```

**Conclusion**: The generalized normal (GNO) distribution is the best-fit distribution for the sample.

## Parameter Estimation

Because L-moments parameter estimation requires stationarity, we use maximum likelihood estimation for nonstationary models.
The `fit_maximum_likelihood` function implements maximum likelihood estimation for both stationary and nonstationary distributions.
It has two required arguments: 

- `data`: The annual maximum series observations. 
- `model`: A three-letter code corresponding to a probability distribution (ex. `"GNO"`).

Since a nonstationary model is used, two additional arguments are required:

- `years`: The corresponding vector of years for the observations in `data`.
- `trend`: The nonstationary `trend` object described above.

```{r}
fit <- fit_maximum_likelihood(
	df$max,
	"GNO",
	years = df$year,
	trend = trend
)

print(fit$params)

print(fit$mll)
```

**Note**: The fitted parameters are: $(\mu_0, \mu_1, \sigma, \kappa)$, where the time-dependent location is modeled as $\mu(t) = \mu_0 + \mu_1 t$.
The covariate $t$ is a linear function of the years: $(\text{years} - 1900) / 100$.

## Uncertainty Quantification

Uncertainty quantification is also essential for nonstationary probability distributions.
In addition to the parametric bootstrap, the framework implements the regula-falsi profile likelihood (RFPL) method for MLE.
The `uncertainty_rfpl` method has two required arguments:

- `data`: The annual maximum series observations. 
- `model`: A three-letter code corresponding to a probability distribution (ex. `"GNO"`).

Since the dataset has a nonstationary trend, three additional arguments are required:

- `years`: The corresponding vector of years for the observations in `data`.
- `trend`: The nonstationary `trend` object described above.
- `slices`: The years at which return levels are computed.

```{r, fig.width = 10, fig.height = 8, fig.align = "center", out.width = "100%"}
uncertainty <- uncertainty_rfpl(
	df$max,
	"GNO",
	years = df$year,
	trend = trend,
	slices = c(1925, 2025)
)

print(uncertainty[[2]]$estimates)

plot_nsffa(uncertainty)
```

**Example Conclusion**: In the year 2025, there is a $1/20$ exceedance probability of a flood of approximately $300\text{m}^3/\text{s}$ or greater.

**Note**: Under nonstationarity, the return period reflects the probability distribution for a fixed year rather than a long-run average.
To clarify this difference from stationary FFA, the phrase "effective return period" is used.
