\documentclass[handout]{beamer} \usepackage[latin1]{inputenc} \usepackage{amsmath} \usepackage{color} \usepackage{xcolor} \usepackage{graphicx} \usepackage{pifont} \usepackage{psfrag} \usepackage{skull} \usepackage{tikz} \usepackage{bm} %\usepackage{enumitem}% http://ctan.org/pkg/enumitem %\usetheme{Warsaw} \usecolortheme{ncl} \usefonttheme[onlymath]{serif} \newcommand{\penaltymark}[1]{\colorbox{nclRed}{\color{white}(-#1)}} \newcommand{\vsp}{\vspace{0.3cm}} \newcommand{\vspl}{\vspace{0.15cm}} \newcommand{\firstsect}{\setbeamercolor{structure}{fg=green!60!black}} \newcommand{\sol}{\setbeamercolor{structure}{fg=purple!80!black}} \newcommand{\normal}{\setbeamercolor{structure}{fg=blue!60!green}} \newcommand{\secondsect}{\setbeamercolor{structure}{fg=black!50!purple}} \newcommand{\thirdsect}{\setbeamercolor{structure}{fg=orange!90!yellow}} \newcommand{\fourthsect}{\setbeamercolor{structure}{fg=green!50!blue}} \newcommand{\fifthsect}{\setbeamercolor{structure}{fg=red!70!blue}} \newcommand{\cover}{\setbeamercolor{structure}{fg=blue!30!black}} \title[]{Predicting Halloween- and \\ Carrington-type solar flares} \author{Lee Fawcett, Holly Park, $\textsf{Thomai Tsiftsi}^{\dagger}$, $\textsf{Shaun Bloomfield}^{*}$} \date{\scriptsize $\textsf{\color{white}{L}}^{\dagger}$Universidad Nacional Aut\'{o}noma de M\'{e}xico, \\ $\textsf{\color{white}{L}}^{*}$Northumbria University} %\\ %\scriptsize $\textsf{\color{white}{L}}^{*}$Northumbria University} \definecolor{dark-gray}{gray}{0.425} \definecolor{light-gray}{gray}{0.125} \setbeamercovered{transparent=10} \begin{document} \begin{frame} \titlepage \end{frame} \begin{frame} \frametitle{Structure of this talk} \begin{itemize} \item Background / Motivation \begin{itemize} \item Carrington event \item Halloween solar storm \item Potential consequences \end{itemize} \pause \vspl \item Data and typical POT analysis \pause \vspl \item Improved analysis \begin{itemize} \item Threshold selection \item Dependence \item Prior specification \item Results / comparison \end{itemize} \end{itemize} \end{frame} \begin{frame}[plain] \frametitle{Carrington solar storm, 1859} \begin{itemize} \item Powerful \textcolor{purple}{\bf geomagnetic solar storm} during solar cycle 10 (1855-1867) \vspl \item Solar \textcolor{purple}{\bf coronal mass ejection} (CME) hit Earth's magnetosphere inducing one of the largest geomagnetic storms on record \vspl \item Observed and recorded by \textcolor{purple}{\bf Richard Carrington} \vspl \item Aurorae visible in the northern hemisphere as far south as Monterrey, Cuba and southern Japan \vspl \item \textcolor{purple}{\bf Telegraph systems} failed across Europe and North America \vspl \item Peak fluxes classified as A, B, C, M, X: each letter indicates a flare an order of magnitude larger than the one preceding it\pause \begin{itemize} \item [--] A1 flare has a peak X-ray flux of $1\times 10^{-8}\textrm{ Wm}^{-2}$ \item [--] X1 flare has a peak X-ray flux of $1\times 10^{-4}\textrm{ Wm}^{-2}$ \item [--] Carrington event: estimated as X45, i.e. $45\times 10^{-4}\textrm{ Wm}^{-2}$ \end{itemize} \end{itemize} \end{frame} \begin{frame} \frametitle{Carrington solar storm, 1859} \begin{center} \includegraphics[width=4.8cm,height=3.4cm]{flares1a.ps} \hspace{0.15cm} \includegraphics[width=4.8cm,height=3.4cm]{flares1c.ps} \vspace{0.15cm} \includegraphics[width=4.8cm,height=3.4cm]{flares1b.ps} \hspace{0.15cm} \includegraphics[width=4.8cm,height=3.4cm]{flares1d.ps} \end{center} \end{frame} \frame{ \frametitle{Halloween solar storms, 2003} \begin{itemize} \item Series of solar flares and CMEs occurring in October/November 2003 during solar cycle 23 (1996-2008) \vspl \item Generated the largest solar flare ever recorded by the \textcolor{purple}{\bf GOES} system, estimated at X28 \vspl \item \textcolor{purple}{\bf Satellite systems} and \textcolor{purple}{\bf communications} affected, numerous \textcolor{purple}{\bf spacecraft damaged}, \textcolor{purple}{\bf power outages} in Sweden and Finland \vspl \item Aurorae observed in Texas and southeast Mediterranean \end{itemize} } \frame{ \frametitle{Near miss: Carrington-like solar storm, 2012} \begin{minipage}{1.6\textwidth} \vspace{-0.13cm} %\hspace{-3.5cm} \hspace{-2cm}\includegraphics[width=17cm]{blast.ps} \end{minipage} } \frame{ \frametitle{What's at stake?} \begin{itemize} \item Disruption to \textcolor{purple}{\bf global positioning systems}. GPS industry estimated to be a \$1 trillion business \vspl \item \textcolor{purple}{\bf Satellite communications} at high risk \vspl \item Huge disruption to the \textcolor{purple}{\bf electrical grid}: areas such as the eastern US particularly vulnerable as power infrastructure highly interconnected, large cities without power for weeks \vspl \item Key questions: When might we see the next Carrington-like solar storm? What is the \textcolor{purple}{\bf return period} of a Carrington-like event? What is the strength of solar flare we might expect to see once per century? \end{itemize} } \frame{ \frametitle{Data} \begin{itemize} \item Extracted from the \textcolor{purple}{\bf Space Weather Prediction Center} (SWPC) online portal \vspl \item Recorded by the \textcolor{purple}{\bf Geostationary Operational Environmental Satellite} (GOES) system, operated by the US National Oceanic and Atmospheric Administration (NOAA) \vspl \item Data represent solar flares normalised from $H_{\alpha}$ spectral observations and X-ray fluxes, spanning a period of 43 years: November 1975-October 2017 \vspl \item Data points fairly regularly-spaced, with around four or five observations per day \vspl \item Data exhibit substantial short-term temporal dependence and cyclic variability \end{itemize} } \frame{ \frametitle{Data} \begin{center} \psfrag{W}[][]{\tiny $\textrm{Wm}^{-2} (\times 1000)$} \psfrag{j}[][]{\tiny \sf J ('03)} \psfrag{M}[][]{\tiny \sf M} \psfrag{F}[][]{\tiny \sf F} \psfrag{A}[][]{\tiny \sf A} \psfrag{J}[][]{\tiny \sf J} \psfrag{A}[][]{\tiny \sf A} \psfrag{S}[][]{\tiny \sf S} \psfrag{N}[][]{\tiny \sf N} \psfrag{jj}[][]{\tiny \sf J ('04)} \psfrag{jjj}[][]{\tiny \sf J ('05)} \psfrag{0.0}[][]{\tiny \sf 0.0} \psfrag{0.5}[][]{\tiny \sf } \psfrag{1.0}[][]{\tiny \sf 1.0} \psfrag{1.5}[][]{\tiny \sf} \psfrag{2.0}[][]{\tiny \sf 2.0} \psfrag{2.5}[][]{\tiny \sf} \psfrag{t}[][]{\tiny $t$} \psfrag{0.002}[][]{\tiny \sf .002} \psfrag{0.004}[][]{\tiny \sf } \psfrag{0.006}[][]{\tiny \sf .006} \psfrag{0.008}[][]{\tiny \sf } \psfrag{0.010}[][]{\tiny \sf .010} \psfrag{0.012}[][]{\tiny \sf } \psfrag{0}[][]{\tiny \sf 0} \psfrag{10000}[][]{\tiny \sf 10000} \psfrag{20000}[][]{\tiny \sf 20000} \psfrag{30000}[][]{\tiny \sf 30000} \psfrag{40000}[][]{\tiny \sf 40000} \psfrag{50000}[][]{\tiny \sf 50000} \psfrag{60000}[][]{\tiny \sf 60000} \includegraphics[width=11.45cm, height=7cm]{newplots1.ps} \end{center} } \frame{ \frametitle{Modelling extremes: Starting point} Let $\{X_{n}\}$ denote a stationary sequence of random variables with common distribution function $F$, and let $M_{n}=\textrm{max}\{X_{1}, \ldots, X_{n}\}$. \vsp It is typically the case that, as $n \to \infty$, \begin{eqnarray} \textrm{Pr}(M_{n} \leq x) &\approx& F^{n\theta}(x), \end{eqnarray} where $\theta \in (0,1)$ is known as the \textcolor{purple}{\bf extremal index}; see, for e.g., \textcolor{brown}{Leadbetter \& Rootz\'{e}n (1988)}. \vsp As $\theta \to 0$ there is increasing dependence in the extremes of the process; for an independent process, $\theta=1$. } \frame{ \frametitle{Modelling extremes: Block maxima} The \textcolor{purple}{\bf Generalised Extreme Value} (GEV) distribution is the limiting model for $F^{n}$, with distribution function (d.f.) \begin{eqnarray*} \mathcal{G}(y; \mu,\varsigma,\xi) &=&\left\{ \begin{array}{l c} \textrm{exp}\left[-\left(1+\xi(y-\mu)/\varsigma\right)^{-1/\xi}\right], &\xi \neq 0,\\ \textrm{exp}\left[-\textrm{exp}\left(-(y-\mu)/\varsigma\right)\right], & \xi=0, \end{array} \right. \end{eqnarray*} defined on $\{y: 1+\xi(y-\mu)/\varsigma>0\}$, where $-\infty < \mu < \infty$, $\varsigma>0$ and $-\infty < \xi < \infty$ are parameters of location, scale and shape. \vspl \begin{itemize} \item It can be shown that $F^{n\theta}$ is also GEV with d.f. $\mathcal{G}(y; \mu^{*},\varsigma^{*},\xi)$, provided long-range dependence is negligible\footnote{Technically, this requires \textcolor{purple}{Leadbetter's $D(u_{n})$ condition} to hold; see \textcolor{brown}{Leadbetter \textit{et al.} (1983)}.} \vspl \item In practice, the GEV is used to model maxima over some convenient calendar unit -- usually years \end{itemize} } \frame{ \frametitle{Modelling extremes: Threshold excesses} \textcolor{brown}{Pickands (1975)} showed that for $\theta=1$ and $u$ large, $(X-u|X>u)$ follows a \textcolor{purple}{\bf Generalised Pareto} (GP) distribution with d.f. \begin{eqnarray*} \mathcal{H}(y) &=& \left\{ \begin{array}{l c} 1-\left(1+\xi y/\sigma\right)^{-1/\xi}, &\xi \neq 0, \\ 1-\textrm{exp}\left[-y/\sigma\right], & \xi =0, \end{array} \right. \end{eqnarray*} defined on $\{y: y>0 \textrm{ \sf and } (1+\xi y/\sigma)>0\}$, where $\sigma=\varsigma+\xi(u-\mu)>0$ and $\xi$ are the scale and shape. \vspl \begin{itemize} \item Graphical diagnostics available for \textcolor{purple}{\bf selection of threshold} $u$, but often difficult to interpret \vspl \item Extremal index $\theta$ can't be absorbed into the GP parameters; careful consideration of \textcolor{purple}{\bf extremal dependence} required \end{itemize} } %An analysis of threshold excesses has the potential to increase estimation precision owing to the inclusion of more extremes.} \frame{ \frametitle{Modelling extremes: Return level estimation} The {\color{purple} $\boldsymbol{r}${\bf -year return level} $\boldsymbol{z_{r}}$} can then be obtained by inversion of $\mathcal{G}(z_{r})$ or $\mathcal{H}^{\theta}(z_{r})$. \pause \vsp For threshold excesses, on equating to $1-r^{-1}$ this gives \begin{eqnarray} z_{r} = \left\{ \begin{array}{l c} u+\sigma \xi^{-1}\left[\left(\lambda_{u}^{-1}w_{r}\right)^{-\xi}-1\right], &\xi \neq 0,\\ u-\sigma\textrm{ log}\left(\lambda_{u}^{-1}w_{r}\right), &\xi=0, \end{array} \right. \end{eqnarray} where $w_{r}=1-\left[1-(rn_{y})^{-1}\right]^{1/\theta}$, $\lambda_{u}$ is the rate of threshold excess and $n_{y}$ is the (average) number of observations per year. \pause \vsp \noindent For seasonally varying extremes, we usually solve \begin{eqnarray*} \prod_{m=1}^{M}\mathcal{H}_{m}^{\theta_{m}}(z_{r}) &=& 1-r^{-1}, \end{eqnarray*} where $\mathcal{H}_{m}$ is the GP distribution function in season $m$ with parameter set $(\lambda_{u_{m}}, \sigma_{m}, \xi_{m})$, and $\theta_{m}$ is the extremal index in season $m$. } \frame{ \frametitle{Cyclic variation} \begin{itemize} \item Close inspection of the time series suggests non-stationarity \textcolor{purple}{\bf inherited from the solar cycle}, affecting both the frequency and the severity of extreme flares \vspl \item After trial-and-improvement, fit sinusoidal waves of the form \[ \sum_{i=1}^{3} A_{i}\sin \omega_{i}t+B_{i}\cos \omega_{i}t + C \] to the series of moving averages, subtracting the estimated function before performing extreme value analyses \end{itemize} \begin{center} \psfrag{W}[][]{\tiny $\textrm{Wm}^{-2} (\times 1000)$} \psfrag{t}[][]{\tiny $t$} \psfrag{0}[][]{\tiny \sf 0} \psfrag{0.002}[][]{\tiny \sf .002} \psfrag{0.004}[][]{\tiny \sf } \psfrag{0.006}[][]{\tiny \sf .006} \psfrag{0.008}[][]{\tiny \sf } \psfrag{0.010}[][]{\tiny \sf .010} \psfrag{0.012}[][]{\tiny \sf } \psfrag{10000}[][]{\tiny \sf 10000} \psfrag{20000}[][]{\tiny \sf 20000} \psfrag{30000}[][]{\tiny \sf 30000} \psfrag{40000}[][]{\tiny \sf 40000} \psfrag{50000}[][]{\tiny \sf 50000} \psfrag{60000}[][]{\tiny \sf 60000} \includegraphics[width=7.6cm, height=3.5cm]{season.ps} \end{center} } \frame{ \frametitle{Typical peaks over threshold (POT) analysis} \begin{itemize} \item Identify independent clusters of solar flare extremes over some \textcolor{purple}{\bf suitably high threshold} $u$ \vspl \item Runs declustering, with \textcolor{purple}{\bf cluster termination interval} $\kappa$, most commonly used \vspl \item Fit the GP distribution to the set of independent cluster peak excesses \vspl \item Substitute $(\hat{\lambda}_{u}, \hat{\sigma}, \hat{\xi}, \theta=1)$ into \textcolor{brown}{(2)} to obtain $\hat{z}_{r} / \hat{r}$: \vspl {\scriptsize \begin{tabular}{|c c|c c|c c|} \hline \multicolumn{2}{|c|}{$(u, \kappa)$} &\multicolumn{2}{c|}{$\hat{z}_{50}$} &\multicolumn{2}{c|}{$\hat{r}$ years} \\ & &Post. mean &95\% CI &Halloween &Carrington \\ \hline $u = 90\%$ &$\kappa = 10$ &\textcolor{red}{X35} &\textcolor{red}{(X31, X50)} &\textcolor{red}{38} &\textcolor{red}{110} \\ quantile &$\kappa = 15$ &X35 &(X31.5, X52) &37.5 &100 \\ &$\kappa = 30$ &X37 &(X32, X65) &33 &96.5 \\ \hline $u = 95\%$ &$\kappa = 10$ &X35 &(X30, X46) &38 &100 \\ quantile &$\kappa = 15$ &X36 &(X31, X47.5) &37 &97 \\ &$\kappa = 30$ &X38.5 &(X30, X72.5) &30 &79 \\ \hline \end{tabular}} {\tiny \textcolor{red}{(As reported in Tsiftsi and De La Luz, 2018)}} \end{itemize} } \frame{ \frametitle{Typical peaks over threshold (POT) analysis} \begin{itemize} \item Choice of threshold difficult: \begin{itemize} \item Just choose a high quantile? \item Graphical diagnostics -- e.g. \textcolor{purple}{\bf mean residual life plot}, \textcolor{purple}{\bf Hill plots}, \textcolor{purple}{\bf parameter stability plots}... \item Once $u$ is chosen, it is treated as fixed \end{itemize} \pause \vspl \item Choice of declustering scheme difficult: \begin{itemize} \item Runs declustering: Choice of $\kappa$? \item Other filtering schemes available \item Wasteful -- \textcolor{purple}{\bf large uncertainty} on return level estimates \end{itemize} \pause \vspl \item Prior specification -- objective priors, missed opportunity? \end{itemize} } \frame{ \frametitle{Improved analysis: Aims} Here, we have three main aims: \pause \vspl \begin{enumerate} \item Estimate the \textcolor{purple}{\bf extremal index} for use in return level estimation via a Markov chain model for threshold excesses \begin{itemize} \item Avoids declustering altogether, and the problems it brings \item Presses into use information on \textit{all} extremes \end{itemize} \pause \vspl \item Capture uncertainty in the \textcolor{purple}{\bf choice of threshold}, and enable a more `automatic' selection of $u$ \pause \vspl \item Get more out of working within the Bayesian framework \begin{itemize} \item Where possible, formulate \textcolor{purple}{\bf informative priors} \item Make use of the \textcolor{purple}{\bf posterior predictive distribution} \end{itemize} \end{enumerate} } \frame{ \frametitle{Improved analysis: Mixture model} Suppose our data comprise a sequence of $n$ identically distributed observations $\boldsymbol{x} = \{x_{i}: i=1, \ldots, n\}$. We assume the following mixture model: \[ f(x) = (1-\lambda_{u}) \cdot f_{1}(x) + \lambda_{u}\cdot f_{2}(x), \] where \begin{eqnarray*} f_{1}(x) &=& \frac{h(x| \phi, \boldsymbol{x})}{H(u| \phi, \boldsymbol{x})}I_{(-\infty, u]}(x) \quad \textsf{and}\\ & & \\ f_{2}(x) &=& g(x| u, \sigma_{u}, \xi)I_{(u, \infty)}(x). \end{eqnarray*} \pause \begin{itemize} \item Observations sub-threshold follow a non-parametric density $h$ \item Excesses above the threshold follow a GP distribution with parameter set $(u, \sigma_{u}, \xi)$ and extremal index $\theta$ \item Possibility of a discontinuity at $u$ \end{itemize} } \frame{ \frametitle{Improved analysis: Bulk distribution} For $f_{1}(x)$ we use a kernel estimator of the form \begin{eqnarray} \hat{f}(x) &=& n^{-1}\sum_{i=1}^{n}K_{\phi}(x-x_{i}), \end{eqnarray} where $K_{\phi}(y)=\phi^{-1}K(y/\phi)$ is chosen to be a unimodal probability density function that is symmetric about zero. \pause \vspl \begin{itemize} \item Kernel function used in \textcolor{brown}{(3)} not critical; we use a zero mean Gaussian kernel \item GP distribution is used for the upper tail -- extrapolation into the upper tails is of no concern \item \textcolor{purple}{\bf Kernel bandwidth likelihood} $L_{\textsf{\tiny K}}(\phi, u|\boldsymbol{x})$, for $A=\{i: x_{i}\leq u\}$: \begin{equation*} \left\{\frac{(1-\lambda_{u})}{n^{-1}\sum_{k=1}^{n} \Phi\left(\frac{u-x_{i}}{\phi}\right)}\right\}^{|A|} \prod_{A}\frac{1}{n-1}\sum_{j=1; j \neq i}^{n}K_{\phi}(x_{i}-x_{j}) \end{equation*} \end{itemize} } \frame{ \frametitle{Improved analysis: Threshold excess distribution} We assume a first-order Markov structure, giving the likelihood: \begin{equation} \prod_{i=1}^{n-1}f(x_{i},x_{i+1};\boldsymbol{\psi})\Bigg/ \prod_{i=2}^{n-1}f(x_{i}; \boldsymbol{\psi}). \end{equation} \pause \vspl \begin{itemize} \item Transformation from GP to unit Fr\'{e}chet margins gives a range of bivariate models to use for the numerator in \textcolor{brown}{(4)}\pause \vspl \item Most commonly-used is the \textcolor{purple}{\bf logistic family} with $\alpha \in [0,1)$: \begin{eqnarray} G(x_{i}, x_{i+1}) &=& \exp\left\{-\left(x_{i}^{-1/\alpha}+x_{i+1}^{-1/\alpha}\right)^{\alpha}\right\}. \end{eqnarray} \pause \vspl \item \textcolor{purple}{\bf Extreme value Markov chain likelihood} $L_{\textsf{\tiny EV}}(u, \sigma_{u}, \xi, \alpha| \boldsymbol{x})$: \begin{itemize} \item Differentiation of \textcolor{brown}{(5)}, with careful censoring, gives pairwise contributions to the numerator in \textcolor{brown}{(4)}\pause \item Univariate contributions to the denominator are given through the GP distribution \end{itemize} \end{itemize} } \frame{ \frametitle{Aside 1: Estimating the extremal index?} \begin{eqnarray*} z_{r} = \left\{ \begin{array}{l c} u+\sigma \xi^{-1}\left[\left(\lambda_{u}^{-1}w_{r}\right)^{-\xi}-1\right], &\xi \neq 0,\\ u-\sigma\textrm{ log}\left(\lambda_{u}^{-1}w_{r}\right), &\xi=0, \end{array} \right. \end{eqnarray*} where $w_{r}=1-\left[1-(rn_{y})^{-1}\right]^{1/\theta}$, $\lambda_{u}$ is the rate of threshold excess and $n_{y}$ is the (average) number of observations per year. } \frame{ \frametitle{Aside 1: Estimating the extremal index?} \noindent Define (arbitrarily) $x_{n}$ such that $F^{n}(x_{n}) = 1/2$ in \textcolor{brown}{(1)}. \pause \vsp \noindent Then, using \textcolor{brown}{(1)} \textcolor{gray}{[$\text{Pr}(M_n \leq x) \approx F^{n\theta}(x)$]}, we can define \begin{eqnarray} \theta_{n} &= &-\frac{\textrm{log } \textrm{Pr}\left(\textrm{max}\{X_{1}, \ldots, X_{n}\} \leq x_{n}\right)}{\textrm{log }2}, \end{eqnarray} \pause and so $\theta_{n} \to \theta$ as $n \to \infty$. \pause \vsp We can use \textcolor{brown}{(6)} to investigate the relationship between $\theta$ and $\alpha$ in the logistic model: \vspl \begin{itemize} \item Simulate $N$ first-order Markov chains, each of length $n$, with logistic dependence $\alpha$ \vspl \item The probability in the numerator of \textcolor{brown}{(6)} is estimated as the proportion of simulated chains whose maximum does not exceed $x_{n}$ \vspl \item We find values for $\theta_{n}$ for $\alpha = \{0.05, 0.10, ..., 0.95\}$, with $N=n=10^{6}$ \end{itemize} } \frame{ \frametitle{Aside 1: Estimating the extremal index?} %\begin{center} \psfrag{0}[][]{\sf \tiny 0} \psfrag{0.0}[][]{\sf \tiny 0.0} \psfrag{0.2}[][]{\sf \tiny 0.2} \psfrag{0.4}[][]{\sf \tiny 0.4} \psfrag{0.6}[][]{\sf \tiny 0.6} \psfrag{0.8}[][]{\sf \tiny 0.8} \psfrag{1.0}[][]{\sf \tiny 1.0} \psfrag{2}[][]{\sf \tiny 2} \psfrag{4}[][]{\sf \tiny 4} \psfrag{6}[][]{\sf \tiny 6} \psfrag{8}[][]{\sf \tiny 8} \psfrag{theta}[][]{\small $\theta$} \psfrag{r}[][]{\small $\rho$} \psfrag{alpha}[][]{\small $\alpha$} \psfrag{beta (alpha=0.6)}[][]{\small $\beta$} \psfrag{Logistic}[][]{\scriptsize \sf Logistic} \psfrag{Negative logistic}[][]{\scriptsize \sf Negative logistic} \psfrag{Bilogistic}[][]{\scriptsize \sf Bilogistic} \psfrag{relationship between alpha and theta}[][]{\sf \small Simulated values of extremal index $\theta$ for given logistic dependence $\alpha$} %\centering \includegraphics[width=12cm, height=5cm]{../PAPERS_2015/figure3.eps} %\end{center} \pause Giving, for example, \[ \theta \approx 0.013 - 0.092\alpha + 1.833\alpha^{2}-0.756\alpha^{3} \] } \frame{ \frametitle{Aside 1: Estimating the extremal index?} \psfrag{Relationship between return level and the extremal index}[][]{\sf \small Relationship between the 50--year and 200--year return level and the extremal index $\theta$} \psfrag{0}[][]{\sf \tiny 0} \psfrag{0.0}[][]{\sf \tiny 0.0} \psfrag{0.2}[][]{\sf \tiny 0.2} \psfrag{0.3}[][]{\sf \tiny 0.3} \psfrag{0.4}[][]{\sf \tiny 0.4} \psfrag{0.5}[][]{\sf \tiny 0.5} \psfrag{0.6}[][]{\sf \tiny 0.6} \psfrag{0.7}[][]{\sf \tiny 0.7} \psfrag{0.8}[][]{\sf \tiny 0.8} \psfrag{0.9}[][]{\sf \tiny 0.9} \psfrag{1.0}[][]{\sf \tiny 1.0} \psfrag{2}[][]{\sf \tiny 2} \psfrag{4}[][]{\sf \tiny 4} \psfrag{6}[][]{\sf \tiny 6} \psfrag{8}[][]{\sf \tiny 8} \psfrag{T}[][]{\scriptsize $\theta(u_{p})$} \psfrag{A}[][]{\tiny $A$} \psfrag{u9}[][]{\scriptsize \sf $\textit{AR}(1)$ with $p=0.9$} \psfrag{u95}[][]{\scriptsize \sf $\textit{AR}(1)$ with $p=0.95$} \psfrag{u99}[][]{\scriptsize \sf $\textit{AR}(1)$ with $p=0.99$} \centering \includegraphics[width=12cm, height=5cm]{../PAPERS_2015/fig8.ps} } \frame{ \frametitle{Improved analysis: Likelihood} Combining the Kernel density model for sub-threshold observations with the GP extreme value Markov chain for threshold excesses gives the \textcolor{purple}{\bf extreme value mixture model} likelihood: \[ L(\boldsymbol{\psi}|\boldsymbol{x}) = L_{\textsf{\tiny K}}(\phi, u|\boldsymbol{x}) \cdot L_{\textsf{\tiny EV}}(u, \sigma_{u}, \xi, \alpha|\boldsymbol{x}), \] where $\boldsymbol{\psi} = (\phi, u, \sigma_{u}, \xi, \alpha)$. } \frame{ \frametitle{Improved analysis: Prior specification} The joint prior for $\boldsymbol{\psi}$ is expressed as \[ \pi(\phi, u, \sigma_{u}, \xi, \alpha) = \pi(\phi) \cdot \pi(u) \cdot \pi(\sigma_{u}, \xi) \cdot \pi(\alpha) \] \begin{enumerate} \item [1. ] For the \textcolor{purple}{\bf kernel density bandwidth parameter} $\phi$ we follow \textcolor{brown}{Brewer (2000)} by specifying the prior for the precision $1/\phi^{2}$ as inverse gamma: \[ \pi(\phi|d_{1}, d_{2}) = \frac{1}{d_{2}^{d_{1}}\Gamma(d_{1})}\left(\frac{1}{\phi^{2}}\right)^{d_{1}-1} \exp\left(-\frac{1}{\phi^{2}d_{2}}\right), \] \pause with $d_{1}$ and $d_{2}$ chosen to give a largely non-informative prior specification. \vspl \item [2. ] In the absence of any expert prior information regarding the \textcolor{purple}{\bf logistic dependence parameter} $\alpha$, we assume $\alpha \sim \mathcal{U}(0,1)$. \end{enumerate} } \frame{ \frametitle{Improved analysis: Prior specification} \begin{enumerate} \item [3. ]Following \textcolor{brown}{Behrens et al. (2004)} we assume a truncated Normal prior for the \textcolor{purple}{\bf threshold} $u$ (truncated below $e_{1}$): \[ \pi(u|\mu_{u}, v_{u}^{2}) = \frac{1}{\sqrt{2 \pi v_{u}^{2}}} \frac{\exp\{-0.5[(u-\mu_{u})/v_{u}]^{2}\}}{\Phi[-(e_{1}-\mu_{u})/v_{u}]}, \] where $\mu_{u}$ is set at some high data percentile; $v_{u}^{2}$ is set sufficiently large to represent a diffuse prior. \end{enumerate} } \frame{ \frametitle{Improved analysis: Prior specification} \begin{enumerate} \item [4. ] For the \textcolor{purple}{\bf GP parameters} $(\sigma_{u}, \xi)$ we attempt to formulate informative priors. \begin{itemize} \item We work with a re-parameterised GP distribution: $\eta = \textrm{log}(\sigma_{u}-\xi u)$ \pause \vspl \item Following \textcolor{brown}{Coles \& Tawn (2005)} we assume \begin{eqnarray*} \tilde{z}_{1} &=&z_{11}-e_{1} \sim \textrm{Ga}(g_{1},h_{1}) \\ \tilde{z}_{2} &=&z_{100}-z_{11} \sim \textrm{Ga}(g_{2},h_{2}), \end{eqnarray*} where $(g_{\cdot},h_{\cdot})$ reflect the opinions of an expert\footnote{\tiny Shaun Bloomfield, Department of Mathematics, Physics and Electrical Engineering, Northumbria University} \item Then: \pause \[ \pi(z_{11},z_{100}) \propto z_{11}^{g_{1}-1}\textrm{e}^{-h_{1}z_{11}}(z_{100}-z_{11})^{g_{2}-1}\textrm{e}^{-h_{2}(z_{100}-z_{11})} \] \item Multiplying by the Jacobian of the transformation $(z_{11},z_{100})\longrightarrow(\eta,\xi)$ gives a prior in terms of $(\eta,\xi)$\pause \end{itemize} \pause \end{enumerate} } \frame{ \frametitle{Improved analysis: Prior specification} {\tiny \begin{center} \begin{tabular}{|l|| l|} \hline Holly &Shaun \\ \hline What classification of flare would you expect & \\ to occur only once per solar cycle (or once every &X20 \\ 11 or so years)? & \\ \hline What range of classifications would you expect & \\ the vast majority (say 95\%) of extreme solar & \\ flares occurring only once in one solar cycle &Confident to give a lower limit \\ to lie in? (i.e in a given solar cycle, about 95\% &of $\sim$X5 but an upper threshold \\ of the once-in-a-cycle flares will lie between &an would be unsubstantiated \\ X... class flare and an X... class flare) & \\ \hline What classification of flare would you & \\ expect to occur only once every hundred years? &X40 \\ \hline What range of classifications would you expect & \\ the vast majority (say 95\%) of extreme solar & \\ flares occurring only once in a hundred years &Again, simply a lower bound would \\ to lie in? (i.e in a one hundred year period, &be most appropriate here, and this \\ about 95\% of the once-in-a-hundred-year flares &can be given as $\sim$X30 \\ will lie between an X... class flare and an & \\ X... class flare) & \\ \hline Could you also answer these questions for the &For 30 years, $\sim$X30, for 150 years, \\ 30-year, 150-year and 300-year scenarios &$\sim$X50 and for 300 year, $\sim$X60 \\ \hline How often would you expect to see an event of & \\ the same calibre as the Halloween event of 2003 &Every $\sim$35-50 years \\ (an X35 event)? & \\ \hline How often would you expect to see an event of & \\ the same calibre as the Carrington event of 1859 &Every $\sim$100-150 years \\ (an X45 event)? & \\ \hline \end{tabular} \end{center} } } \frame{ \frametitle{Improved analysis: Prior specification} \textbf{Constructing gamma priors:} {\small \begin{center} \begin{tabular}{|c c |c c| c c|}\hline $i$ &$r_{i}$ &Expected &95\% Lower Bound &$g_i$ &$h_i$ \\ \hline 1 &11 &20 &5 &2.40 &1012.80 \\ 2 &100 &40 &30 &36.33 &9001.01 \\ \hline \end{tabular} \end{center} } \vspace{1cm} \textbf{Checking coherence:} \vspace{-0.3cm} {\small \begin{center} \begin{tabular}{|c |c|c c c|}\hline Return Period &Elicited Returns &\multicolumn{3}{c|}{Derived Returns} \\ (years) &Expected &Estimate &Predictive &95\% CI\\ \hline 30 &30 &25.74 &25.78 &(20.81, 31.31) \\ 150 &50 &54.60 &55.14 &(41.21, 73.52) \\ 300 &60 &74.90 &76.30 &(53.84, 105.24) \\ \hline \end{tabular} \end{center} } } \frame{ \frametitle{Improved analysis: Posterior sampling scheme} We use a simple Metropolis-Hastings scheme with random walk updates, tuning the innovation variances to optimise efficiency, giving successive draws \[ \boldsymbol{\psi}^{[j]} = (\phi^{[j]},u^{[j]}, \eta^{[j]}, \xi^{[j]}, \alpha^{[j]})^{T}, \quad j=1, \ldots, 10^{6}. \] \begin{itemize} \item Chains thinned by a factor of 50 to remove substantial autocorrelation \vspl \item Convergence assessed using some standard diagnostics as discussed in \textcolor{brown}{Gelman and Rubin (1992)} \vspl \item For threshold excesses, posterior draws for $(u, \sigma_{u}=f(u, \eta, \xi), \xi, \theta=f(\alpha))^{[j]}$ used to obtain posterior draws $z_{r}^{[j]}$ via \textcolor{brown}{(2)} \vspl \item This leads to standard reporting of, e.g., \textcolor{purple}{\bf posterior mean} plus \textcolor{purple}{\bf posterior st. dev.} / \textcolor{purple}{\bf credible interval} \end{itemize} } \frame{ \frametitle{Improved analysis: Predictive inference} \begin{itemize} \item Bayesian framework offers natural extension to prediction through the \textcolor{purple}{\bf posterior predictive distribution} \pause \vspl \item Let $Y$ denote a future extreme of our solar flares process. Then we can write \begin{eqnarray} \textsf{Pr}(Y\leq y|\boldsymbol{x}) &=& \int_{\boldsymbol{\Psi}}\textsf{Pr}(Y \leq y|\boldsymbol{\psi})\pi(\boldsymbol{\psi}|\boldsymbol{x}) d\boldsymbol{\psi} \end{eqnarray} \item Equating to $1-r^{-1}$ and solving for $y=z_{r, \textsf{\tiny pred}}$ gives a return level estimate which incorporates uncertainty due to parameter estimation \pause \vspl \item \textcolor{brown}{(7)} analytically intractable; approximate using \[ \frac{1}{S}\sum_{j=1}^{S} \textsf{Pr}\left\{Y\leq z_{r, \textsf{\tiny pred}}|\boldsymbol{\psi}^{[j]}\right\} = \pause \frac{1}{S} \sum_{j=1}^{S} \mathcal{H}^{\theta^{[j]}}(z_{r, \textsf{\tiny pred}}|u^{[j]}\sigma_{u}^{[j]}, \xi^{[j]}) \] where $S$ is the number of MCMC iterations after burn-in \end{itemize} } \frame{ \frametitle{Aside 2: Summarising our return level inference?} \psfrag{a0.3}[][]{\sf \tiny $\alpha=0.3$} \psfrag{a0.5}[][]{\sf \tiny $\alpha=0.5$} \psfrag{a0.9}[][]{\sf \tiny $\alpha=0.9$} \psfrag{A0.7}[][]{\sf \tiny $A=0.7$} \psfrag{A0.5}[][]{\sf \tiny $A=0.5$} \psfrag{A0.3}[][]{\sf \tiny $A=0.3$} \psfrag{A}[][]{\sf \tiny $A$} \psfrag{r}[][]{\sf \tiny $r$} \psfrag{P}[][]{\sf \tiny $p_{*}$} \psfrag{0.00}[][]{\sf \tiny 0.00} \psfrag{0.01}[][]{\sf \tiny 0.01} \psfrag{0.02}[][]{\sf \tiny 0.02} \psfrag{0.03}[][]{\sf \tiny 0.03} \psfrag{0.04}[][]{\sf \tiny 0.04} \psfrag{0.05}[][]{\sf \tiny 0.05} \psfrag{100}[][]{\sf \tiny 100} \psfrag{1000}[][]{\sf \tiny 1000} \psfrag{10000}[][]{\sf \tiny 10000} \includegraphics[width=11.7cm, height=8.25cm]{../PAPERS_2015/fig5.ps} } \frame{ \frametitle{Improved analysis: Results} \begin{center} \begin{tabular}{|c|c|c|} \hline &$z_{50}$ &$z_{1000}$ \\ \hline Post. mean &X35 &X57 \\ &\textcolor{gray}{\it X35} &\textcolor{gray}{\it X54} \\ 95\% CI &(X32, X41) &(X49, X66) \\ &\textcolor{gray}{\it (X31, X50)} &\textcolor{gray}{\it (X40, X88)} \\ \hline \textcolor{red}{Predictive} &\textcolor{red}{X39} &\textcolor{red}{X79} \\ \hline \hline &Halloween &Carrington \\ \hline Post. mean &37 years &93 years \\ &\textcolor{gray}{\it 38 years} &\textcolor{gray}{\it 110 years} \\ 95\% CI &(31, 45) years &(80, 114) years \\ &\textcolor{gray}{\it (22, 67) years} &\textcolor{gray}{\it (83, 239) years} \\ \hline \textcolor{red}{Predictive} &\textcolor{red}{43 years} &\textcolor{red}{112 years} \\ \hline \end{tabular} \end{center} } \frame{ \frametitle{Conclusions} \begin{itemize} \item Fitted a model which \textcolor{purple}{\bf doesn't require a choice of threshold / declustering scheme} \begin{itemize} \item Graphical diagnostics for threshold choice are often difficult to interpret / contradictory \item Uncertainty in threshold choice is accounted for \item Return level estimates can be sensitive to the choice of declustering scheme \item Using declustered extremes can result in low precision on return level/period estimates \end{itemize} \vspl \item Is the model for dependence the `best'? \begin{itemize} \item We used a model which implies asymptotic dependence \item Even if this a valid assumption, there are several choices of bivariate extreme value model \item Model order? \end{itemize} \vspl \item Unsure about the effect of the bulk distribution \pause \vspl \item \textcolor{purple}{\bf Directional modelling}? \end{itemize} } \end{document}