From b0f5e12d8d6644aeee2457962ec02101f3d30772 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 30 Jan 2026 16:42:42 +0000 Subject: [PATCH 1/4] Initial plan From b2c208d4077e7d3effe47ff6614f246bbcfd176d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 30 Jan 2026 16:51:50 +0000 Subject: [PATCH 2/4] Add comprehensive tsEVA 2.0 function reference documentation Co-authored-by: menta78 <19606593+menta78@users.noreply.github.com> --- docs/0_README.md | 87 +++ docs/2_Function_Reference.md | 1331 ++++++++++++++++++++++++++++++++++ 2 files changed, 1418 insertions(+) create mode 100644 docs/0_README.md create mode 100644 docs/2_Function_Reference.md diff --git a/docs/0_README.md b/docs/0_README.md new file mode 100644 index 0000000..3cf889a --- /dev/null +++ b/docs/0_README.md @@ -0,0 +1,87 @@ +# tsEVA 2.0 Documentation Navigation Guide + +**Purpose**: This document serves as the routing map for the tsEVA 2.0 MATLAB Expert AI Assistant (Tessa M). It defines what documentation must be read at the start of every session and what should be consulted on-demand. + +## Source of Truth + +This file (`0_README.md`) is the **authoritative routing document**. If any guidance conflicts between documentation files, follow the instructions in this file. + +## Mandatory Reading (Start of Every Session) + +Every new chat session MUST read these files first: + +1. **[1_Core_Methodology.md](1_Core_Methodology.md)** - The Transformed-Stationary (TS) approach, fundamental concepts, and theoretical foundations +2. **[2_Function_Reference.md](2_Function_Reference.md)** - Complete list of documented tsEVA 2.0 functions with signatures and descriptions +3. **[3_Workflow_Patterns.md](3_Workflow_Patterns.md)** - Standard analysis workflows and common patterns + +## Read Only When Needed + +Consult these documents when the topic is relevant to the user's question: + +### Monovariate Analysis +- **[4_Monovariate_Examples.md](4_Monovariate_Examples.md)** - Detailed walkthroughs of monovariate EVA examples (GEV, GPD, stationary, non-stationary) + +### Multivariate Analysis +- **[5_Copula_Examples.md](5_Copula_Examples.md)** - Detailed walkthroughs of copula-based multivariate analysis +- **[6_Case_Studies.md](6_Case_Studies.md)** - Real-world applications from Bahmanpour et al., 2025 + +### Specialized Topics +- **[7_Visualization_Functions.md](7_Visualization_Functions.md)** - Plotting and visualization functions +- **[8_Data_Preparation.md](8_Data_Preparation.md)** - Data loading, transformation, and preparation utilities +- **[9_Advanced_Topics.md](9_Advanced_Topics.md)** - Ensemble analysis, GOF testing, special configurations + +## Critical Constraints + +### Function Reference Policy +**You may ONLY reference functions documented in [2_Function_Reference.md](2_Function_Reference.md).** + +- Never invent, improvise, or assume functions exist +- If a function is not in the documentation, you MUST NOT suggest it +- If uncertain, verify in [2_Function_Reference.md](2_Function_Reference.md) before suggesting +- This constraint protects users from errors and maintains scientific integrity + +### Copula Family Support +- **Multivariate**: Gaussian, Gumbel +- **Bivariate only**: Frank +- Never suggest Frank for multivariate (>2 variables) analysis + +### Example-Driven Development +- Base ALL code suggestions on documented examples +- Adapt patterns from [4_Monovariate_Examples.md](4_Monovariate_Examples.md) or [5_Copula_Examples.md](5_Copula_Examples.md) +- Cite specific examples when recommending approaches (e.g., "similar to caseStudy01.m") + +## Document Update Procedures + +When tsEVA 2.0 evolves: + +1. Add new functions to [2_Function_Reference.md](2_Function_Reference.md) +2. Add new examples to appropriate example documents +3. Update [3_Workflow_Patterns.md](3_Workflow_Patterns.md) if new patterns emerge +4. **Never remove** documented functions unless officially deprecated + +## Repository Context + +- **GitHub Repository**: https://github.com/menta78/tsEva_dvlp/tree/multivariateArchimedeanCopula +- **Main Branch**: multivariateArchimedeanCopula +- **Language**: MATLAB +- **Dependencies**: None (standalone toolbox) + +## Key References + +These papers define the tsEVA methodology: + +- **Mentaschi et al. (2016)**: Transformed-Stationary approach, *Hydrol. Earth Syst. Sci.*, 20, 3527-3547 +- **Bahmanpour et al. (2025)**: Transformed-Stationary EVA 2.0: A Generalized Framework for Non-Stationary Joint Extremes Analysis (under review) + +## Session Workflow for AI Assistant + +1. **On session start**: Read mandatory files (1-3) +2. **On user question**: Determine topic, consult relevant optional documents (4-9) +3. **On code suggestion**: Verify all functions in [2_Function_Reference.md](2_Function_Reference.md) +4. **On uncertainty**: Check this routing document, then consult appropriate references + +--- + +**Last Updated**: 2026-01-30 +**Documentation Version**: 2.0 +**Maintained By**: tsEVA Development Team diff --git a/docs/2_Function_Reference.md b/docs/2_Function_Reference.md new file mode 100644 index 0000000..e2a887c --- /dev/null +++ b/docs/2_Function_Reference.md @@ -0,0 +1,1331 @@ +# tsEVA 2.0 Function Reference +This document lists ALL documented functions in tsEVA 2.0. **Only functions listed here may be referenced by the AI assistant.** +**Total Functions:** 98 +--- +## Table of Contents +1. [Core Analysis Functions](#core-analysis-functions) + - [Monovariate EVA](#monovariate-eva) (25 functions) + - [Copula Analysis](#copula-analysis) (23 functions) +2. [Plotting & Visualization](#plotting--visualization) (23 functions) +3. [Data Preparation & Utilities](#data-preparation--utilities) (27 functions) +4. [Function Details](#function-details) + +--- +## Quick Summary +| Category | Count | Description | +|----------|-------|-------------| +| Monovariate EVA | 25 | GEV, GPD, stationary, non-stationary analysis | +| Copula Analysis | 23 | Copula fitting, simulation, multivariate analysis | +| Plotting & Visualization | 23 | All plotting and graphical output functions | +| Data Preparation & Utilities | 27 | Data transformation, helper functions, utilities | + +--- +## Core Analysis Functions +### Monovariate EVA +Functions for univariate extreme value analysis using GEV and GPD distributions. + +| Function | Description | +|----------|-------------| +| `Modified_MannKendall_test` | % FUNCTION INPUTS AND OUTPUTS | +| `tsApproxP` | No description available | +| `tsEstimateAverageSeasonality` | estimating the first 2 fourier components | +| `tsEvaComputeAnnualMaxima` | No description available | +| `tsEvaComputeAnnualMaximaMtx` | No description available | +| `tsEvaComputeMonthlyMaxima` | No description available | +| `tsEvaComputeReturnLevelsGEV` | tsEvaComputeReturnLevelsGEV: returns the return levels given the gev parameters | +| `tsEvaComputeReturnLevelsGEVFromAnalysisObj` | No description available | +| `tsEvaComputeReturnLevelsGPD` | reference: Stuart Coles 2001, pag 81. sampleTimeHorizon and returnPeriods must b | +| `tsEvaComputeReturnLevelsGPDFromAnalysisObj` | percentile = nonStationaryEvaParams(2).parameters.percentile; dtSample = nonStat | +| `tsEvaDetrendTimeSeries` | No description available | +| `tsEvaNanRunningMean` | at both extremeties of the series, half windowSize is used which gradually incre | +| `tsEvaNanRunningPercentile` | tsEvaNanRunningPercentile: computes a runnig percentile for a given series, usin | +| `tsEvaNanRunningStatistics` | No description available | +| `tsEvaNanRunningVariance` | !!! series must be 0 averaged!! | +| `tsEvaNonStationary` | tsEvaNonStationary: performs the TS EVA analysis as described by Mentaschi et al | +| `tsEvaReduceOutputObjSize` | redStatTransData.stationarySeries = redStatTransData.stationarySeries(tsIndxs); | +| `tsEvaRunningMeanTrend` | No description available | +| `tsEvaStationary` | tsEvaStationary: executes a regular stationary EVA on timeAndSeries. stationaryE | +| `tsEvaTransformSeriesToStatSeasonal_ciPercentile` | this function decomposes the series into a season-dependent trend and a season-d | +| `tsEvaTransformSeriesToStationaryMultiplicativeSeasonality` | this function decomposes the series into a season-dependent trend and a season-d | +| `tsEvaTransformSeriesToStationaryTrendLinear` | this function first calculates linear trend of the series and also return linear | +| `tsEvaTransformSeriesToStationaryTrendOnly` | further smoothing | +| `tsEvaTransformSeriesToStationaryTrendOnly_ciPercentile` | normalizing to standard deviation (just to be able to make acceptable graphs wit | +| `tsGpdNegShapeFit` | TSGPDNEGSHAPEFIT Fit GPD with negative shape constraint (k < 0) [paramEsts, para | + +### Copula Analysis +Functions for multivariate dependence modeling using copulas. + +| Function | Description | +|----------|-------------| +| `tsCopulaCdfFromSamples` | tsCopulaCdfFromSamples Empirical copula CDF from sample points u : (q x d | +| `tsCopulaComputeBivarRP` | tsCopulaComputeBivarRP computing of bivariate return period of type "AND" | +| `tsCopulaComputeandPlotBivarRP` | tsCopulaComputeandPlotBivarRP computing and plotting of bivariate return period | +| `tsCopulaExtremes` | tsCopulaExtremes joint distribution of non-stationary compound extremes | +| `tsCopulaFit` | Replace negatives in the lower triangle | +| `tsCopulaGOFNonStat` | tsCopulaGOFNonStat estimation of copula goodness-of-fit and other battery of sta | +| `tsCopulaGetFamilyFromId` | No description available | +| `tsCopulaGetFamilyId` | No description available | +| `tsCopulaMontecarlo` | tsCopulaCompoundGPDMontecarlo pefrom Monte-Carlo simulation (resampling) from a | +| `tsCopulaPeakExtrPlotSctrBivar` | No description available | +| `tsCopulaPlotBivariate` | tsCopulaPlotBivariate plotting of joint peaks fitted by a copula | +| `tsCopulaPlotJointReturnPeriod` | plots the multivariate return period according to AND/OR Scenarios see https://d | +| `tsCopulaPlotTrivariate` | tsCopulaPlotTrivariate plotting joint peaks and Monte-Carlo resampled values | +| `tsCopulaPlotTrivariateWithMap` | tsCopulaPlotTrivariate plotting joint peaks and Monte-Carlo resampled values | +| `tsCopulaRnd` | generates the random vector from a copula distribution. if multivariate t or gau | +| `tsCopulaSampleJointPeaksMultiVariatePruning` | tsCopulaSampleJointPeaksMultiVariatePruning multivariate peak-over-thresho | +| `tsCopulaYearExtrDistribution` | normal copula | +| `tsCopulaYearExtrFit` | this is a stationary set of return levels | +| `tsCopulaYearExtrGetMltvrtRetPeriod` | computes the multivariate return period according to Slavadori and De Michele 20 | +| `tsCopulaYearExtrPlotJdistTrivar` | No description available | +| `tsCopulaYearExtrPlotSctrBivar` | No description available | +| `tsCopulaYearExtrPlotSctrTrivar` | No description available | +| `tsCopulaYearExtrRnd` | No description available | + +## Plotting & Visualization +Functions for creating plots and visualizations. + +| Function | Description | +|----------|-------------| +| `plotGEV3D` | No description available | +| `plotGPD3D` | No description available | +| `plotGPD3DFromAnalysisObj` | No description available | +| `plotReturnLevelsGEVStationary` | nonStationaryEvaParams, stationaryTransformData are as these returned by functio | +| `plotReturnLevelsGPDStationary` | nonStationaryEvaParams, stationaryTransformData are as these returned by functio | +| `tsEvaPlotGEV3DFromAnalysisObj` | No description available | +| `tsEvaPlotGEVImageSc` | npdf = (year(maxTS) - year(minTS) + 1)*args.nPlottedTimesByYear; | +| `tsEvaPlotGEVImageScFromAnalysisObj` | No description available | +| `tsEvaPlotGPDImageSc` | npdf = (year(maxTS) - year(minTS) + 1)*args.nPlottedTimesByYear; | +| `tsEvaPlotGPDImageScFromAnalysisObj` | No description available | +| `tsEvaPlotReturnLevelsGEV` | No description available | +| `tsEvaPlotReturnLevelsGEVFromAnalysisObj` | nonStationaryEvaParams, stationaryTransformData are as these returned by functio | +| `tsEvaPlotReturnLevelsGPD` | No description available | +| `tsEvaPlotReturnLevelsGPDFromAnalysisObj` | nonStationaryEvaParams, stationaryTransformData are as these returned by functio | +| `tsEvaPlotSeasonalityGev` | tsEvaPlotSeasonalityGev plots a single year of data adding the series of monthly | +| `tsEvaPlotSeasonalityGevFromAnalysisObj` | No description available | +| `tsEvaPlotSeriesTrendStdDev` | No description available | +| `tsEvaPlotSeriesTrendStdDevFromAnalyisObj` | No description available | +| `tsEvaPlotTransfToStat` | stdThirdMom = third root of the third statistical momentum stdFouthMom = fourth | +| `tsEvaPlotTransfToStatFromAnalysisObj` | std3mom = nthroot(stationaryTransformData.statSer3Mom, 3.); std4mom = nthroot(st | +| `tsLcSubplotManager` | No description available | +| `tsPlotBivarReturnPeriod` | tsCopulaPlotJointReturnPeriod plotting of multivariate return periods | +| `tsPlotSeriesPotGPDRetLevFromAnalysisObj` | computing the return levels | + +## Data Preparation & Utilities +Helper functions for data manipulation and transformation. + +| Function | Description | +|----------|-------------| +| `cvineOrder` | Root-first C-vine order from a Gumbel-θ matrix (alpha). Heuristic: pick node wit | +| `tsEVstatistics` | Evangelos Voukouvalas, Michalis Vousdoukas 2015 gevMaxima can be annual or month | +| `tsEasyParseNamedArgs` | No description available | +| `tsEmpirical` | tsEmpirical Empirical copula [ C ] = tsEmpirical( U ) retuns a variable C contai | +| `tsEnsemble` | calls tsEnsembleEvaParams and tsEnsembleStatTransfData. ASSUMING THAT ALL THE SE | +| `tsEnsembleEvaParams` | From a cell array of non nonStatEvaParams computes the average nonStatEvaParamsA | +| `tsEnsembleStatTransfData` | From a cell array of non stationaryTransformData computes the average stationary | +| `tsEstimateConfidenceIntervalOfRL` | % tsEstimateConfidenceIntervalOfRL: estimates the confidence interval of the ret | +| `tsEvaFillSeries` | ensuring monotonic time vector and no dubplicates. In case of dubplicate time st | +| `tsEvaGetReturnPeriodOfLevelGEV` | GEV | +| `tsEvaGetReturnPeriodOfLevelGPD` | No description available | +| `tsEvaGetTimeStep` | No description available | +| `tsEvaSampleData` | args.pcts = [50 70 85:2:99 99.1:0.1:99.9 99.91:0.01:99.99]; | +| `tsGet` | No description available | +| `tsGetNumberPerYear` | function nperYear=tsGetNumberPerYear(ms,locs) Gives number of events per year fr | +| `tsGetPOT` | function POTdata=tsGetPOT(ms,pcts,desiredEventsPerYear) Gets POT using an automa | +| `tsGetReturnPeriodOfLevel` | given a list of retPeriod and corresponding retLevel with error, estimates the r | +| `tsInterp1Extrap` | function nvals=interp1Extrap(Tr,values,RP,logExtrap) This function interpolates | +| `tsLinearExtrapolation` | Linearly extrapolates time series x,y to points extrax, according to the slope o | +| `tsLoglogExtrapolation` | Extrapolates exponentiallytime series x,y to points extrax, according to the slo | +| `tsPseudoObservations` | No description available | +| `tsRankmax` | tsRankmax Returns vector of one-based ranks for each element | +| `tsRemoveConstantSubseries` | No description available | +| `tsRoundSDate` | function [sdround,dvec,sdunique,dvunique]=tsRoundSDate(sd,sd_precision) ROunds s | +| `tsSameValuesSegmentation` | function [inds,rinds]=tsSameValuesSegmentation(iii,val) separates segments of sa | +| `tsTimeSeriesToPointData` | tsTimeSeriesToPointData: given a ms produces a structure pointData like the one | +| `tsYear` | No description available | + +--- +## Function Details +Comprehensive documentation for each function. + +### Monovariate EVA +#### Modified_MannKendall_test +- **Signature**: `[tau, z, p, H] = Modified_MannKendall_test(t, X, alpha, alpha_ac)` +- **Description**: % FUNCTION INPUTS AND OUTPUTS +- **Category**: Monovariate EVA +- **Inputs**: + - `t` + - `X` + - `alpha` + - `alpha_ac` +- **Outputs**: + - `tau` + - `z` + - `p` + - `H` + +#### tsApproxP +- **Signature**: `[Pval]=tsApproxP(N,copulaFamily,rho,nu,snSample,s2Sample)` +- **Category**: Monovariate EVA +- **Inputs**: + - `N` + - `copulaFamily` + - `rho` + - `nu` + - `snSample` + - `s2Sample` +- **Outputs**: + - `Pval` + +#### tsEstimateAverageSeasonality +- **Signature**: `averageSeasonalitySeries = tsEstimateAverageSeasonality( timeStamps, seasonalitySeries )` +- **Description**: estimating the first 2 fourier components +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `seasonalitySeries` +- **Outputs**: + - `averageSeasonalitySeries` + +#### tsEvaComputeAnnualMaxima +- **Signature**: `[annualMax, annualMaxDate, annualMaxIndx] = tsEvaComputeAnnualMaxima(timeAndSeries)` +- **Category**: Monovariate EVA +- **Inputs**: + - `timeAndSeries` +- **Outputs**: + - `annualMax` + - `annualMaxDate` + - `annualMaxIndx` + +#### tsEvaComputeAnnualMaximaMtx +- **Signature**: `[annualMax, annualMaxDate, annualMaxIndx] = tsEvaComputeAnnualMaximaMtx(timeStamps, srs)` +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `srs` +- **Outputs**: + - `annualMax` + - `annualMaxDate` + - `annualMaxIndx` + +#### tsEvaComputeMonthlyMaxima +- **Signature**: `[monthlyMax, monthlyMaxDate, monthlyMaxIndx] = tsEvaComputeMonthlyMaxima(timeAndSeries)` +- **Category**: Monovariate EVA +- **Inputs**: + - `timeAndSeries` +- **Outputs**: + - `monthlyMax` + - `monthlyMaxDate` + - `monthlyMaxIndx` + +#### tsEvaComputeReturnLevelsGEV +- **Signature**: `[returnLevels, returnLevelsErr] = tsEvaComputeReturnLevelsGEV( epsilon, sigma, mu, epsilonStdErr, sigmaStdErr, muStdErr, returnPeriodsInDts, varargin )` +- **Description**: tsEvaComputeReturnLevelsGEV: returns the return levels given the gev parameters and their standard error. The parameter returnPeriodsInDts contains the return period expressed in a time unit that corresponds to the size of the time segments where we +- **Category**: Monovariate EVA +- **Inputs**: + - `epsilon` + - `sigma` + - `mu` + - `epsilonStdErr` + - `sigmaStdErr` + - `muStdErr` + - `returnPeriodsInDts` + - `varargin` +- **Outputs**: + - `returnLevels` + - `returnLevelsErr` + +#### tsEvaComputeReturnLevelsGEVFromAnalysisObj +- **Signature**: `[returnLevels, returnLevelsErr, returnLevelsErrFit, returnLevelsErrTransf] = tsEvaComputeReturnLevelsGEVFromAnalysisObj(nonStationaryEvaParams, returnPeriodsInYears, varargin)` +- **Category**: Monovariate EVA +- **Inputs**: + - `nonStationaryEvaParams` + - `returnPeriodsInYears` + - `varargin` +- **Outputs**: + - `returnLevels` + - `returnLevelsErr` + - `returnLevelsErrFit` + - `returnLevelsErrTransf` + +#### tsEvaComputeReturnLevelsGPD +- **Signature**: `[returnLevels, returnLevelsErr] = tsEvaComputeReturnLevelsGPD( epsilon, sigma, threshold, epsilonStdErr, sigmaStdErr, thresholdStdErr, nPeaks, sampleTimeHorizon, returnPeriods )` +- **Description**: reference: Stuart Coles 2001, pag 81. sampleTimeHorizon and returnPeriods must be in the same units, e.g. years +- **Category**: Monovariate EVA +- **Inputs**: + - `epsilon` + - `sigma` + - `threshold` + - `epsilonStdErr` + - `sigmaStdErr` + - `thresholdStdErr` + - `nPeaks` + - `sampleTimeHorizon` + - `returnPeriods` +- **Outputs**: + - `returnLevels` + - `returnLevelsErr` + +#### tsEvaComputeReturnLevelsGPDFromAnalysisObj +- **Signature**: `[returnLevels, returnLevelsErr, returnLevelsErrFit, returnLevelsErrTransf] = tsEvaComputeReturnLevelsGPDFromAnalysisObj(nonStationaryEvaParams, returnPeriodsInYears, varargin )` +- **Description**: percentile = nonStationaryEvaParams(2).parameters.percentile; dtSample = nonStationaryEvaParams(2).parameters.timeDeltaYears; +- **Category**: Monovariate EVA +- **Inputs**: + - `nonStationaryEvaParams` + - `returnPeriodsInYears` + - `varargin` +- **Outputs**: + - `returnLevels` + - `returnLevelsErr` + - `returnLevelsErrFit` + - `returnLevelsErrTransf` + +#### tsEvaDetrendTimeSeries +- **Signature**: `[ detrendSeries, trendSeries, filledTimeStamps, filledSeries, nRunMn ] = tsEvaDetrendTimeSeries( timeStamps, series, timeWindow, varargin )` +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `series` + - `timeWindow` + - `varargin` +- **Outputs**: + - `detrendSeries` + - `trendSeries` + - `filledTimeStamps` + - `filledSeries` + - `nRunMn` + +#### tsEvaNanRunningMean +- **Signature**: `[ rnmn ] = tsEvaNanRunningMean( series, windowSize )` +- **Description**: at both extremeties of the series, half windowSize is used which gradually increases to reach windowSize; once windowSize is reached, windowSize is rolled throghout the series +- **Category**: Monovariate EVA +- **Inputs**: + - `series` + - `windowSize` +- **Outputs**: + - `rnmn` + +#### tsEvaNanRunningPercentile +- **Signature**: `[ rnprcnt, stdError ] = tsEvaNanRunningPercentile( series, windowSize, percent, varargin )` +- **Description**: tsEvaNanRunningPercentile: computes a runnig percentile for a given series, using a window with a size given by windowSize. +- **Category**: Monovariate EVA +- **Inputs**: + - `windowSize`: size of the window for the running percentile. Cannot be < 1000 + - `percent`: percent level to which the percentile is compute. + - `percentDelta`: delta for the computation of a percentile interval + - `parameter`: around the requested percentage. If for example + - `parameter`: percent==90 and percentDelta==1, then the 89th, 90th and + - `parameter`: 91st percentiles are computed. Default value: 1 if + - `parameter`: windowSize > 2000, 2 if 2000 > windowsize > 1000. + - `nLowLimit`: minimum number of non nan elements for a window for + - `parameter`: percentile computation. + - `rnprcnt`: approximated running percentile. + - `parameter`: How it works: + - `parameter`: let's suppose that percent == 90. + - `parameter`: For the first window we compute the right percentile using matlab + - `parameter`: function prctile, for percentages 89, 90, 91. + - `parameter`: Then for each step, we update these percentages on the basis + - `parameter`: of the quitting values and incoming values, + - `parameter`: and interpolate an approximated percentile for the requested percentage. +- **Outputs**: + - `rnprcnt` + - `stdError` + +#### tsEvaNanRunningStatistics +- **Signature**: `[ rnmn, rnvar, rn3mom, rn4mom ] = tsEvaNanRunningStatistics( series, windowSize )` +- **Category**: Monovariate EVA +- **Inputs**: + - `series` + - `windowSize` +- **Outputs**: + - `output`: to the forth. + - `rnmn`: running mean + - `rnvar`: running variance + - `rn3mom`: running third statistical momentum + - `rn4mom`: running fourth statistical momentum + +#### tsEvaNanRunningVariance +- **Signature**: `[ rnmn ] = tsEvaNanRunningVariance( series, windowSize )` +- **Description**: !!! series must be 0 averaged!! +- **Category**: Monovariate EVA +- **Inputs**: + - `series` + - `windowSize` +- **Outputs**: + - `rnmn` + +#### tsEvaNonStationary +- **Signature**: `[nonStationaryEvaParams, stationaryTransformData, isValid] = tsEvaNonStationary( timeAndSeries, timeWindow, varargin )` +- **Description**: tsEvaNonStationary: performs the TS EVA analysis as described by Mentaschi et al 2016. +- **Category**: Monovariate EVA +- **Inputs**: + - `timeAndSeries`: array with shape nx2, with the time stamps in the first + - `parameter`: column and the values in the second. + - `timeWindow`: time window for the transformation expressed in days. + - `transfType`: can assume values + - `parameter`: 1) 'trend': long term variability. The trend is computed + - `parameter`: with a running mean, the ci with the running standard deviation. + - `parameter`: 2) 'seasonal': long term + seasonal variability. The trend is computed + - `parameter`: with a running mean, the ci with the running + - `parameter`: standard deviation. + - `parameter`: 3) 'trendCIPercentile': long term variability. The trend is computed + - `parameter`: with a running mean, the ci with the running xx percentile. + - `parameter`: Using this option the argument ciPercentile is + - `parameter`: mandatory. + - `parameter`: 4) 'trendlinear':long term variability. the trend is + - `parameter`: computed with a linear fit, the ci with a linear fit + - `parameter`: of percentile set by the user + - `parameter`: % sample calls + - `parameter`: nonStatEvaParams = tsEvaNonStationary(ms, timeWindow, 'potPercentiles',[95], 'minPeakDistanceInDays', 3) + - `parameter`: samples POT data using a fixed 95 percentile threshold, with peaks at + - `parameter`: a minimum distance of 3 days, looking for a threshold so that we have an average + - `parameter`: of 5 events every year. + - `parameter`: nonStatEvaParams = tsEvaNonStationary(ms, timeWindow, 'minPeakDistanceInDays', 3, 'desiredeventsperyear', 6) + - `parameter`: samples POT data looking for a threshold so that we have an average + - `parameter`: of 6 events every year. + - `parameter`: nonStatEvaParams = tsEvaNonStationary(ms, timeWindow, 'minPeakDistanceInDays', 3, 'trasftype', 'trendCIPercentile', 'ciPercentile', 99) + - `parameter`: for the transformation uses instead of the moving standard deviation, + - `parameter`: the moving 99th percentile. + - `parameter`: % %%%%%%%%%%%%% +- **Outputs**: + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `isValid` + +#### tsEvaReduceOutputObjSize +- **Signature**: `[ redNonStatEvaParams, redStatTransData ] = tsEvaReduceOutputObjSize( nonStationaryEvaParams, stationaryTransformData, newTimeStamps, varargin )` +- **Description**: redStatTransData.stationarySeries = redStatTransData.stationarySeries(tsIndxs); redStatTransData.nonStatSeries = redStatTransData.nonStatSeries(tsIndxs); +- **Category**: Monovariate EVA +- **Inputs**: + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `newTimeStamps` + - `varargin` +- **Outputs**: + - `redNonStatEvaParams` + - `redStatTransData` + +#### tsEvaRunningMeanTrend +- **Signature**: `[ trendSeries, filledTimeStamps, filledSeries, nRunMn ] = tsEvaRunningMeanTrend( timeStamps, series, timeWindow)` +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `series` + - `timeWindow` +- **Outputs**: + - `trendSeries` + - `filledTimeStamps` + - `filledSeries` + - `nRunMn` + +#### tsEvaStationary +- **Signature**: `[stationaryEvaParams, isValid] = tsEvaStationary( timeAndSeries, varargin )` +- **Description**: tsEvaStationary: executes a regular stationary EVA on timeAndSeries. stationaryEvaParams includes the parameters estimated for GEV and GPD % %%%%%%%%%%%%% +- **Category**: Monovariate EVA +- **Inputs**: + - `timeAndSeries` + - `varargin` +- **Outputs**: + - `stationaryEvaParams` + - `isValid` + +#### tsEvaTransformSeriesToStatSeasonal_ciPercentile +- **Signature**: `[trasfData] = tsEvaTransformSeriesToStatSeasonal_ciPercentile( timeStamps, series, timeWindow, percentile, varargin )` +- **Description**: this function decomposes the series into a season-dependent trend and a season-dependent standard deviation. The season-dependent standard deviation is given by a seasonal factor multiplied by a slowly varying standard deviation. +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `series` + - `timeWindow` + - `percentile` + - `varargin` +- **Outputs**: + - `trasfData` + +#### tsEvaTransformSeriesToStationaryMultiplicativeSeasonality +- **Signature**: `[trasfData] = tsEvaTransformSeriesToStationaryMultiplicativeSeasonality( timeStamps, series, timeWindow, varargin )` +- **Description**: this function decomposes the series into a season-dependent trend and a season-dependent standard deviation. The season-dependent standard deviation is given by a seasonal factor multiplied by a slowly varying standard deviation. +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `series` + - `timeWindow` + - `varargin` +- **Outputs**: + - `trasfData` + +#### tsEvaTransformSeriesToStationaryTrendLinear +- **Signature**: `[trasfData] = tsEvaTransformSeriesToStationaryTrendLinear( timeStamps, series, timeWindow, percentile, varargin )` +- **Description**: this function first calculates linear trend of the series and also return linear trend of percetile series (at any particular percentile level set by the user) through a fitting algorithm +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `series` + - `timeWindow` + - `percentile` + - `varargin` +- **Outputs**: + - `trasfData` + +#### tsEvaTransformSeriesToStationaryTrendOnly +- **Signature**: `[trasfData] = tsEvaTransformSeriesToStationaryTrendOnly( timeStamps, series, timeWindow, varargin )` +- **Description**: further smoothing +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `series` + - `timeWindow` + - `varargin` +- **Outputs**: + - `trasfData` + +#### tsEvaTransformSeriesToStationaryTrendOnly_ciPercentile +- **Signature**: `[trasfData] = tsEvaTransformSeriesToStationaryTrendOnly_ciPercentile( timeStamps, series, timeWindow, percentile, varargin )` +- **Description**: normalizing to standard deviation (just to be able to make acceptable graphs with the scripts of this library) +- **Category**: Monovariate EVA +- **Inputs**: + - `timeStamps` + - `series` + - `timeWindow` + - `percentile` + - `varargin` +- **Outputs**: + - `trasfData` + +#### tsGpdNegShapeFit +- **Signature**: `[paramEsts, paramCIs] = tsGpdNegShapeFit(data, alphaCI)` +- **Description**: TSGPDNEGSHAPEFIT Fit GPD with negative shape constraint (k < 0) [paramEsts, paramCIs] = tsGpdNegShapeFit(data, alphaCI) mimics gpfit but constrains the shape parameter to be negative. +- **Category**: Monovariate EVA +- **Inputs**: + - `data` + - `alphaCI` +- **Outputs**: + - `paramEsts` + - `paramCIs` + +### Copula Analysis +#### tsCopulaCdfFromSamples +- **Signature**: `[C, se] = tsCopulaCdfFromSamples(u, Usample)` +- **Description**: tsCopulaCdfFromSamples Empirical copula CDF from sample points u : (q x d) query points in [0,1]^d Usample : (M x d) sample points from the fitted copula C : (q x 1) empirical CDF estimates se : (q x 1) standard errors sqrt(C(1-C)/M) [optional] +- **Category**: Copula Analysis +- **Inputs**: + - `u` + - `Usample` +- **Outputs**: + - `C` + - `se` + +#### tsCopulaComputeBivarRP +- **Signature**: ` [rpAnalysis] = tsCopulaComputeBivarRP(copulaAnalysis, monteCarloAnalysis, varargin)` +- **Description**: tsCopulaComputeBivarRP computing of bivariate return period of type "AND" +- **Category**: Copula Analysis +- **Inputs**: + - `copulaAnalysis` + - `monteCarloAnalysis` + - `varargin` +- **Outputs**: + - `rpAnalysis` + +#### tsCopulaComputeandPlotBivarRP +- **Signature**: ` [rpAnalysis,hFig]=tsCopulaComputeandPlotBivarRP(copulaAnalysis,varargin)` +- **Description**: tsCopulaComputeandPlotBivarRP computing and plotting of bivariate return period of type "AND" +- **Category**: Copula Analysis +- **Inputs**: + - `copulaAnalysis` + - `varargin` +- **Outputs**: + - `rpAnalysis` + - `hFig` + +#### tsCopulaExtremes +- **Signature**: `[CopulaAnalysis] = tsCopulaExtremes(inputtimestamps,inputtimeseries, varargin)` +- **Description**: tsCopulaExtremes joint distribution of non-stationary compound extremes +- **Category**: Copula Analysis +- **Inputs**: + - `inputtimestamps` + - `inputtimeseries` + - `varargin` +- **Outputs**: + - `CopulaAnalysis` + +#### tsCopulaFit +- **Signature**: `copulaParam = tsCopulaFit(copulaFamily, uProb)` +- **Description**: Replace negatives in the lower triangle +- **Category**: Copula Analysis +- **Inputs**: + - `copulaFamily` + - `uProb` +- **Outputs**: + - `copulaParam` + +#### tsCopulaGOFNonStat +- **Signature**: `[gofStatistics] = tsCopulaGOFNonStat(copulaAnalysis, monteCarloAnalysis, varargin)` +- **Description**: tsCopulaGOFNonStat estimation of copula goodness-of-fit and other battery of statistics +- **Category**: Copula Analysis +- **Inputs**: + - `copulaAnalysis` + - `monteCarloAnalysis` + - `varargin` +- **Outputs**: + - `gofStatistics` + +#### tsCopulaGetFamilyFromId +- **Signature**: `copulaFamily = tsCopulaGetFamilyFromId(familyId)` +- **Category**: Copula Analysis +- **Inputs**: + - `familyId` +- **Outputs**: + - `copulaFamily` + +#### tsCopulaGetFamilyId +- **Signature**: `familyId = tsCopulaGetFamilyId(copulaFamily)` +- **Category**: Copula Analysis +- **Inputs**: + - `copulaFamily` +- **Outputs**: + - `familyId` + +#### tsCopulaMontecarlo +- **Signature**: `[monteCarloAnalysis] = tsCopulaMontecarlo(copulaAnalysis, varargin)` +- **Description**: tsCopulaCompoundGPDMontecarlo pefrom Monte-Carlo simulation (resampling) from a pre-determined copula function +- **Category**: Copula Analysis +- **Inputs**: + - `copulaAnalysis` + - `varargin` +- **Outputs**: + - `monteCarloAnalysis` + +#### tsCopulaPeakExtrPlotSctrBivar +- **Signature**: `[handles] = tsCopulaPeakExtrPlotSctrBivar(monteCarloRsmpl, yMaxLevel, varargin)` +- **Category**: Copula Analysis +- **Inputs**: + - `monteCarloRsmpl` + - `yMaxLevel` + - `varargin` +- **Outputs**: + - `handles` + +#### tsCopulaPlotBivariate +- **Signature**: `[axxArray] = tsCopulaPlotBivariate(copulaAnalysis, monteCarloAnalysis, varargin)` +- **Description**: tsCopulaPlotBivariate plotting of joint peaks fitted by a copula +- **Category**: Copula Analysis +- **Inputs**: + - `copulaAnalysis` + - `monteCarloAnalysis` + - `varargin` +- **Outputs**: + - `axxArray` + +#### tsCopulaPlotJointReturnPeriod +- **Signature**: ` tsCopulaPlotJointReturnPeriod(copulaAnalysis,varargin)` +- **Description**: plots the multivariate return period according to AND/OR Scenarios see https://doi.org/10.1002/2015WR017225 Parts of the code were reworked from MvCAT toolbox https://doi.org/10.1002/2016WR020242 Bahmanpour, M.H., 2023 +- **Category**: Copula Analysis +- **Inputs**: + - `copulaAnalysis` + - `varargin` + +#### tsCopulaPlotTrivariate +- **Signature**: `[axxArray] = tsCopulaPlotTrivariate(copulaAnalysis, monteCarloAnalysis, varargin)` +- **Description**: tsCopulaPlotTrivariate plotting joint peaks and Monte-Carlo resampled values +- **Category**: Copula Analysis +- **Inputs**: + - `copulaAnalysis` + - `monteCarloAnalysis` + - `varargin` +- **Outputs**: + - `axxArray` + +#### tsCopulaPlotTrivariateWithMap +- **Signature**: `[axxArray] = tsCopulaPlotTrivariateWithMap(copulaAnalysis,gofStatistics,varargin)` +- **Description**: tsCopulaPlotTrivariate plotting joint peaks and Monte-Carlo resampled values +- **Category**: Copula Analysis +- **Inputs**: + - `copulaAnalysis` + - `gofStatistics` + - `varargin` +- **Outputs**: + - `axxArray` + +#### tsCopulaRnd +- **Signature**: `u = tsCopulaRnd(family, copulaPar, N, uProb)` +- **Description**: generates the random vector from a copula distribution. if multivariate t or gaussian or bivariate archimedean, it uses copularnd else (this would be the case of multivariate archimedean) it implements a simple C-Vine copula +- **Category**: Copula Analysis +- **Inputs**: + - `family`: copula family + - `copulaPar`: rho for gaussian, alpha for archimedean + - `N`: sample size + - `uProb`: pseudo-observation data (that is, transformed in uniform probability space) + - `parameter`: from which the copula was built. Necessary in case the + - `c`: vine copula must be constructed +- **Outputs**: + - `u` + +#### tsCopulaSampleJointPeaksMultiVariatePruning +- **Signature**: `[samplingAnalysis] = tsCopulaSampleJointPeaksMultiVariatePruning(inputtimestamps,inputtimeseries,varargin)` +- **Description**: tsCopulaSampleJointPeaksMultiVariatePruning multivariate peak-over-threshold sampling of compound events +- **Category**: Copula Analysis +- **Inputs**: + - `inputtimestamps` + - `inputtimeseries` + - `varargin` +- **Outputs**: + - `samplingAnalysis` + +#### tsCopulaYearExtrDistribution +- **Signature**: `[jpdf, jcdf] = tsCopulaYearExtrDistribution(retPeriod, copulaParam, varargin)` +- **Description**: normal copula +- **Category**: Copula Analysis +- **Inputs**: + - `retPeriod` + - `copulaParam` + - `varargin` +- **Outputs**: + - `jpdf` + - `jcdf` + +#### tsCopulaYearExtrFit +- **Signature**: `[retLev, copulaParam, yRetPer, yProb] = tsCopulaYearExtrFit(retPeriod, retLev, yMax, varargin)` +- **Description**: this is a stationary set of return levels +- **Category**: Copula Analysis +- **Inputs**: + - `retPeriod` + - `retLev` + - `yMax` + - `varargin` +- **Outputs**: + - `retLev` + - `copulaParam` + - `yRetPer` + - `yProb` + +#### tsCopulaYearExtrGetMltvrtRetPeriod +- **Signature**: `[returnPeriod, prob] = tsCopulaYearExtrGetMltvrtRetPeriod(randomSample, level)` +- **Description**: computes the multivariate return period according to Slavadori and De Michele 2004, Salvadori et al. 2011, used by Zscheischler et al. 2017 +- **Category**: Copula Analysis +- **Inputs**: + - `randomSample` + - `level` +- **Outputs**: + - `returnPeriod` + - `prob` + +#### tsCopulaYearExtrPlotJdistTrivar +- **Signature**: `handles = tsCopulaYearExtrPlotJdistTrivar( retLev, jdist, varargin )` +- **Category**: Copula Analysis +- **Inputs**: + - `retLev` + - `jdist` + - `varargin` +- **Outputs**: + - `handles` + +#### tsCopulaYearExtrPlotSctrBivar +- **Signature**: `[handles] = tsCopulaYearExtrPlotSctrBivar(monteCarloRsmpl, yMaxLevel, varargin)` +- **Category**: Copula Analysis +- **Inputs**: + - `monteCarloRsmpl` + - `yMaxLevel` + - `varargin` +- **Outputs**: + - `handles` + +#### tsCopulaYearExtrPlotSctrTrivar +- **Signature**: `handles = tsCopulaYearExtrPlotSctrTrivar(monteCarloRsmpl, yMaxLevel, varargin)` +- **Category**: Copula Analysis +- **Inputs**: + - `monteCarloRsmpl` + - `yMaxLevel` + - `varargin` +- **Outputs**: + - `handles` + +#### tsCopulaYearExtrRnd +- **Signature**: `[monteCarloRsmpl, resampleProb, resampleRetPer] = tsCopulaYearExtrRnd(retPeriod, retLev, copulaParam, nResample, varargin)` +- **Category**: Copula Analysis +- **Inputs**: + - `retPeriod` + - `retLev` + - `copulaParam` + - `nResample` + - `varargin` +- **Outputs**: + - `monteCarloRsmpl` + - `resampleProb` + - `resampleRetPer` + +### Plotting & Visualization +#### plotGEV3D +- **Signature**: `phandles = plotGEV3D( X, timeStamps, epsilon, sigma, mu, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `X` + - `timeStamps` + - `epsilon` + - `sigma` + - `mu` + - `varargin` +- **Outputs**: + - `phandles` + +#### plotGPD3D +- **Signature**: `phandles = plotGPD3D( X, timeStamps, epsilon, sigma, threshold, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `X` + - `timeStamps` + - `epsilon` + - `sigma` + - `threshold` + - `varargin` +- **Outputs**: + - `phandles` + +#### plotGPD3DFromAnalysisObj +- **Signature**: `phandles = plotGPD3DFromAnalysisObj( X, nonStationaryEvaParams, stationaryTransformData )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `X` + - `nonStationaryEvaParams` + - `stationaryTransformData` +- **Outputs**: + - `phandles` + +#### plotReturnLevelsGEVStationary +- **Signature**: `phandles = plotReturnLevelsGEVStationary( nonStationaryEvaParams, varargin )` +- **Description**: nonStationaryEvaParams, stationaryTransformData are as these returned by function nonStationaryEvaJRCApproach timeIndex is the index at which the time varying analysis should be estimated. +- **Category**: Plotting & Visualization +- **Inputs**: + - `nonStationaryEvaParams` + - `varargin` +- **Outputs**: + - `phandles` + +#### plotReturnLevelsGPDStationary +- **Signature**: `phandles = plotReturnLevelsGPDStationary( nonStationaryEvaParams, varargin )` +- **Description**: nonStationaryEvaParams, stationaryTransformData are as these returned by function nonStationaryEvaJRCApproach timeIndex is the index at which the time varying analysis should be estimated. +- **Category**: Plotting & Visualization +- **Inputs**: + - `nonStationaryEvaParams` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotGEV3DFromAnalysisObj +- **Signature**: `phandles = tsEvaPlotGEV3DFromAnalysisObj( X, nonStationaryEvaParams, stationaryTransformData, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `X` + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotGEVImageSc +- **Signature**: `phandles = tsEvaPlotGEVImageSc( Y, timeStamps, epsilon, sigma, mu, varargin )` +- **Description**: npdf = (year(maxTS) - year(minTS) + 1)*args.nPlottedTimesByYear; +- **Category**: Plotting & Visualization +- **Inputs**: + - `Y` + - `timeStamps` + - `epsilon` + - `sigma` + - `mu` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotGEVImageScFromAnalysisObj +- **Signature**: `phandles = tsEvaPlotGEVImageScFromAnalysisObj( X, nonStationaryEvaParams, stationaryTransformData, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `X` + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotGPDImageSc +- **Signature**: `phandles = tsEvaPlotGPDImageSc( Y, timeStamps, epsilon, sigma, threshold, varargin )` +- **Description**: npdf = (year(maxTS) - year(minTS) + 1)*args.nPlottedTimesByYear; +- **Category**: Plotting & Visualization +- **Inputs**: + - `Y` + - `timeStamps` + - `epsilon` + - `sigma` + - `threshold` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotGPDImageScFromAnalysisObj +- **Signature**: `phandles = tsEvaPlotGPDImageScFromAnalysisObj( Y, nonStationaryEvaParams, stationaryTransformData, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `Y` + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotReturnLevelsGEV +- **Signature**: `phandles = tsEvaPlotReturnLevelsGEV( epsilon, sigma, mu, epsilonStdErr, sigmaStdErr, muStdErr, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `epsilon` + - `sigma` + - `mu` + - `epsilonStdErr` + - `sigmaStdErr` + - `muStdErr` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotReturnLevelsGEVFromAnalysisObj +- **Signature**: `phandles = tsEvaPlotReturnLevelsGEVFromAnalysisObj( nonStationaryEvaParams, timeIndex, varargin )` +- **Description**: nonStationaryEvaParams, stationaryTransformData are as these returned by function nonStationaryEvaJRCApproach timeIndex is the index at which the time varying analysis should be estimated. +- **Category**: Plotting & Visualization +- **Inputs**: + - `nonStationaryEvaParams` + - `timeIndex` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotReturnLevelsGPD +- **Signature**: `[phandles, returnPeriods, returnLevels, retrunLevelsErrs] = tsEvaPlotReturnLevelsGPD( epsilon, sigma, threshold, epsilonStdErr, sigmaStdErr, thresholdStdErr, nPeaks, timeHorizonInYears, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `epsilon` + - `sigma` + - `threshold` + - `epsilonStdErr` + - `sigmaStdErr` + - `thresholdStdErr` + - `nPeaks` + - `timeHorizonInYears` + - `varargin` +- **Outputs**: + - `phandles` + - `returnPeriods` + - `returnLevels` + - `retrunLevelsErrs` + +#### tsEvaPlotReturnLevelsGPDFromAnalysisObj +- **Signature**: `[phandles, returnPeriods, returnLevels, retrunLevelsErrs] = tsEvaPlotReturnLevelsGPDFromAnalysisObj( nonStationaryEvaParams, timeIndex, varargin )` +- **Description**: nonStationaryEvaParams, stationaryTransformData are as these returned by function nonStationaryEvaJRCApproach timeIndex is the index at which the time varying analysis should be estimated. +- **Category**: Plotting & Visualization +- **Inputs**: + - `nonStationaryEvaParams` + - `timeIndex` + - `varargin` +- **Outputs**: + - `phandles` + - `returnPeriods` + - `returnLevels` + - `retrunLevelsErrs` + +#### tsEvaPlotSeasonalityGev +- **Signature**: `phandles = tsEvaPlotSeasonalityGev(extremesRange, referenceYear, timeStamps, epsilon, sigma, mu, monthlyMaxIndexes, series, trend, stddev, varargin)` +- **Description**: tsEvaPlotSeasonalityGev plots a single year of data adding the series of monthly maxima renormalized to the considered year and superimposed to the series. I furtherly plots the time varying location and scale parameters. +- **Category**: Plotting & Visualization +- **Inputs**: + - `extremesRange` + - `referenceYear` + - `timeStamps` + - `epsilon` + - `sigma` + - `mu` + - `monthlyMaxIndexes` + - `series` + - `trend` + - `stddev` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotSeasonalityGevFromAnalysisObj +- **Signature**: `phandles = tsEvaPlotSeasonalityGevFromAnalysisObj( extremesRange, referenceYear, nonStationaryEvaParams, stationaryTransformData, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `extremesRange` + - `referenceYear` + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotSeriesTrendStdDev +- **Signature**: `phandles = tsEvaPlotSeriesTrendStdDev( timeStamps, series, trend, stdDev, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `timeStamps` + - `series` + - `trend` + - `stdDev` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotSeriesTrendStdDevFromAnalyisObj +- **Signature**: `phandles = tsEvaPlotSeriesTrendStdDevFromAnalyisObj( nonStationaryEvaParams, stationaryTransformData, varargin )` +- **Category**: Plotting & Visualization +- **Inputs**: + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotTransfToStat +- **Signature**: `phandles = tsEvaPlotTransfToStat( timeStamps, statSeries, srsmean, stdDev, thirdMom, fourthMom, varargin )` +- **Description**: stdThirdMom = third root of the third statistical momentum stdFouthMom = fourth root of the fourth statistical momentum +- **Category**: Plotting & Visualization +- **Inputs**: + - `timeStamps` + - `statSeries` + - `srsmean` + - `stdDev` + - `thirdMom` + - `fourthMom` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsEvaPlotTransfToStatFromAnalysisObj +- **Signature**: `phandles = tsEvaPlotTransfToStatFromAnalysisObj( nonStationaryEvaParams, stationaryTransformData, varargin )` +- **Description**: std3mom = nthroot(stationaryTransformData.statSer3Mom, 3.); std4mom = nthroot(stationaryTransformData.statSer4Mom, 4.); +- **Category**: Plotting & Visualization +- **Inputs**: + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `varargin` +- **Outputs**: + - `phandles` + +#### tsLcSubplotManager +- **Signature**: `obj = tsLcSubplotManager(N, M, varargin)` +- **Category**: Plotting & Visualization +- **Inputs**: + - `N` + - `M` + - `varargin` +- **Outputs**: + - `obj` + +#### tsPlotBivarReturnPeriod +- **Signature**: ` [copulaAnalysis,hFig1] = tsPlotBivarReturnPeriod(copulaAnalysis,axxArray,varargin)` +- **Description**: tsCopulaPlotJointReturnPeriod plotting of multivariate return periods +- **Category**: Plotting & Visualization +- **Inputs**: + - `copulaAnalysis` + - `axxArray` + - `varargin` +- **Outputs**: + - `copulaAnalysis` + - `hFig1` + +#### tsPlotSeriesPotGPDRetLevFromAnalysisObj +- **Signature**: `phandles = tsPlotSeriesPotGPDRetLevFromAnalysisObj( nonStationaryEvaParams, stationaryTransformData, varargin )` +- **Description**: computing the return levels +- **Category**: Plotting & Visualization +- **Inputs**: + - `nonStationaryEvaParams` + - `stationaryTransformData` + - `varargin` +- **Outputs**: + - `phandles` + +### Data Preparation & Utilities +#### cvineOrder +- **Signature**: `order = cvineOrder(alpha)` +- **Description**: Root-first C-vine order from a Gumbel-θ matrix (alpha). Heuristic: pick node with max total |tau| as root, then greedily add the next. +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `alpha` +- **Outputs**: + - `order` + +#### tsEVstatistics +- **Signature**: `[EVmeta,EVdata,isValid] = tsEVstatistics(pointData, varargin)` +- **Description**: Evangelos Voukouvalas, Michalis Vousdoukas 2015 gevMaxima can be annual or monthly. annual by default +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `pointData` + - `varargin` +- **Outputs**: + - `EVmeta` + - `EVdata` + - `isValid` + +#### tsEasyParseNamedArgs +- **Signature**: `[argStruct] = tsEasyParseNamedArgs(args, argStruct)` +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `args` + - `argStruct` +- **Outputs**: + - `argStruct` + +#### tsEmpirical +- **Signature**: `[ C ] = tsEmpirical( U )` +- **Description**: tsEmpirical Empirical copula [ C ] = tsEmpirical( U ) retuns a variable C containing the empirical copula corresponding with the U variable including the uniform variates +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `U` +- **Outputs**: + - `C` + +#### tsEnsemble +- **Signature**: `[nonStationaryEvaParamsEns, stationaryTransformDataEns] = tsEnsemble( nonStatEvaParamsArray, stationaryTransformDataArray )` +- **Description**: calls tsEnsembleEvaParams and tsEnsembleStatTransfData. ASSUMING THAT ALL THE SERIES HAVE THE SAME LENGTH AND UNIFORM TIME STAMPS +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `nonStatEvaParamsArray` + - `stationaryTransformDataArray` +- **Outputs**: + - `nonStationaryEvaParamsEns` + - `stationaryTransformDataEns` + +#### tsEnsembleEvaParams +- **Signature**: `nonStatEvaParamsEnsemble = tsEnsembleEvaParams( nonStatEvaParamsArray )` +- **Description**: From a cell array of non nonStatEvaParams computes the average nonStatEvaParamsArray is a cell array of nonStatEvaParams, the object type returned by tsEvaNonStationary ASSUMING THAT ALL THE SERIES HAVE THE SAME LENGTH AND UNIFORM TIME STAMPS +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `nonStatEvaParamsArray` +- **Outputs**: + - `nonStatEvaParamsEnsemble` + +#### tsEnsembleStatTransfData +- **Signature**: `[ stationaryTransformDataEnsemble ] = tsEnsembleStatTransfData( stationaryTransformDataArray )` +- **Description**: From a cell array of non stationaryTransformData computes the average stationaryTransformDataArray is a cell array of stationaryTransformData, the transformation object type returned by tsEvaNonStationary ASSUMING THAT ALL THE SERIES HAVE THE SAME LENGTH AND UNIFORM TIME STAMPS +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `stationaryTransformDataArray` +- **Outputs**: + - `stationaryTransformDataEnsemble` + +#### tsEstimateConfidenceIntervalOfRL +- **Signature**: `[ lowCI, highCI ] = tsEstimateConfidenceIntervalOfRL( rl, stdErr, p )` +- **Description**: % tsEstimateConfidenceIntervalOfRL: estimates the confidence interval of the return level (rl), assuming that the uncertainty distribution around rl is a lognormal with standard deviation stdErr % +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `rl` + - `stdErr` + - `p` +- **Outputs**: + - `lowCI` + - `highCI` + +#### tsEvaFillSeries +- **Signature**: `[ filledTimeStamps, filledSeries, dt ] = tsEvaFillSeries( timeStamps, series )` +- **Description**: ensuring monotonic time vector and no dubplicates. In case of dubplicate time stamps the highest value is considered +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `timeStamps` + - `series` +- **Outputs**: + - `filledTimeStamps` + - `filledSeries` + - `dt` + +#### tsEvaGetReturnPeriodOfLevelGEV +- **Signature**: `[retPer, exceedProb] = tsEvaGetReturnPeriodOfLevelGEV( epsilon, sigma, mu, retLev)` +- **Description**: GEV +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `epsilon` + - `sigma` + - `mu` + - `retLev` +- **Outputs**: + - `retPer` + - `exceedProb` + +#### tsEvaGetReturnPeriodOfLevelGPD +- **Signature**: `[retPer, exceedProb] = tsEvaGetReturnPeriodOfLevelGPD( epsilon, sigma, pPeak, retLev)` +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `epsilon` + - `sigma` + - `pPeak` + - `retLev` +- **Outputs**: + - `retPer` + - `exceedProb` + +#### tsEvaGetTimeStep +- **Signature**: `dt = tsEvaGetTimeStep( times )` +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `times` +- **Outputs**: + - `dt` + +#### tsEvaSampleData +- **Signature**: `pointData = tsEvaSampleData(ms, varargin)` +- **Description**: args.pcts = [50 70 85:2:99 99.1:0.1:99.9 99.91:0.01:99.99]; +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `ms` + - `varargin` +- **Outputs**: + - `pointData` + +#### tsGet +- **Signature**: `rho = tsGet()` +- **Category**: Data Preparation & Utilities +- **Outputs**: + - `rho` + +#### tsGetNumberPerYear +- **Signature**: `nperYear=tsGetNumberPerYear(ms,locs)` +- **Description**: function nperYear=tsGetNumberPerYear(ms,locs) Gives number of events per year from a time series ms (sd and values vertical vectors) and sporadic indices of events (locs) Michalis Vousdoukas 2015 +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `ms` + - `locs` +- **Outputs**: + - `nperYear` + +#### tsGetPOT +- **Signature**: `[POTdata]=tsGetPOT(ms,pcts,desiredEventsPerYear, varargin)` +- **Description**: function POTdata=tsGetPOT(ms,pcts,desiredEventsPerYear) Gets POT using an automatic threshold such that the mean number of events per year is equal to desiredEventsPerYear +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `parameter`: ms data in two columns sdate and values. IT IS ASSUMED + - `parameter`: pcts vestor of percentiles tested + - `parameter`: desiredEventsPerYear mean number of events per year + - `parameter`: Michalis Vousdoukas, Evangelos Voukouvalas, Lorenzo Mentaschi 2015 +- **Outputs**: + - `POTdata` + +#### tsGetReturnPeriodOfLevel +- **Signature**: `[ myRetPeriod, myRetPeriodCISup, myRetPeriodCIInf ] = tsGetReturnPeriodOfLevel( retPeriod, retLevel, retLevError, myLevel, varargin )` +- **Description**: given a list of retPeriod and corresponding retLevel with error, estimates the return period for a level myLevel, and the related confidence interval. +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `retPeriod` + - `retLevel` + - `retLevError` + - `myLevel` + - `varargin` +- **Outputs**: + - `myRetPeriod` + - `myRetPeriodCISup` + - `myRetPeriodCIInf` + +#### tsInterp1Extrap +- **Signature**: `nvals = tsInterp1Extrap(X,V,Xq,logExtrap)` +- **Description**: function nvals=interp1Extrap(Tr,values,RP,logExtrap) This function interpolates a series and extrapolates linearly or exponeniantly if necessary. Can also handle matrices (only for dependent variable values) Will also sort the Tr series if not monotonically increasing Tr- dependent variable +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `X` + - `V` + - `Xq` + - `logExtrap` +- **Outputs**: + - `nvals` + +#### tsLinearExtrapolation +- **Signature**: `[nx,ny,addedy]=tsLinearExtrapolation(x,y,extrax,npoints)` +- **Description**: Linearly extrapolates time series x,y to points extrax, according to the slope of the first/last number of points of the series (npoints) Michalis Vousdoukas 2016 +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `x` + - `y` + - `extrax` + - `npoints` +- **Outputs**: + - `nx` + - `ny` + - `addedy` + +#### tsLoglogExtrapolation +- **Signature**: `[nx,ny,addedy]=tsLoglogExtrapolation(x,y,extrax,npoints)` +- **Description**: Extrapolates exponentiallytime series x,y to points extrax, according to the slope of the first/last number of points of the series (npoints) Michalis Vousdoukas 2016 +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `x` + - `y` + - `extrax` + - `npoints` +- **Outputs**: + - `nx` + - `ny` + - `addedy` + +#### tsPseudoObservations +- **Signature**: `[ U ] = tsPseudoObservations( X )` +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `X` +- **Outputs**: + - `U` + +#### tsRankmax +- **Signature**: `[ R ] = tsRankmax( X )` +- **Description**: tsRankmax Returns vector of one-based ranks for each element +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `X` +- **Outputs**: + - `R` + +#### tsRemoveConstantSubseries +- **Signature**: `cleaned_series = tsRemoveConstantSubseries( srs, stackedValuesCount )` +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `srs` + - `stackedValuesCount` +- **Outputs**: + - `cleaned_series` + +#### tsRoundSDate +- **Signature**: `[sdround,dvecm,sdunique,dvunique]=tsRoundSDate(sd,sd_precision)` +- **Description**: function [sdround,dvec,sdunique,dvunique]=tsRoundSDate(sd,sd_precision) ROunds sd values according to desired precision, years, months, days, etc : 1-4 Michalis Vousdoukas 2015 +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `sd` + - `sd_precision` +- **Outputs**: + - `sdround` + - `dvecm` + - `sdunique` + - `dvunique` + +#### tsSameValuesSegmentation +- **Signature**: `[inds,rinds]=tsSameValuesSegmentation(iii,varargin)` +- **Description**: function [inds,rinds]=tsSameValuesSegmentation(iii,val) separates segments of same value val +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `iii`: the series + - `parameter`: optional variable-the value +- **Outputs**: + - `inds`: cell with the values of each continuous set + - `rinds`: cell with indexes of each continuous set + - `output`: Michalis Vousdoukas 2009 + +#### tsTimeSeriesToPointData +- **Signature**: `pointData = tsTimeSeriesToPointData( ms, potThreshold, potThresholdError )` +- **Description**: tsTimeSeriesToPointData: given a ms produces a structure pointData like the one produced by tsEvaSampleData +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `ms` + - `potThreshold` + - `potThresholdError` +- **Outputs**: + - `pointData` + +#### tsYear +- **Signature**: `year = tsYear(timeStamp)` +- **Category**: Data Preparation & Utilities +- **Inputs**: + - `timeStamp` +- **Outputs**: + - `year` + From 9f21a22897adf5f9a79503fe9f7ba5a7fa373f30 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 30 Jan 2026 16:58:04 +0000 Subject: [PATCH 3/4] Add core documentation: routing, methodology, functions, and workflows Co-authored-by: menta78 <19606593+menta78@users.noreply.github.com> --- docs/1_Core_Methodology.md | 269 ++++++++++++++++++++++ docs/3_Workflow_Patterns.md | 444 ++++++++++++++++++++++++++++++++++++ 2 files changed, 713 insertions(+) create mode 100644 docs/1_Core_Methodology.md create mode 100644 docs/3_Workflow_Patterns.md diff --git a/docs/1_Core_Methodology.md b/docs/1_Core_Methodology.md new file mode 100644 index 0000000..012adae --- /dev/null +++ b/docs/1_Core_Methodology.md @@ -0,0 +1,269 @@ +# tsEVA 2.0 Core Methodology + +## Overview + +tsEVA 2.0 implements the **Transformed-Stationary (TS) approach** for non-stationary extreme value analysis (EVA). This methodology enables rigorous statistical analysis of extremes even under changing climate conditions. + +## The Transformed-Stationary Paradigm + +### Fundamental Concept + +The TS approach is NOT a detrending technique—it is a fully-fledged non-stationary methodology. Mentaschi et al. (2016) demonstrated that **for each time-varying extreme value distribution, there exists a family of distributions that transform a non-stationary signal into a supposed-stationary one**. + +### Three-Step Process + +1. **Transform**: Non-stationary time series → Stationary time series +2. **Analyze**: Apply stationary EVA theory (GEV, GPD) or copula modeling +3. **Reverse-Transform**: Stationary results → Non-stationary extreme value distribution + +This framework allows us to leverage the well-established theory of stationary extremes while properly accounting for non-stationarity in the data. + +## Key Distributions + +### Generalized Extreme Value (GEV) Distribution + +Used for analyzing block maxima (e.g., annual maxima). + +**Parameters:** +- **Location (μ)**: Center of the distribution +- **Scale (σ)**: Spread of the distribution +- **Shape (ξ)**: Tail behavior + - ξ > 0: Heavy-tailed (Fréchet family) + - ξ = 0: Exponential tail (Gumbel family) + - ξ < 0: Bounded tail (Weibull family) + +### Generalized Pareto Distribution (GPD) + +Used for analyzing peaks over threshold (POT). + +**Parameters:** +- **Scale (σ)**: Spread of exceedances +- **Shape (ξ)**: Tail behavior (same interpretation as GEV) +- **Threshold (u)**: The threshold above which peaks are selected + +## Non-Stationarity Detection + +### Time Window Approach + +tsEVA uses a **moving window** to detect and model non-stationarity: + +- **Time Window**: The minimum period over which statistics are considered stationary +- **Typical Values**: + - Short-term variability: 5-10 years + - Long-term climate trends: 30-50 years +- **Trade-off**: Smaller windows capture more variability but require more data + +### Transformation Types + +tsEVA supports multiple transformation approaches: + +1. **Trend Only** (`'transfType', 'trend'`) + - Models time-varying mean and standard deviation + - No seasonal component + - Best for: Data without strong seasonal cycles + +2. **Seasonal** (`'transfType', 'seasonal'`) + - Models both trend and seasonality + - Captures annual cycles in extremes + - Best for: Temperature, precipitation, coastal data with seasonal patterns + +3. **Percentile-based CI** (`'ciPercentile'`) + - Uses moving percentile instead of moving standard deviation + - More sensitive to changes in extremes + - Trade-off: Broader confidence intervals + +## Multivariate Analysis: Copula Framework + +### Purpose + +Copulas model the **dependence structure** between multiple variables while preserving their individual (marginal) distributions. + +### Supported Copula Families + +#### Gaussian Copula +- **Support**: Multivariate (2+ variables) +- **Dependence**: Symmetric, full correlation matrix +- **Best for**: Variables with symmetric, linear-like dependence + +#### Gumbel Copula +- **Support**: Multivariate (2+ variables) +- **Dependence**: Upper tail dependence (positive association) +- **Best for**: Variables that tend to be extreme together (e.g., storm surge + wave height) + +#### Frank Copula +- **Support**: Bivariate ONLY (exactly 2 variables) +- **Dependence**: Symmetric, no tail dependence +- **Best for**: Variables with symmetric but not necessarily linear dependence +- **Critical**: Never suggest Frank for multivariate (>2 variables) analysis + +### Joint Return Periods + +For multivariate extremes, tsEVA computes: + +- **AND Return Period**: Time until ALL variables exceed thresholds simultaneously +- **OR Return Period**: Time until ANY variable exceeds its threshold +- **Conditional Return Periods**: Given one variable, what's the return period of others + +## Return Levels and Return Periods + +### Definitions + +- **Return Level (RL)**: The level expected to be exceeded on average once every T time units +- **Return Period (RP)**: The average time interval between exceedances of a given level +- **Confidence Intervals**: Quantify uncertainty in return level estimates + +### Typical Applications + +- **Infrastructure Design**: 100-year, 500-year return levels +- **Risk Assessment**: Probability of exceedance over project lifetime +- **Climate Projections**: Future changes in extreme event frequency + +## Data Requirements + +### Monovariate Analysis + +**Minimum Requirements:** +- Sufficient data to populate the time window (e.g., 6-10 years of data for 6-year window) +- Regular or irregular time series (tsEVA handles both) +- Quality-controlled data (gaps are acceptable but should be documented) + +**Optimal Data:** +- Multiple decades for detecting long-term trends +- High temporal resolution (daily to sub-daily) +- Metadata on data collection and quality flags + +### Multivariate Analysis + +**Additional Requirements:** +- **Simultaneous observations** across all variables +- **Sufficient joint extremes**: Typically 50+ joint peak events +- **Peak separation**: Define minimum separation between independent events +- **Temporal alignment**: Define maximum lag for joint occurrence + +## Goodness-of-Fit Assessment + +### Visual Diagnostics + +- **Probability plots**: Compare empirical vs. theoretical distributions +- **Return level plots**: Assess fit across return periods +- **Transformation plots**: Verify stationarity after transformation + +### Statistical Tests + +- **Kolmogorov-Smirnov**: Test distributional fit +- **Anderson-Darling**: Emphasizes tail fit +- **Chi-square**: Goodness-of-fit for copulas + +## Uncertainty Quantification + +### Sources of Uncertainty + +1. **Parameter Uncertainty**: Limited sample size +2. **Model Uncertainty**: Choice of distribution family +3. **Climate Variability**: Natural variability vs. forced trends +4. **Threshold Selection**: For GPD analysis + +### Methods + +- **Bootstrap resampling**: Empirical confidence intervals +- **Profile likelihood**: Asymptotic confidence intervals +- **Ensemble approaches**: Multiple time windows, thresholds + +## Best Practices + +### Peak Selection + +- **Independence**: Ensure peaks are independent events + - Typical separation: 2-5 days for sub-daily data, 1-2 months for monthly maxima +- **Threshold Selection**: Balance bias (too low) vs. variance (too high) + - Rule of thumb: Select threshold giving 3-5 events per year + +### Non-Stationarity Assessment + +1. **Visual inspection**: Plot time series with running statistics +2. **Trend tests**: Mann-Kendall, modified Mann-Kendall +3. **Compare stationary vs. non-stationary**: Use likelihood ratio tests + +### Model Selection + +- **Start simple**: Try stationary model first +- **Add complexity**: Only if non-stationarity is clear +- **Physical plausibility**: Ensure trends align with known processes +- **Parsimony**: Avoid over-fitting with too many parameters + +## Common Pitfalls + +### Over-fitting Non-Stationarity + +**Problem**: Fitting trends to natural variability +**Solution**: Use sufficiently long time windows, validate with independent data + +### Under-estimating Uncertainty + +**Problem**: Ignoring model and parameter uncertainty +**Solution**: Report confidence intervals, consider multiple models + +### Inappropriate Copula Choice + +**Problem**: Frank copula for multivariate analysis +**Solution**: Use Gaussian or Gumbel for >2 variables + +### Ignoring Physical Context + +**Problem**: Statistical fit without physical interpretation +**Solution**: Ground analysis in climate science, validate against known processes + +## Theoretical Foundations + +### Key Papers + +1. **Mentaschi et al. (2016)**: Original TS approach for monovariate analysis + - *Hydrology and Earth System Sciences*, 20, 3527-3547 + - DOI: 10.5194/hess-20-3527-2016 + +2. **Bahmanpour et al. (2025)**: Extension to multivariate analysis with copulas + - *Hydrology and Earth System Sciences* (under review) + - Introduces time-varying copula framework + +### Extreme Value Theory Background + +- **Fisher-Tippett-Gnedenko theorem**: Justifies GEV for block maxima +- **Pickands-Balkema-de Haan theorem**: Justifies GPD for threshold exceedances +- **Sklar's theorem**: Justifies copula decomposition for multivariate distributions + +## Relation to Other Approaches + +### What tsEVA Is NOT + +- **Not simple detrending**: TS transformation preserves full distribution structure +- **Not stationary EVA on residuals**: Accounts for time-varying parameters properly +- **Not limited to linear trends**: Handles arbitrary time-varying statistics + +### Advantages Over Alternatives + +- **No parametric assumptions** about trend form (uses data-driven moving statistics) +- **Maintains physical interpretability** (parameters remain in original units) +- **Computationally efficient** (no complex likelihood optimization for trend parameters) +- **Robust to model misspecification** (transformation based on empirical statistics) + +### When to Use Other Methods + +- **Bayesian Hierarchical Models**: When prior information is available and important +- **GAMLSS (Generalized Additive Models)**: When covariates (not just time) drive non-stationarity +- **Point Process Models**: When exact timing of events matters + +## Applications + +tsEVA has been used in numerous peer-reviewed climate studies: + +- **Coastal flooding**: Global projections of extreme sea levels (Vousdoukas et al., Nature Communications) +- **Heat waves**: European extremes under warming scenarios (Dosio et al., Environmental Research Letters) +- **Droughts**: Economic impacts of drought intensification (Naumann et al., Nature Climate Change) +- **Wave extremes**: Changes in coastal wave energy (Mentaschi et al., Geophysical Research Letters) +- **River flooding**: Cost-effective adaptation strategies (Dottori et al., Nature Climate Change) + +These applications demonstrate tsEVA's versatility and reliability for high-stakes climate risk assessment. + +--- + +**Key Takeaway**: tsEVA enables rigorous, physically interpretable extreme value analysis under non-stationary conditions by transforming to stationarity, applying established EVA theory, and reverse-transforming results. This approach balances statistical rigor with computational efficiency and physical interpretability. diff --git a/docs/3_Workflow_Patterns.md b/docs/3_Workflow_Patterns.md new file mode 100644 index 0000000..481ef65 --- /dev/null +++ b/docs/3_Workflow_Patterns.md @@ -0,0 +1,444 @@ +# tsEVA 2.0 Workflow Patterns + +This document describes standard analysis workflows and common patterns for using tsEVA 2.0. All code examples use only documented functions. + +## Standard Workflows + +### Workflow 1: Basic Non-Stationary GEV Analysis + +**Use Case**: Analyzing annual maxima with time-varying statistics (trend only, no seasonality) + +```matlab +% 1. Load data +load('timeAndSeriesData.mat'); % Assumes variable timeAndSeries [time, values] + +% 2. Define parameters +timeWindow = 365.25 * 6; % 6 years in days +minPeakDistanceInDays = 3; % Minimum separation between peaks + +% 3. Perform non-stationary EVA +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trend', ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +% 4. Compute return levels +returnPeriods = [10, 50, 100, 500]; % Years +timeInstants = linspace(min(timeAndSeries(:,1)), max(timeAndSeries(:,1)), 100); +[rl] = tsEvaComputeReturnLevelsGEVFromAnalysisObj(returnPeriods, ... + timeInstants, nonStatEvaParams, statTransfData); + +% 5. Visualize results +tsEvaPlotGEVImageScFromAnalysisObj(wr, nonStatEvaParams, statTransfData); +tsEvaPlotReturnLevelsGEVFromAnalysisObj(returnPeriods, nonStatEvaParams, statTransfData); +``` + +**Key Pattern**: `NonStationary → ComputeReturnLevels → Plot` + +--- + +### Workflow 2: Seasonal Non-Stationary Analysis + +**Use Case**: Data with clear seasonal cycle (e.g., sea levels, temperature) + +```matlab +% 1. Load data +load('timeAndSeriesData.mat'); + +% 2. Define parameters +timeWindow = 365.25 * 6; % 6 years +minPeakDistanceInDays = 3; + +% 3. Non-stationary EVA with seasonality +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'seasonal', ... % KEY: Include seasonality + 'minPeakDistanceInDays', minPeakDistanceInDays); + +% 4. Plot seasonal component +tsEvaPlotSeasonalityGevFromAnalysisObj(nonStatEvaParams, statTransfData); + +% 5. Compute and plot return levels +returnPeriods = [10, 50, 100]; +tsEvaPlotReturnLevelsGEVFromAnalysisObj(returnPeriods, nonStatEvaParams, statTransfData); +``` + +**Key Pattern**: Use `'transfType', 'seasonal'` for seasonal data + +--- + +### Workflow 3: GPD (Peaks Over Threshold) + +**Use Case**: When events are not limited to annual maxima (e.g., all storms above threshold) + +```matlab +% 1. Load data +load('timeAndSeriesData.mat'); + +% 2. Define parameters +timeWindow = 365.25 * 50; % 50 years for long-term trends +minPeakDistanceInDays = 60; % Longer separation for independence +thresholdQuantile = 0.95; % Use 95th percentile as threshold + +% 3. Non-stationary EVA (GPD approach) +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trend', ... + 'minPeakDistanceInDays', minPeakDistanceInDays, ... + 'thresholdQuantile', thresholdQuantile); + +% 4. Compute GPD return levels +returnPeriods = [10, 50, 100, 500]; +timeInstants = linspace(min(timeAndSeries(:,1)), max(timeAndSeries(:,1)), 50); +[rl] = tsEvaComputeReturnLevelsGPDFromAnalysisObj(returnPeriods, ... + timeInstants, nonStatEvaParams, statTransfData); + +% 5. Visualize GPD results +tsEvaPlotGPDImageScFromAnalysisObj(wr, nonStatEvaParams, statTransfData); +tsEvaPlotReturnLevelsGPDFromAnalysisObj(returnPeriods, nonStatEvaParams, statTransfData); +``` + +**Key Pattern**: Set `thresholdQuantile` for POT approach, use GPD-specific plotting functions + +--- + +### Workflow 4: Percentile-Based Confidence Interval + +**Use Case**: When extremes change more than the overall distribution (more sensitive to tail behavior) + +```matlab +% 1. Load data +load('timeAndSeriesData.mat'); + +% 2. Define parameters +timeWindow = 365.25 * 6; +minPeakDistanceInDays = 3; +ciPercentile = 98.5; % Use 98.5th percentile for CI estimation + +% 3. Transformation using percentile-based CI +[stSeries] = tsEvaTransformSeriesToStat_ciPercentile(timeAndSeries, ... + timeWindow, ciPercentile); + +% 4. Perform stationary EVA on transformed series +[evaParams] = tsEvaStationary(stSeries, 'minPeakDistanceInDays', minPeakDistanceInDays); + +% 5. Back-transform and compute return levels +% (Manual back-transformation required - consult exampleGenerateSeriesEVAGraphs_ciPercentile.m) +``` + +**Key Pattern**: Use `tsEvaTransformSeriesToStat_ciPercentile` when extremes are changing faster than mean + +--- + +### Workflow 5: Stationary EVA (Baseline) + +**Use Case**: Data without significant trends (or for comparison with non-stationary approach) + +```matlab +% 1. Load data +load('timeAndSeriesData.mat'); + +% 2. Stationary EVA +[evaParams] = tsEvaStationary(timeAndSeries, ... + 'minPeakDistanceInDays', 3); + +% 3. Compute return levels (stationary) +returnPeriods = [10, 50, 100, 500]; +[rl] = tsEvaComputeReturnLevelsGEV(returnPeriods, evaParams.gevParams); + +% 4. Plot (stationary version) +tsEvaPlotReturnLevelsGEVStationary(returnPeriods, evaParams); +``` + +**Key Pattern**: Use `tsEvaStationary` instead of `tsEvaNonStationary` for baseline analysis + +--- + +### Workflow 6: Bivariate Copula Analysis (Non-Stationary) + +**Use Case**: Joint analysis of two related variables (e.g., wave height + storm surge) + +```matlab +% 1. Load data +load('bivariateData.mat'); % timeAndSeries1, timeAndSeries2 + +% 2. Define parameters +timeWindow = 365.25 * 40; % 40 years +minPeakDistanceInDays = 30; +thresholdQuantile1 = 0.95; % 95th percentile for variable 1 +thresholdQuantile2 = 0.99; % 99th percentile for variable 2 +maxTimeLagInDays = 45; % Maximum lag for joint occurrence + +% 3. Sample joint peaks +[peaks] = tsCopulaSampleJointPeaksMultiVariatePruning(... + {timeAndSeries1, timeAndSeries2}, ... + minPeakDistanceInDays, ... + [thresholdQuantile1, thresholdQuantile2], ... + maxTimeLagInDays); + +% 4. Fit copula extremes +copulaFamily = 'Gumbel'; % Or 'Gaussian' or 'Frank' (bivariate only) +[copulaParams, marginalParams] = tsCopulaExtremes(peaks, timeWindow, ... + 'copulaFamily', copulaFamily, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +% 5. Compute bivariate return periods +returnPeriods = [10, 50, 100]; +[rpMatrix] = tsCopulaComputeBivarRP(returnPeriods, copulaParams, marginalParams); + +% 6. Visualize +tsCopulaPlotBivariate(copulaParams, marginalParams); +tsPlotBivarReturnPeriod(rpMatrix, marginalParams); +``` + +**Key Pattern**: +1. Sample joint peaks with `tsCopulaSampleJointPeaksMultiVariatePruning` +2. Fit with `tsCopulaExtremes` +3. Compute return periods with `tsCopulaComputeBivarRP` +4. Visualize with copula-specific plotting functions + +--- + +### Workflow 7: Trivariate Copula Analysis + +**Use Case**: Joint analysis of three variables (e.g., significant wave height, peak period, direction) + +```matlab +% 1. Load data +load('trivariateData.mat'); % timeAndSeries1, timeAndSeries2, timeAndSeries3 + +% 2. Define parameters +timeWindow = 365.25 * 40; +minPeakDistanceInDays = 30; +thresholds = [0.95, 0.95, 0.90]; % Quantile thresholds for each variable +maxTimeLagInDays = 45; + +% 3. Sample joint peaks (trivariate) +[peaks] = tsCopulaSampleJointPeaksMultiVariatePruning(... + {timeAndSeries1, timeAndSeries2, timeAndSeries3}, ... + minPeakDistanceInDays, ... + thresholds, ... + maxTimeLagInDays); + +% 4. Fit copula extremes (Gaussian or Gumbel for multivariate) +copulaFamily = 'Gumbel'; % MUST be Gaussian or Gumbel (not Frank) +[copulaParams, marginalParams] = tsCopulaExtremes(peaks, timeWindow, ... + 'copulaFamily', copulaFamily, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +% 5. Visualize trivariate results +tsCopulaPlotTrivariate(copulaParams, marginalParams); +``` + +**Key Pattern**: +- Use Gaussian or Gumbel copula (NOT Frank) for trivariate analysis +- Adjust thresholds to ensure sufficient joint events (50+ recommended) + +--- + +### Workflow 8: Year Extremes Copula (Annual Maxima Joint Analysis) + +**Use Case**: Analyzing joint distribution of annual maxima across multiple variables + +```matlab +% 1. Load data +load('bivariateData.mat'); + +% 2. Compute annual maxima for each variable +annualMaxima1 = tsEvaComputeAnnualMaxima(timeAndSeries1); +annualMaxima2 = tsEvaComputeAnnualMaxima(timeAndSeries2); + +% 3. Fit year extremes copula +copulaFamily = 'Gumbel'; +[copulaFit] = tsCopulaYearExtrFit({annualMaxima1, annualMaxima2}, copulaFamily); + +% 4. Generate random samples from fitted copula +nSamples = 1000; +[samples] = tsCopulaYearExtrRnd(copulaFit, nSamples); + +% 5. Visualize +tsCopulaYearExtrPlotSctrBivar(copulaFit); +``` + +**Key Pattern**: Use `YearExtr` functions for annual maxima joint analysis + +--- + +## Common Patterns + +### Pattern A: Data Preparation + +```matlab +% Load and inspect +load('rawData.mat'); +timeAndSeries = [time(:), values(:)]; % Ensure column vectors + +% Remove NaNs or fill gaps +timeAndSeries = timeAndSeries(~isnan(timeAndSeries(:,2)), :); +% OR +[filledSeries] = tsEvaFillSeries(timeAndSeries); + +% Check time step +dt = tsEvaGetTimeStep(timeAndSeries); +fprintf('Time step: %.2f days\n', dt); + +% Compute annual maxima if needed +annualMaxima = tsEvaComputeAnnualMaxima(timeAndSeries); +``` + +--- + +### Pattern B: Parameter Selection Guidelines + +```matlab +% Time window selection +dataSpan = max(timeAndSeries(:,1)) - min(timeAndSeries(:,1)); +timeWindow = dataSpan / 5; % Start with 1/5 of data span +% Typical: 5-10 years for short series, 30-50 years for climate projections + +% Peak distance selection +dt = tsEvaGetTimeStep(timeAndSeries); +minPeakDistanceInDays = max(3, 3*dt); % At least 3 time steps +% For monthly data: 30-60 days +% For daily data: 2-5 days +% For sub-daily: 1-3 days + +% Threshold quantile +thresholdQuantile = 0.95; % Gives ~18 events/year for daily data +% Lower (0.90): More events, less bias, more variance +% Higher (0.98): Fewer events, more bias, less variance +``` + +--- + +### Pattern C: Ensemble Analysis + +**Use Case**: Quantify uncertainty by varying parameters + +```matlab +% Define parameter ranges +timeWindows = [365.25*5, 365.25*6, 365.25*8]; % Multiple time windows +thresholdQuantiles = [0.93, 0.95, 0.97]; % Multiple thresholds + +% Run ensemble +[ensembleParams] = tsEnsembleEvaParams(timeAndSeries, ... + timeWindows, thresholdQuantiles); + +% Analyze spread +meanReturnLevel = mean(ensembleParams.returnLevels, 3); +stdReturnLevel = std(ensembleParams.returnLevels, 0, 3); +``` + +--- + +### Pattern D: Model Comparison (Stationary vs Non-Stationary) + +```matlab +% Fit both models +[statParams] = tsEvaStationary(timeAndSeries, 'minPeakDistanceInDays', 3); +[nonStatParams, statTransfData] = tsEvaNonStationary(timeAndSeries, timeWindow, ... + 'transfType', 'trend', 'minPeakDistanceInDays', 3); + +% Compare return levels +returnPeriods = [10, 50, 100]; +timeInstant = max(timeAndSeries(:,1)); % Present day + +rlStat = tsEvaComputeReturnLevelsGEV(returnPeriods, statParams.gevParams); +rlNonStat = tsEvaComputeReturnLevelsGEVFromAnalysisObj(returnPeriods, ... + timeInstant, nonStatParams, statTransfData); + +fprintf('Return Level Comparison (Present):\n'); +fprintf('RP\tStationary\tNon-Stationary\tDifference\n'); +for i = 1:length(returnPeriods) + fprintf('%d\t%.3f\t\t%.3f\t\t%.3f\n', ... + returnPeriods(i), rlStat(i), rlNonStat(i), rlNonStat(i) - rlStat(i)); +end +``` + +--- + +### Pattern E: Goodness-of-Fit Testing + +```matlab +% For copulas +[gofResults] = tsCopulaGOFNonStat(copulaParams, marginalParams, peaks); + +% Visual inspection +% - Check transformation to stationarity +tsEvaPlotTransfToStatFromAnalysisObj(nonStatEvaParams, statTransfData); + +% - Check trend and variance +tsEvaPlotSeriesTrendStdDevFromAnalyisObj(nonStatEvaParams, statTransfData); +``` + +--- + +## Function Call Sequences + +### Monovariate EVA +1. `tsEvaNonStationary` OR `tsEvaStationary` +2. `tsEvaComputeReturnLevels[GEV|GPD]FromAnalysisObj` OR direct GEV/GPD functions +3. `tsEvaPlot*FromAnalysisObj` OR direct plotting functions + +### Copula EVA +1. `tsCopulaSampleJointPeaksMultiVariatePruning` +2. `tsCopulaExtremes` OR `tsCopulaYearExtrFit` +3. `tsCopulaComputeBivarRP` OR `tsCopulaYearExtrGetMltvrtRetPeriod` +4. `tsCopulaPlot*` OR `tsPlotBivarReturnPeriod` + +### Data Transformation (Advanced) +1. `tsEvaTransformSeriesToStationaryTrend*` (if manual transformation needed) +2. `tsEvaStationary` on transformed data +3. Manual back-transformation using transformation parameters + +--- + +## Best Practices + +### DO +- Start with visual inspection of time series +- Test for stationarity before assuming non-stationarity +- Use appropriate time windows relative to data span +- Ensure peak independence through adequate separation +- Report confidence intervals and uncertainty +- Validate results against physical understanding + +### DON'T +- Use Frank copula for >2 variables (multivariate) +- Set time window longer than 1/3 of data span +- Ignore seasonal patterns in seasonal data +- Over-interpret small samples (< 30 years for long-term trends) +- Extrapolate beyond reasonable return periods (max ~3x data span) + +--- + +## Troubleshooting + +### Issue: "Not enough peaks" +**Cause**: Threshold too high or peak separation too long +**Solution**: Lower `thresholdQuantile` or reduce `minPeakDistanceInDays` + +### Issue: "Transformation doesn't look stationary" +**Cause**: Time window too short or data has multiple scales of variability +**Solution**: Increase time window or use seasonal transformation + +### Issue: "Copula fit fails" +**Cause**: Insufficient joint events or incompatible copula family +**Solution**: Check joint peak count (need 50+), verify copula family for dimensionality + +### Issue: "Return levels seem unrealistic" +**Cause**: Extrapolation too far, poor fit, or data quality issues +**Solution**: Check data quality, assess goodness-of-fit, limit return period range + +--- + +## Example Script References + +For complete working examples, see: +- **Monovariate**: `examplesMonovariate/exampleGenerateSeriesEVAGraphs.m` +- **Seasonal**: `examplesMonovariate/exampleGenerateSeriesEVAGraphs_ciPercentile.m` +- **GPD**: `examplesMonovariate/exampleSPISeries.m` +- **Bivariate Copula**: `examplesCopula/exampleCopulaBivariateNonStationary.m` +- **Case Studies**: `examplesCopula/caseStudy01.m`, `caseStudy02.m`, `caseStudy03.m` + +Always adapt these examples to your specific data and research questions. From e362ebf3aed555a3101fbaa41a28cc26fafd67ac Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 30 Jan 2026 17:07:18 +0000 Subject: [PATCH 4/4] Add comprehensive monovariate examples documentation Co-authored-by: menta78 <19606593+menta78@users.noreply.github.com> --- docs/4_Monovariate_Examples.md | 1179 ++++++++++++++++++++++++++++++++ 1 file changed, 1179 insertions(+) create mode 100644 docs/4_Monovariate_Examples.md diff --git a/docs/4_Monovariate_Examples.md b/docs/4_Monovariate_Examples.md new file mode 100644 index 0000000..8e15d35 --- /dev/null +++ b/docs/4_Monovariate_Examples.md @@ -0,0 +1,1179 @@ +# tsEVA 2.0 Monovariate Examples + +This document provides comprehensive walkthroughs of monovariate Extreme Value Analysis (EVA) examples in tsEVA. Each example demonstrates practical applications with real-world data, explaining both the implementation steps and the scientific reasoning behind parameter choices. + +## Overview of Monovariate EVA in tsEVA + +Monovariate Extreme Value Analysis focuses on understanding the statistical behavior of extreme values in a single variable time series. The tsEVA toolbox provides two main approaches: + +1. **GEV (Generalized Extreme Value)**: For block maxima (typically annual maxima) +2. **GPD (Generalized Pareto Distribution)**: For Peaks Over Threshold (POT) analysis + +Both approaches can handle: +- **Stationary data**: Statistics remain constant over time +- **Non-stationary data**: Statistics vary with time (trend and/or seasonality) + +### Key Concepts + +- **Time Window**: The period over which moving statistics are computed (typically 5-30 years) +- **Transformation Types**: + - `'trend'`: Accounts for long-term trends using moving standard deviation + - `'trendCiPercentile'`: Uses moving percentiles (better for extremes, more uncertain) + - `'trendLinear'`: Assumes linear trend in the data + - `'seasonal'`: Includes seasonal cycles in addition to trends + - `'seasonalCiPercentile'`: Seasonal analysis with percentile-based transformation +- **Return Levels**: Values expected to be exceeded with a given probability +- **Return Periods**: Average time between exceedances of a given level + +--- + +## Example 1: Basic Non-Stationary GEV Analysis +### `exampleGenerateSeriesEVAGraphs.m` + +**Purpose**: Comprehensive analysis of residual water levels with time-varying statistics and seasonality. + +**Use Case**: Ocean water levels, wave heights, or any environmental variable with both long-term trends and seasonal patterns. + +### Data Description + +- **Dataset**: `timeAndSeriesHebrides.mat` - 30 years of residual water levels at the Hebrides islands +- **Format**: `timeAndSeries` = [timestamps, water levels (m)] +- **Characteristics**: Non-stationary with trends and seasonal variations + +### Key Parameters + +```matlab +timeWindow = 365.25*6; % 6 years: minimum window for non-stationarity detection +minPeakDistanceInDays = 3; % Ensure peak independence (3 days for water levels) +extremesRange = [.2 1.2]; % Range for plotting extremes (meters) +seasonalExtrRange = [.1 1.1]; % Range for seasonal extremes +``` + +**Why these values?** +- **6-year window**: Below this, series is considered quasi-stationary; above this, detects long-term changes +- **3-day separation**: Water level storms typically last 1-2 days; 3 days ensures independence +- **Ranges**: Set based on physical bounds of the phenomenon + +### Step-by-Step Walkthrough + +#### Part 1: Trend-Only Analysis + +```matlab +% Load data +addpath('../'); +load('timeAndSeriesHebrides.mat'); +timeAndSeries = timeAndSeriesHebrides; + +% Define parameters +timeWindow = 365.25*6; +minPeakDistanceInDays = 3; + +% Perform non-stationary EVA with trend transformation +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trend', ... + 'minPeakDistanceInDays', minPeakDistanceInDays); +``` + +**What happens here?** +1. Time series is divided into overlapping windows +2. For each window: compute mean and standard deviation +3. Transform to stationarity by removing trend +4. Fit GEV and GPD to transformed extremes +5. Back-transform to original scale + +#### Visualization: Series with Trend + +```matlab +% Plot original series with detected trend and variance +hndl = tsEvaPlotSeriesTrendStdDevFromAnalyisObj(nonStatEvaParams, ... + statTransfData, ... + 'ylabel', 'Lvl (m)', ... + 'title', 'Hebrides', ... + 'dateformat', 'yy', ... + 'xtick', tickTmStmp); +saveas(hndl{1}, 'seriesTrendOnly.png'); +``` + +**Key Insight**: This plot shows the original data overlaid with: +- Moving mean (captures long-term trends) +- Moving standard deviation bands (captures changing variability) + +#### Visualization: Time-Varying GEV and GPD + +```matlab +% Define range for extreme values +wr = linspace(min(extremesRange), max(extremesRange), 1501); + +% Plot 2D GEV (time on x-axis, level on y-axis, probability as color) +hndl = tsEvaPlotGEVImageScFromAnalysisObj(wr, nonStatEvaParams, ... + statTransfData, ... + 'ylabel', 'Lvl (m)', ... + 'dateformat', 'yy', ... + 'xtick', tickTmStmp); +title('GEV'); +saveas(hndl{1}, 'GEV2DTrendOnly.png'); + +% Plot 2D GPD +hndl = tsEvaPlotGPDImageScFromAnalysisObj(wr, nonStatEvaParams, ... + statTransfData, ... + 'ylabel', 'Lvl (m)', ... + 'dateformat', 'yy', ... + 'xtick', tickTmStmp); +title('GPD'); +saveas(hndl{1}, 'GPD2DTrendOnly.png'); +``` + +**Interpretation**: +- **GEV**: Shows probability distribution of annual maxima changing over time +- **GPD**: Shows probability distribution of threshold exceedances changing over time +- Color indicates probability density +- Horizontal slices = distribution at specific time +- Vertical slices = how a specific level's probability changes over time + +#### Computing Return Levels + +```matlab +% Choose a time instant (index 1000 in the time series) +timeIndex = 1000; +timeStamps = statTransfData.timeStamps; + +% Plot GEV return levels at this time instant +hndl = tsEvaPlotReturnLevelsGEVFromAnalysisObj(nonStatEvaParams, ... + timeIndex, ... + 'ylim', [.5 1.5]); +saveas(hndl{1}, 'GEV_ReturnLevels.png'); + +% Plot GPD return levels +hndl = tsEvaPlotReturnLevelsGPDFromAnalysisObj(nonStatEvaParams, ... + timeIndex, ... + 'ylim', [.5 1.5]); +saveas(hndl{1}, 'GPD_ReturnLevels.png'); +``` + +**Understanding Return Levels**: +- X-axis: Return period in years (e.g., 10, 50, 100) +- Y-axis: Water level (m) +- Shaded area: Confidence interval +- A 100-year return level of 1.2m means: this level is exceeded on average once every 100 years + +**Note**: "For GEV the sample is small and the confidence interval is broad" - annual maxima give fewer data points than POT approach, leading to wider uncertainty. + +#### Transformation Diagnostic + +```matlab +% Plot the transformed stationary series +hndl = tsEvaPlotTransfToStatFromAnalysisObj(nonStatEvaParams, ... + statTransfData, ... + 'dateformat', 'yy', ... + 'xtick', tickTmStmp); +saveas(hndl{1}, 'statSeriesTrendOnly.png'); +``` + +**Purpose**: Verify transformation quality - transformed series should appear stationary (constant mean and variance). + +#### Part 2: Seasonal Analysis + +```matlab +% Perform non-stationary EVA with seasonal transformation +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'seasonal', ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +% Define new range for seasonal extremes +wr = linspace(min(seasonalExtrRange), max(seasonalExtrRange), 1501); +``` + +**Why seasonal analysis?** +- Water levels often show seasonal patterns (winter storms vs. summer calm) +- Seasonal transformation removes both trend AND annual cycle +- Provides more accurate extreme value estimates + +```matlab +% Plot a time slice (1988-1993) +slice = {1988, 1993}; + +% Plot series for this period +hndl = tsEvaPlotSeriesTrendStdDevFromAnalyisObj(nonStatEvaParams, ... + statTransfData, ... + 'ylabel', 'Lvl (m)', ... + 'title', '1988-1993', ... + 'minyear', slice{1}, ... + 'maxyear', slice{2}); +saveas(hndl{1}, 'seriesSeasonal.png'); + +% Plot 3D GEV for this period +hndl = tsEvaPlotGEV3DFromAnalysisObj(wr, nonStatEvaParams, ... + statTransfData, ... + 'xlabel', 'Lvl (m)', ... + 'minyear', slice{1}, ... + 'maxyear', slice{2}); +title('GEV 3D, 1988-1993'); +saveas(hndl{1}, 'GEV3DSeasonal.png'); +``` + +**3D Plot Interpretation**: +- X-axis: Water level +- Y-axis: Time +- Z-axis: Probability density +- Shows how distribution evolves through time and seasons + +### Expected Outputs + +1. **seriesTrendOnly.png**: Original data with trend line and variability bands +2. **GEV2DTrendOnly.png**: Time-varying GEV distribution +3. **GPD2DTrendOnly.png**: Time-varying GPD distribution +4. **GEV_ReturnLevels.png**: Return level plot for GEV at specific time +5. **GPD_ReturnLevels.png**: Return level plot for GPD at specific time +6. **statSeriesTrendOnly.png**: Transformed stationary series +7. **seriesSeasonal.png**: Series with seasonal component for 1988-1993 +8. **GEV3DSeasonal.png**: 3D visualization of seasonal GEV +9. **GEV2DSeasonal.png**: 2D GEV with seasonality +10. **GPD2DSeasonal.png**: 2D GPD with seasonality + +### Key Takeaways + +1. **Two-stage approach**: First analyze trend-only, then add seasonality +2. **Always plot diagnostics**: Check that transformation produces stationary series +3. **GEV vs GPD tradeoff**: GEV has wider confidence intervals but is theoretically rigorous for maxima; GPD uses more data but requires threshold selection +4. **Time windows matter**: Too small = noisy estimates; too large = miss real changes +5. **Return levels are time-dependent**: In non-stationary analysis, a 100-year event today may be different from a 100-year event 20 years from now + +--- + +## Example 2: Percentile-Based Confidence Intervals +### `exampleGenerateSeriesEVAGraphs_ciPercentile.m` + +**Purpose**: Use moving percentiles instead of moving standard deviation to capture extreme value variability. + +**Use Case**: When extremes show different variability patterns than the bulk of the data, or when data has heavy tails. + +### Why Percentiles? + +Standard deviation captures variability of the entire distribution, but extremes may behave differently. Using a high percentile (e.g., 98th) directly captures extreme value behavior. + +**Trade-off**: +- ✅ Better models extreme variability +- ⚠️ Higher uncertainty (fewer data points at high percentiles) + +### Key Parameters + +```matlab +ciPercentile = 98; % Use 98th percentile to capture extreme variability +``` + +**Choice of percentile**: +- Too low (e.g., 90%): Not capturing true extremes +- Too high (e.g., 99.5%): Too few points, high uncertainty +- Typical range: 95-99% + +### Code Walkthrough + +```matlab +% Trend-only analysis with percentile-based CI +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trendCiPercentile', ... + 'ciPercentile', ciPercentile, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); +``` + +**What's different**: +- `'transfType', 'trendCiPercentile'` instead of `'trend'` +- Must specify `'ciPercentile'` parameter +- Transformation uses moving 98th percentile instead of moving std dev + +### Computing and Displaying Return Levels + +```matlab +timeIndex = 1000; +timeStamps = statTransfData.timeStamps; + +% Compute return levels explicitly +[rlevGEV, rlevGEVErr] = tsEvaComputeReturnLevelsGEVFromAnalysisObj(... + nonStatEvaParams, [10, 20, 50, 100], ... + 'timeindex', timeIndex); +rlevGEV + +% Plot with custom range +hndl = tsEvaPlotReturnLevelsGEVFromAnalysisObj(nonStatEvaParams, ... + timeIndex, ... + 'ylim', rlRange); +title(['GEV return levels for ' datestr(tmstmpref)]); +saveas(hndl{1}, 'GEV_ReturnLevels_ciPercentile.png'); +``` + +**Note**: This example also computes return levels explicitly before plotting, allowing you to inspect the numerical values. + +### Seasonal Analysis with Percentiles + +```matlab +% Seasonal analysis with percentile-based transformation +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'seasonalCiPercentile', ... + 'ciPercentile', ciPercentile, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); +``` + +Combines both seasonality AND percentile-based variability estimation. + +### Expected Outputs + +1. **seriesTrendOnly_ciPercentile.png**: Series with 98th percentile bands (instead of std dev) +2. **GEV2DTrendOnly_ciPercentile.png**: GEV with percentile transformation +3. **GEV_ReturnLevels_ciPercentile.png**: Return levels with percentile-based uncertainty +4. Similar plots for GPD and seasonal analysis + +### Key Takeaways + +1. **Choose transformation based on data**: If extremes show different variability patterns, use percentile approach +2. **Expect wider uncertainty**: Percentile-based methods have larger confidence intervals +3. **Compare approaches**: Run both standard deviation and percentile methods to understand sensitivity +4. **Percentile selection**: Test different percentiles (95, 98, 99) to assess robustness + +--- + +## Example 3: Comparing Different CI Approaches +### `exampleCompareDifferentCI.m` + +**Purpose**: Systematically compare moving standard deviation vs. multiple percentile-based transformations. + +**Use Case**: Determining which transformation approach is most appropriate for your data. + +### Data Description + +- **Dataset**: `timeAndSeries_waves_015_220E_055_509N.mat` - Wave height time series +- Can also use alternative wave datasets (commented in code) + +### Key Parameters + +```matlab +timeWindow = 30*365.25; % 30 years (longer for wave data) +minPeakDistanceInDays = 3; +``` + +**Why 30 years?** Wave climate changes occur over longer timescales than water level variations. + +### Step-by-Step Comparison + +#### Method 1: Standard Deviation + +```matlab +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +% Define plotting range based on data +minext = (max(statTransfData.nonStatSeries) + 3*min(statTransfData.nonStatSeries))/4; +maxext = max(statTransfData.nonStatSeries)*1.2; +xext = minext:.01:maxext; + +% Plot GEV +tsEvaPlotGEVImageScFromAnalysisObj(xext, nonStatEvaParams, ... + statTransfData, 'ylabel', 'Hs (m)'); +title('standard deviation'); + +% Transformation diagnostic +tsEvaPlotTransfToStatFromAnalysisObj(nonStatEvaParams, statTransfData); +title('standard deviation'); +``` + +#### Method 2: 98th Percentile + +```matlab +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trendCiPercentile', ... + 'ciPercentile', 98, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +tsEvaPlotGEVImageScFromAnalysisObj(xext, nonStatEvaParams, ... + statTransfData, 'ylabel', 'Hs (m)'); +title('98th percentile'); + +tsEvaPlotTransfToStatFromAnalysisObj(nonStatEvaParams, statTransfData); +title('98th percentile'); +``` + +#### Method 3: 98.5th Percentile + +```matlab +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trendCiPercentile', ... + 'ciPercentile', 98.5, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +tsEvaPlotGEVImageScFromAnalysisObj(xext, nonStatEvaParams, ... + statTransfData, 'ylabel', 'Hs (m)'); +title('98.5th percentile'); + +tsEvaPlotTransfToStatFromAnalysisObj(nonStatEvaParams, statTransfData); +title('98.5th percentile'); +``` + +#### Method 4: 99th Percentile + +```matlab +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trendCiPercentile', ... + 'ciPercentile', 99, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +tsEvaPlotGEVImageScFromAnalysisObj(xext, nonStatEvaParams, ... + statTransfData, 'ylabel', 'Hs (m)'); +title('99th percentile'); + +tsEvaPlotTransfToStatFromAnalysisObj(nonStatEvaParams, statTransfData); +title('99th percentile'); +``` + +### How to Compare Results + +For each method, examine: + +1. **GEV plots**: Are the patterns similar or different? + - Similar patterns → results are robust + - Very different patterns → data sensitive to transformation choice + +2. **Transformation diagnostic plots**: Is the transformed series truly stationary? + - Should show constant mean and variance + - If not, transformation is inadequate + +3. **Width of extreme value bands**: How do uncertainty estimates compare? + - Percentile methods typically show wider bands + - Very narrow bands may underestimate uncertainty + +### Expected Outputs + +- 4 pairs of plots (GEV + transformation diagnostic) for each method +- Visual comparison of how transformation choice affects results + +### Key Takeaways + +1. **No universal best method**: Choice depends on data characteristics +2. **Robustness check**: If results dramatically change with percentile choice, be cautious +3. **Document your choice**: Explain why you selected a particular transformation +4. **Consider physics**: Choose based on understanding of the phenomenon, not just statistics +5. **Standard deviation is safer**: If unsure, start with standard deviation (more stable) +6. **Percentiles for heavy tails**: Use percentile approach when data has very heavy tails or extreme outliers + +--- + +## Example 4: GPD/POT for Drought Analysis +### `exampleSPISeries.m` + +**Purpose**: Analyze Standardized Precipitation Index (SPI) using GPD only (no GEV), with long peak separation. + +**Use Case**: Phenomena where: +- Events are rare and widely separated +- Annual maxima approach is inappropriate +- GPD/POT is more suitable than GEV + +### Data Description + +- **Dataset**: `timeAndSeries_SPI_179_750E_-16.750N.mat` - SPI time series +- **Variable**: SPI (Standardized Precipitation Index, inverted to analyze droughts as peaks) +- **Special characteristic**: Drought events are separated by months, not days + +### Why SPI is Different + +```matlab +% Invert series: analyze droughts (negative SPI) as positive peaks +timeAndSeries(:,2) = -timeAndSeries(:,2); +``` + +SPI is already normalized and standardized. Negative values indicate drought. By inverting, we analyze drought severity as if they were peaks. + +### Key Parameters + +```matlab +timeWindow = 50*315.25; % 50 years: very long for climate analysis +minPeakDistanceInDays = 5*30.2; % ~5 months: droughts are long-lasting events +returnPeriodsInYears = [10 20 50 100]; +``` + +**Critical choices**: +- **5-month separation**: Drought events typically last months; ensures independence +- **50-year window**: Climate trends emerge over multi-decadal timescales +- **GPD only**: With such long peak separation, "peaks per year" is very small; POT/GPD approach is more appropriate than annual maxima + +### Code Walkthrough + +```matlab +[nonStationaryEvaParams, stationaryTransformData, isValid] = tsEvaNonStationary(... + timeAndSeries, ... + timeWindow, ... + 'minPeakDistanceInDays', minPeakDistanceInDays, ... + 'transfType', 'trendCIPercentile', ... + 'cipercentile', 80, ... + 'potPercentiles', 80, ... % Use 80th percentile as threshold + 'evdType', 'GPD'); % ONLY fit GPD, not GEV +``` + +**Key parameters explained**: +- `'evdType', 'GPD'`: Skip GEV fitting entirely +- `'potPercentiles', 80`: Use only one threshold (80th percentile) for POT +- `'cipercentile', 80`: Match transformation percentile to POT threshold + +**Why 80th percentile?** +- Lower than typical (95-98%) because SPI is already normalized +- Captures sufficient drought events for GPD fitting +- Higher values would give too few events + +### Visualization + +```matlab +tsEvaPlotSeriesTrendStdDevFromAnalyisObj(nonStationaryEvaParams, ... + stationaryTransformData, ... + 'plotpercentile', 95., ... + 'ylabel', '-SPI', ... + 'legendLocation', 'southwest'); +``` + +**Note**: `'plotpercentile', 95.` adds 95th percentile line to plot for reference, independent of the 80th percentile used in analysis. + +### Computing Return Levels + +```matlab +[returnLevels, returnLevelsErr] = tsEvaComputeReturnLevelsGPDFromAnalysisObj(... + nonStationaryEvaParams, returnPeriodsInYears); + +% Invert back to original SPI scale +returnLevels = returnLevels * -1 +``` + +**Interpretation**: A 100-year return level of -2.5 means a drought of SPI = -2.5 is expected once per century. + +### Expected Outputs + +- Series plot showing inverted SPI with trend and percentile bands +- Return level values for 10, 20, 50, 100-year drought events + +### Key Takeaways + +1. **Adapt method to data**: Long-separated events need different parameters +2. **GPD without GEV**: Valid when annual maxima concept doesn't apply +3. **Single threshold**: Can use just one percentile for POT (not the typical multi-percentile exploration) +4. **Transform for interpretation**: Remember to invert results back to original scale +5. **Climate timescales**: Use very long time windows (50+ years) for climate data + +--- + +## Example 5: GEV-Only Analysis for Temperature Extremes +### `exampleTASMaxSeries.m` + +**Purpose**: Analyze yearly maximum temperatures (heat waves) using GEV exclusively. + +**Use Case**: When data consists of one value per year (annual maxima) - GPD/POT is meaningless. + +### Data Description + +- **Dataset**: `timeAndSeriesTASMax.mat` - Annual maximum air surface temperature +- **Format**: One value per year (1850-2100 or similar) +- **Application**: Climate change impact on heat waves + +### Why GEV Only? + +```matlab +'evdType', 'GEV' % Only fit GEV, skip GPD +``` + +When you have exactly one value per year (annual maxima): +- GEV is theoretically appropriate (models maxima directly) +- GPD requires threshold exceedances within years → not applicable here +- No need to select peaks or thresholds + +### Key Parameters + +```matlab +timeWindow = 50*315.25; % 50 years +minPeakDistanceInDays = 5*30.2; % Not really used (already annual maxima) +returnPeriodsInYears = [20 50 100 300]; +``` + +### Code Walkthrough + +```matlab +[nonStationaryEvaParams, stationaryTransformData, isValid] = tsEvaNonStationary(... + timeAndSeries, ... + timeWindow, ... + 'minPeakDistanceInDays', minPeakDistanceInDays, ... + 'extremeLowThreshold', .1, ... % Exclude very low values (data quality) + 'evdType', 'GEV'); % GEV only +``` + +**New parameter**: `'extremeLowThreshold', .1` - excludes the lowest 10% of values, useful for filtering erroneous or missing data. + +### Visualization: Time Series + +```matlab +tsEvaPlotSeriesTrendStdDevFromAnalyisObj(nonStationaryEvaParams, ... + stationaryTransformData, ... + 'ylabel', 'TAS', ... + 'legendLocation', 'northwest'); +ylim([0 40]); +text(datenum(2060, 1, 1), 35, 'Series and trends', 'fontsize', 25); +``` + +Shows temperature data with trends - clear for climate change visualization. + +### Visualization: Time-Varying GEV + +```matlab +tsEvaPlotGEVImageScFromAnalysisObj((0:.001:40)', ... + nonStationaryEvaParams, ... + stationaryTransformData, ... + 'ylabel', 'TAS'); +text(datenum(1980, 1, 1), 35, 'Time varying GEV', 'fontsize', 30); +``` + +Shows how the probability of extreme heat changes over time (climate change impact). + +### Return Levels at Different Time Points + +```matlab +% Early period (1995) +timeIndex = 26; +rlRange = [0 14]; +hndl = tsEvaPlotReturnLevelsGEVFromAnalysisObj(nonStationaryEvaParams, ... + timeIndex, ... + 'ylim', rlRange, ... + 'ylabel', 'return levels (TAS)'); +ax = gca; +ax.YTick = 0:2:16; +text(7, 12.5, 'Return level 1995', 'fontsize', 30); + +% Late period (2095) +timeIndex = size(timeAndSeries, 1) - 4; +rlRange = [0 70]; +hndl = tsEvaPlotReturnLevelsGEVFromAnalysisObj(nonStationaryEvaParams, ... + timeIndex, ... + 'ylim', rlRange, ... + 'ylabel', 'return levels (TAS)'); +text(7, 65, 'Return level 2095', 'fontsize', 30); +``` + +**Key insight**: Note the dramatic difference in y-axis ranges: +- 1995: 0-14°C +- 2095: 0-70°C + +This illustrates the increasing severity of heat wave extremes under climate change. + +### Computing Return Levels + +```matlab +[returnLevels, returnLevelsErr] = tsEvaComputeReturnLevelsGEVFromAnalysisObj(... + nonStationaryEvaParams, returnPeriodsInYears); +``` + +Provides numerical values for 20, 50, 100, 300-year heat wave temperatures. + +### Expected Outputs + +1. Series plot showing temperature trends +2. Time-varying GEV distribution image +3. Return level plots for 1995 and 2095 showing dramatic increase +4. Numerical return level values + +### Key Takeaways + +1. **Annual maxima = GEV**: When data is already annual maxima, use GEV exclusively +2. **Skip POT**: No need for threshold selection or peak finding +3. **Climate change visualization**: Time-varying return levels clearly show changing risk +4. **Long time windows**: Essential for capturing climate trends (50+ years) +5. **Projections**: Can analyze future projections if data includes model outputs + +--- + +## Example 6: Stationary Analysis +### `exampleEVAStationary.m` + +**Purpose**: Apply traditional stationary EVA when time-varying statistics are not needed. + +**Use Case**: +- Short time series where non-stationarity can't be reliably estimated +- Exploratory analysis before attempting non-stationary approach +- Data with no apparent trends + +### Key Difference + +```matlab +statEvaParams = tsEvaStationary(timeAndSeries, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); +``` + +Uses `tsEvaStationary` instead of `tsEvaNonStationary` - fits a single GEV and GPD to entire series. + +### Code Walkthrough + +#### Basic Stationary Analysis + +```matlab +% Fit stationary GEV and GPD +statEvaParams = tsEvaStationary(timeAndSeries, ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +% Compute return levels (single set of values, not time-varying) +[rlevGEV, rlevGEVErr] = tsEvaComputeReturnLevelsGEVFromAnalysisObj(... + statEvaParams, [10, 20, 50, 100]); +rlevGEV + +% Plot return levels (timeIndex=1 because there's only one time point) +hndl = tsEvaPlotReturnLevelsGEVFromAnalysisObj(statEvaParams, 1, ... + 'ylim', [.5 1.5]); +title('GEV'); +saveas(hndl{1}, 'GEV_ReturnLevels_STATIONARY.png'); +``` + +#### Fixed Threshold POT + +```matlab +% Manually specify POT threshold +potThreshold = prctile(timeAndSeries(:,2), 98); + +statEvaParams = tsEvaStationary(timeAndSeries, ... + 'minPeakDistanceInDays', minPeakDistanceInDays, ... + 'doSampleData', false, ... % Skip automatic threshold selection + 'potThreshold', potThreshold); % Use fixed threshold + +[rlevGPD, rlevGPDErr] = tsEvaComputeReturnLevelsGPDFromAnalysisObj(... + statEvaParams, [10, 20, 50, 100]); +rlevGPD +``` + +**Why fixed threshold?** +- Gives full control over threshold selection +- Useful for sensitivity analysis +- Can align with physical thresholds (e.g., flood stage) + +#### Gumbel Distribution + +```matlab +% Fit Gumbel (special case of GEV with shape parameter = 0) +statEvaParams = tsEvaStationary(timeAndSeries, ... + 'minPeakDistanceInDays', minPeakDistanceInDays, ... + 'gevtype', 'gumbel'); + +[rlevGEV, rlevGEVErr] = tsEvaComputeReturnLevelsGEVFromAnalysisObj(... + statEvaParams, [10, 20, 50, 100]); +rlevGEV + +hndl = tsEvaPlotReturnLevelsGEVFromAnalysisObj(statEvaParams, 1, ... + 'ylim', [.5 1.5]); +title('Gumbel'); +``` + +**When to use Gumbel?** +- When shape parameter is close to zero +- For simplicity (fewer parameters to estimate) +- When data suggests exponential tail behavior + +### Expected Outputs + +- Return level plots showing single curve with confidence intervals +- Numerical return level values for GEV, GPD, and Gumbel + +### Key Takeaways + +1. **Simpler but limited**: Stationary analysis is easier but assumes no temporal changes +2. **Useful for validation**: Compare stationary vs non-stationary to assess if complexity is justified +3. **Fixed thresholds**: Give more control in POT analysis +4. **Gumbel option**: Reduces parameters when appropriate +5. **Return levels constant**: Unlike non-stationary case, return levels don't vary with time + +--- + +## Example 7: Linear Trend Analysis +### `exampleGenerateSeriesEVAGraphs_trendLinear.m` + +**Purpose**: Analyze extremes when a linear trend is appropriate (common in climate change studies). + +**Use Case**: Data showing steady increase/decrease over time without complex variations. + +### Key Parameters + +```matlab +timeWindow = 365.25*30; % 30 years +ciPercentile = 99; % 99th percentile for transformation +``` + +### Code Walkthrough + +```matlab +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trendlinear', ... + 'ciPercentile', ciPercentile, ... + 'potPercentiles', [97:0.5:99], ... % Multiple thresholds for robustness + 'minPeakDistanceInDays', minPeakDistanceInDays); +``` + +**What's different**: +- `'transfType', 'trendlinear'`: Assumes linear change in mean +- `'potPercentiles', [97:0.5:99]`: Tests multiple thresholds (97%, 97.5%, 98%, ..., 99%) + +### Extracting Statistics + +```matlab +epsilon = nonStatEvaParams(2).parameters.epsilon; % GPD shape parameter +pvalue = statTransfData.pValueChange; % Significance of trend +``` + +Useful for assessing: +- **epsilon**: Tail behavior (negative = bounded, zero = exponential, positive = heavy) +- **pvalue**: Is the detected trend statistically significant? + +### Multiple Time Point Analysis + +```matlab +% Beginning of series +timeIndex = 1; +[rlevGPD, rlevGPDErr] = tsEvaComputeReturnLevelsGPDFromAnalysisObj(... + nonStatEvaParams, [5, 10, 30, 100], 'timeindex', timeIndex); +tsEvaPlotReturnLevelsGPDFromAnalysisObj(nonStatEvaParams, timeIndex, ... + 'ylim', rlRange, 'maxReturnPeriodYears', 200); +title('GPD return levels-begin'); + +% End of series +timeIndex = size(timeAndSeries, 1); +[rlevGPD, rlevGPDErr] = tsEvaComputeReturnLevelsGPDFromAnalysisObj(... + nonStatEvaParams, [5, 10, 30, 100], 'timeindex', timeIndex); +tsEvaPlotReturnLevelsGPDFromAnalysisObj(nonStatEvaParams, timeIndex, ... + 'ylim', rlRange, 'maxReturnPeriodYears', 200); +title('GPD return levels-end'); +``` + +Directly compares how return levels have changed from beginning to end of time series. + +### Visualization of POT + +```matlab +tsPlotSeriesPotGPDRetLevFromAnalysisObj(nonStatEvaParams, statTransfData); +title('PotAndReturnLevelsLinearTrend'); +``` + +Shows original series with: +- POT threshold(s) +- Selected peaks +- Return level curves + +### Expected Outputs + +1. Series plot with linear trend +2. POT visualization with threshold and peaks +3. 2D GPD distribution over time +4. Return level plots for beginning and end of series +5. Numerical epsilon and p-value for trend assessment + +### Key Takeaways + +1. **Linear trend assumption**: Simpler than general trend but may not capture complex variations +2. **Multiple thresholds**: Using a range of percentiles improves robustness +3. **Shape parameter**: Check epsilon to understand tail behavior +4. **Statistical significance**: Use p-value to validate that trend is real +5. **Temporal comparison**: Plot return levels at multiple times to quantify change + +--- + +## Example 8: GPD with Negative Shape Parameter +### `exampleGenerateSeriesEVAGraphs_gpdNegShapeParam.m` + +**Purpose**: Handle cases where GPD shape parameter is negative (bounded upper tail). + +**Use Case**: Variables with physical upper bounds or when fitting suggests negative shape. + +### Key Parameters + +```matlab +timeWindow = 365.25*15; % 15 years +minPeakDistanceInDays = 14; % 2 weeks (longer than typical) +ciPercentile = 99; +potPercentiles = [97:0.5:99]; % Multiple thresholds +``` + +### Code Walkthrough + +```matlab +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(timeAndSeries, ... + timeWindow, ... + 'transfType', 'trendlinear', ... + 'ciPercentile', ciPercentile, ... + 'potPercentiles', [97:0.5:99], ... + 'minPeakDistanceInDays', minPeakDistanceInDays); + +% Extract shape parameter +epsilon = nonStatEvaParams(2).parameters.epsilon; +``` + +**Negative epsilon interpretation**: +- Distribution has an upper bound +- Upper bound = location + scale/abs(epsilon) +- Return levels eventually plateau +- Extrapolation beyond observed range is unreliable + +### Visualization with Shape Parameter + +```matlab +tsEvaPlotReturnLevelsGPDFromAnalysisObj(nonStatEvaParams, timeIndex, ... + 'ylim', rlRange, 'maxReturnPeriodYears', 200); +hold on; +legend(sprintf('negative shapeparemeter = %.3e', epsilon), ... + 'Location', 'northwest'); +title('GPD return levels-begin'); +``` + +Annotates plot with actual shape parameter value for transparency. + +### Comparing Beginning and End + +```matlab +% Beginning (timeIndex = 1) +epsilon_begin = nonStatEvaParams(2).parameters.epsilon; +[rlevGPD_begin, rlevGPDErr_begin] = tsEvaComputeReturnLevelsGPDFromAnalysisObj(... + nonStatEvaParams, [5, 10, 30, 100], 'timeindex', 1); + +% End (timeIndex = last) +epsilon_end = nonStatEvaParams(2).parameters.epsilon; +[rlevGPD_end, rlevGPDErr_end] = tsEvaComputeReturnLevelsGPDFromAnalysisObj(... + nonStatEvaParams, [5, 10, 30, 100], 'timeindex', timeIndex); +``` + +If negative shape parameter, check if: +- Upper bound is changing over time +- Bound is consistent with physical constraints + +### Expected Outputs + +1. Series plot showing bounded behavior +2. POT and return level visualization +3. 2D GPD showing upper bound constraint +4. Return level plots with shape parameter annotated +5. Comparison of bounds at beginning and end + +### Key Takeaways + +1. **Negative shape is valid**: Not an error, indicates bounded distribution +2. **Physical interpretation**: Makes sense for variables with natural upper limits +3. **Extrapolation caution**: Don't extrapolate far beyond observed data with negative shape +4. **Upper bound formula**: Can compute exact upper bound from parameters +5. **Time-varying bounds**: In non-stationary case, upper bound may change over time + +--- + +## Example 9: Gumbel for SPI +### `exampleSPISeries_Gumbel.m` + +**Purpose**: Use Gumbel distribution (GEV with shape=0) for SPI drought analysis. + +**Use Case**: When shape parameter is close to zero or when simpler model is preferred. + +### Key Difference from Example 4 + +```matlab +[nonStationaryEvaParams, stationaryTransformData, isValid] = tsEvaNonStationary(... + timeAndSeries, ... + timeWindow, ... + 'minPeakDistanceInDays', minPeakDistanceInDays, ... + 'transfType', 'trendCIPercentile', ... + 'cipercentile', 80, ... + 'gevType', 'gumbel', ... % Use Gumbel instead of full GEV + 'evdType', 'GEV'); % Fit GEV (Gumbel variant), not GPD +``` + +**Gumbel vs. full GEV**: +- Gumbel: 2 parameters (location, scale) +- GEV: 3 parameters (location, scale, shape) +- Gumbel: Simpler, more stable with limited data +- GEV: More flexible, can capture heavier or lighter tails + +### When to Choose Gumbel + +1. **Limited data**: Fewer parameters = more stable estimates +2. **Shape ≈ 0**: If full GEV gives shape parameter near zero, Gumbel is appropriate +3. **Theoretical grounds**: Some phenomena have exponential tail behavior +4. **Simplicity**: When interpretability is important + +### Code Walkthrough + +```matlab +% Invert SPI for drought analysis +timeAndSeries(:,2) = -timeAndSeries(:,2); + +% Fit Gumbel with percentile transformation +[nonStationaryEvaParams, stationaryTransformData, isValid] = tsEvaNonStationary(... + timeAndSeries, timeWindow, ... + 'minPeakDistanceInDays', 5*30.2, ... + 'transfType', 'trendCIPercentile', ... + 'cipercentile', 80, ... + 'gevType', 'gumbel', ... + 'evdType', 'GEV'); + +% Plot series +tsEvaPlotSeriesTrendStdDevFromAnalyisObj(nonStationaryEvaParams, ... + stationaryTransformData, ... + 'plotpercentile', 95., ... + 'ylabel', '-SPI', ... + 'legendLocation', 'southwest'); + +% Compute return levels (Gumbel-based) +[returnLevels, returnLevelsErr] = tsEvaComputeReturnLevelsGEVFromAnalysisObj(... + nonStationaryEvaParams, returnPeriodsInYears); + +% Invert back to SPI scale +returnLevels = returnLevels * -1 +``` + +### Expected Outputs + +- Series plot with Gumbel-based transformation +- Return level values for drought events (inverted SPI) + +### Key Takeaways + +1. **Gumbel = simplified GEV**: Appropriate when shape parameter is zero +2. **Stability**: More stable parameter estimation with limited data +3. **Compare with full GEV**: Check if simplification is justified +4. **SPI application**: Particularly relevant for normalized indices like SPI +5. **Use GEV functions**: Even though it's Gumbel, use GEV computation/plotting functions + +--- + +## Common Patterns Across Examples + +### 1. Data Loading and Setup + +```matlab +addpath('../'); % Add tsEVA functions to path +load('datafile.mat'); % Load time series data +timeAndSeries = loadedVariable; % Ensure format: [time, values] +``` + +### 2. Parameter Selection Logic + +| Parameter | Typical Range | Selection Criteria | +|-----------|---------------|-------------------| +| `timeWindow` | 5-50 years | Longer for climate data, shorter for engineering | +| `minPeakDistanceInDays` | 1-150 days | Based on event duration and independence | +| `ciPercentile` | 95-99% | Higher for better extreme capture, but more uncertainty | +| `potPercentiles` | 95-99% or range | Single value or range for robustness testing | +| `evdType` | 'GEV', 'GPD', or both | GEV for annual maxima, GPD for POT, both for comparison | + +### 3. Analysis Flow + +```matlab +% 1. Non-stationary analysis +[nonStatEvaParams, statTransfData] = tsEvaNonStationary(...); + +% 2. Plot series with trends +tsEvaPlotSeriesTrendStdDevFromAnalyisObj(...); + +% 3. Visualize distributions +tsEvaPlotGEVImageScFromAnalysisObj(...); % and/or GPD + +% 4. Compute return levels +[returnLevels, returnLevelsErr] = tsEvaComputeReturnLevelsGEVFromAnalysisObj(...); + +% 5. Plot return levels +tsEvaPlotReturnLevelsGEVFromAnalysisObj(...); + +% 6. Check transformation +tsEvaPlotTransfToStatFromAnalysisObj(...); +``` + +### 4. Scientific Reasoning + +**Time Window Selection**: +- Too small → Noisy estimates, can't detect long-term trends +- Too large → Over-smoothing, miss real variations +- Rule of thumb: 5-10 years for engineering, 30-50 years for climate + +**Peak Distance Selection**: +- Must ensure independence of extreme events +- Typical event duration + safety margin +- Examples: 3 days for water levels, 14 days for large storms, 150 days for droughts + +**Transformation Type**: +- Start with `'trend'` (standard deviation based) +- Use `'trendCiPercentile'` if extremes show different variability pattern +- Add seasonality (`'seasonal'`) if clear annual cycle +- Use `'trendLinear'` for simple monotonic trends + +**GEV vs GPD**: +- GEV: Theoretically sound for annual maxima, but fewer data points +- GPD: Uses more data (all threshold exceedances), but requires threshold selection +- Both: Run both and compare for robustness + +### 5. Diagnostic Checks + +Always examine: + +1. **Transformation plot**: Transformed series should be stationary +2. **Shape parameter**: Check if reasonable (negative → bounded, positive → heavy tails) +3. **P-values**: Assess statistical significance of trends +4. **Confidence intervals**: Wide intervals indicate high uncertainty +5. **Sensitivity**: Test different parameters to assess robustness + +### 6. Common Pitfalls + +❌ **Avoid**: +- Time window smaller than parameter of interest (e.g., 5-year window for 100-year return level) +- Ignoring seasonality in seasonal data +- Over-interpreting results with wide confidence intervals +- Extrapolating far beyond observed data range +- Using GEV for threshold exceedance or GPD for annual maxima + +✅ **Best Practices**: +- Plot diagnostics before believing results +- Compare multiple approaches (GEV vs GPD, std dev vs percentile) +- Document all parameter choices and reasoning +- Report confidence intervals, not just point estimates +- Validate with physical understanding of the phenomenon + +--- + +## Summary Table of Examples + +| Example | Transformation | Distribution | Data Type | Key Learning | +|---------|---------------|--------------|-----------|--------------| +| 1. Basic Non-Stationary | Trend + Seasonal | GEV + GPD | Water levels | Complete workflow, both approaches | +| 2. Percentile CI | Trend/Seasonal CI Percentile | GEV + GPD | Water levels | Alternative variability estimation | +| 3. Compare CI | Multiple methods | GEV + GPD | Wave heights | Systematic comparison of methods | +| 4. SPI/GPD | Trend CI Percentile | GPD only | Drought index | Long peaks, GPD-only approach | +| 5. TAS/GEV | Trend | GEV only | Temperature | Annual maxima, GEV-only approach | +| 6. Stationary | None (stationary) | GEV + GPD + Gumbel | Water levels | Traditional EVA, no time-variation | +| 7. Linear Trend | Trend Linear | GEV + GPD | Various | Linear trend assumption, multiple thresholds | +| 8. Negative Shape | Trend Linear | GPD with negative ε | Various | Bounded distributions | +| 9. Gumbel SPI | Trend CI Percentile | Gumbel (GEV) | Drought index | Simplified GEV model | + +--- + +## Next Steps + +After working through these examples: + +1. **Apply to your data**: Adapt the examples to your specific dataset +2. **Explore parameter sensitivity**: Test different time windows, peak distances, percentiles +3. **Compare methods**: Run multiple approaches and understand differences +4. **Validate results**: Use physical reasoning and domain knowledge +5. **Read copula examples**: For multivariate extremes, see copula documentation +6. **Consult function reference**: See `2_Function_Reference.md` for detailed parameter descriptions + +--- + +## References + +For theoretical background on these methods, see: +- `1_Core_Methodology.md` - Mathematical foundations +- `3_Workflow_Patterns.md` - Standard analysis patterns +- Coles, S. (2001). *An Introduction to Statistical Modeling of Extreme Values*. Springer. +- Extreme value theory literature for GEV and GPD distributions + +--- + +*Document version: 1.0* +*Last updated: 2024* +*Part of tsEVA 2.0 documentation suite*