%**************************************************************************
%Illustration of non-stationary extreme value modelling in MATLAB accompanying  
%Threshold modeling of nonstationary extremes
%Paul Northrop, Philip Jonathan and David Randell
%in 
%Extreme Value Modeling and Risk Analysis: Methods and Applications
%Dipak Dey and Jun Yan (Editors)
%2015
%**************************************************************************

%**************************************************************************
%Use this MATLAB script as a "driver" for execution of extreme value analysis
%**************************************************************************

%**************************************************************************
clc; clear; close all;
set(0,'DefaultAxesFontSize',18,'defaultFigureColor',[1,1,1]);
set(0, 'DefaultAxesBox', 'on');
set(0,'DefaultAxesXGrid','on','DefaultAxesYGrid','on');
%**************************************************************************

%**************************************************************************
%ANALYSIS PARAMETERS
%
%Spline basis
Prm.BssDmn=[0 360]';                       %2 x 1   spline basis domain (for covariate; must be in degrees for periodic)
Prm.BssIsPrd=1;                            %1 x 1   is periodic (==1 periodic, ==0 aperiodic)
Prm.BssnSgm=10;                            %1 x 1   number of (b-sline) segments
Prm.BssSplDgr=3;                           %1 x 1   spline degree (typically 3)
Prm.BssPnlOrd=1;                           %1 x 1   spline order (typically 1)
Prm.BssnDsp=360;                           %1 x 1   number of display locations (equally spaced in covariate domain) for parameter estimates
%Estimation of non-stationary value threshold
Prm.QRNEP=0.5;                             %1 x 1   non-exceedance probability for quantile regression estimation
%Estimation of Poisson rate of occurrence
Prm.PssnBin=40;                            %1 x 1   number of bins (for Poisson count estimation)
%Estimation of Poisson rate of occurrence
Prm.GPHlfWdt=50;                           %1 x 1   number of observations (sorted by covariate) to include in local estimate for GP starting solution (e.g. 30)
Prm.GPSclFct=3;                            %1 x 1   scale factor for roughness of GP scale (e.g. 3 means that the roughness of scale is 3^2 times the roughness of shape)
%Simulation of return values
Prm.RVnSct=8;                              %1 x 1   number of covariate sectors for display of extreme value analsis (e.g. 8 octants for directional covariate)
Prm.RVRtrPrdFct=1;                        %1 x 1   number of spline degree (typically 3)
Prm.RVnRls=100;                            %1 x 1   number of realisations for estimation of return value cumulative distribution functions
%Roughness penalisation, cross-validation and bootstrapping
Prm.Rgh=logspace(-6,6,10)';                %? x 1   values of parameter roughness to consider for all of QR, Pss and GP (specify single value to avoid cross-validation
Prm.CVnGrp=10;                             %1 x 1   "k" in k-fold cross-validation for QR, Pss and GP (e.g. 10; if LgrRgh is a scalar, cross-validation not preformed)
Prm.BSnRsm=10;                             %1 x 1   Number of bootstrap resamples (involving full model re-estimation) 
Prm.BSnRls=100;                            %1 x 1   Number of realisations of each bootstrap resample (used in return values simulation only) 
%**************************************************************************

%**************************************************************************
%GENERATE DATA (or use test deck provided)
%
if 0; %Generate toy data sample
    Dat=MakDat(Prm);
    save Dat Dat Prm;
else;
    Dat=load('Demo.mat');
end;
%
%Sample data structure "Dat" contains fields:
%  n       1 x 1 number of points
%  Y       n x 1 values of response
%  X       n x 1 values of covariate (all values must be in BssDmn)
%
%**************************************************************************

%**************************************************************************
%RUN FULL ANALYSIS
%
%Estimate extreme value model and simulate under estimated model
Rsl=KrnBS(Dat,Prm);
%Structure Rsl (result) has fields
%Rsl.X             q x 1  covariate values for display of QR, Pss and GP estimates (q=Prm.BssnDsp)
%Rsl.QR.MLE        q x 1  extreme value threshold maximum likelihood estimate
%Rsl.QR.BS         q x 3  extreme value threshold bootstrap estimates for median, lower 2.5%ile and upper 97.5%ile
%Rsl.Pss.MLE       q x 1  exceedance rate maximum likelihood estimate
%Rsl.Pss.BS        q x 3  exceedance rate bootstrap estimates for median, lower 2.5%ile and upper 97.5%ile
%Rsl.GP.Shp.MLE    q x 1  generalised Pareto shape maximum likelihood estimate
%Rsl.GP.Shp.BS     q x 3  generalised Pareto shape bootstrap estimates for median, lower 2.5%ile and upper 97.5%ile
%Rsl.RV.Omn        r x 1  (Not distributed) return value cumulative distribution function over full covariate domain (r=Prm.BSnRsm x Prm.BSnRpl)
%Rsl.RV.Sct        r x s  (Not distributed) return value cumulative distribution functions for covariate sectors (s=Prm.RVnSct)
%**************************************************************************

%**************************************************************************
%FIGURES
%
%Figure 1: Test deck visualisation
%Figure 2: Optimal parameter estimates
%Figure 3: Optimal roughness penalty coefficients
%Figure 4: (Not distributed) Return value cumulative distribution function for period of sample
%Figure 5: (Not distributed) Return value cumulative distribution function for period of sample
%**************************************************************************