---
title: "Bayesian Likelihood-Based Inference for Time Series Extremes"
author: "Paul Northrop"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Bayesian Likelihood-Based Inference for Time Series Extremes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: lite.bib
csl: taylor-and-francis-chicago-author-date.csl
---

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

required <- c("exdex", "revdbayes")

if (!all(unlist(lapply(required, function(pkg) requireNamespace(pkg, quietly = TRUE)))))
  knitr::opts_chunk$set(eval = FALSE)
```

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

For information about the model underlying the inferences performed in this vignette, including the (adjusted) likelihood for $(p_u, \sigma_u, \xi, \theta)$ and the Cheeseboro wind gusts dataset used below, see the vignette [Frequentist Likelihood-Based Inference for Time Series Extremes](lite-1-frequentist.html). 

## Bayesian inferences for model parameters

We perform Bayesian inference for $(p_u, \sigma_u, \xi, \theta)$, combining a likelihood with a prior distribution for these parameters. For $\theta$ the likelihood is based on the $K$-gaps model. For $(p_u, \sigma_u, \xi)$ the likelihood is based on vertically-adjusted log-likelihoods for $p_u$ and  $(\sigma_u, \xi)$.  This latter aspect is an example of Bayesian inference using a composite likelihood (@RCD2012).  

Currently, **lite** allows only priors where $p_u$, $(\sigma_u, \xi$) and $\theta$ are independent a priori.  In the example below, we use the `blite` function's default priors for the exceedance probability $p_u$ (the Jeffreys' prior), the generalised Pareto parameters $(\sigma_u, \xi)$ (and maximal data information prior) and $\theta$ (a uniform prior on [0,1]).  Different priors can be set using the respective arguments `gp_prior`, `b_prior` and `theta_prior_pars`.

The `blite` function samples from the posterior distribution for $(p_u, \sigma_u, \xi)$ based on vertically-adjusted log-likelihoods for $p_u$ and $(\sigma_u, \xi)$.  To use unadjusted log-likelihoods set the argument `type = "none"`.

```{r cheeseboro}
library(lite)
cdata <- exdex::cheeseboro
# Each column of the matrix cdata corresponds to data from a different year
# blite() sets cluster automatically to correspond to column (year)
cpost <- blite(cdata, u = 45, k = 3, ny = 31 * 24, n = 10000)
```

The `summary` and `plot` methods produce numerical and graphical marginal summaries of the posterior samples.

```{r vertical, fig.width = 7, fig.height = 5, message = FALSE}
summary(cpost)
plot(cpost)
```

## Predictive inference

We perform posterior predictive inference for the largest value $M_N$ to be observed over a future time period of length $N$ years.  Objects returned from the `blite` function have a `predict` method based on the `predict.evpost` method in the **revdbayes** package [@revdbayes]. See the [Posterior Predictive Extreme Value Inference using the revdbayes Package](https://cran.r-project.org/package=revdbayes) vignette for information. The effect of adjusting for the values of the extremal index \eqn{\theta} in the posterior sample is to change the effective time horizon from $N$ to $\theta N$. 

The function `predict.blite` can estimate predictive intervals, quantiles, distribution and density functions for $M_N$ and simulate from the predictive distribution for $M_N$.  For example, the following code estimates a 95\% highest predictive density (HPD) interval for $M_{100}$ and plots the predictive densities of $M_{100}$ and $M_{1000}$.

```{r predinterval}
# Interval estimation
predict(cpost, hpd = TRUE)$short
```

```{r preddensity, fig.width = 7, fig.height = 5, message = FALSE}
# Density function
plot(predict(cpost, type = "d", n_years = c(100, 1000)))
```

## References

<script type="text/x-mathjax-config">
   MathJax.Hub.Config({  "HTML-CSS": { minScaleAdjust: 125, availableFonts: [] }  });
</script>