StormSim PST

StormSim: PST/SST Extreme Value Engine

Source files: StormSim_Library/Extreme_Value_Analysis/StormSim_PST.m (entry point) plus StormSim_Library/Extreme_Value_Analysis/PST/{StormSim_MRL.m, StormSim_POT.m, StormSim_PST_Fit.m, StormSim_PST_plot.m, KernReg_LocalMean.m, Monotonic_adjustment.m, ecdf_boot.m}. Sole in-repo caller: StormSim_Library/Extreme_Value_Analysis/call_hazard_curve_builder.m (XC branch), itself called only from stormsim_pros.m (the extratropical branch of PROS). Workflow context: config.workflow = 1/4 (PROS/PROS-FB); gated by config.storm_sampling containing 'XC' (or 'CC', which combines this engine's output with JPM's — see StormSim JPM and StormSim PROS).

What PST is

Despite its two conflicting names — "Probabilistic Simulation Technique (PST)" in its own header, "Stochastic Simulation Technique (SST)" everywhere it's called from — this engine is not a stochastic/Monte Carlo storm-resampling method at all. It is a parametric-bootstrap Peaks-Over-Threshold (POT) / Generalized Pareto Distribution (GPD) extreme value analysis: a standard, well-established extreme value theory (EVT) technique for fitting a distribution's tail from a sample of storm peaks, distinct in kind from a Monte Carlo storm-simulation method like LCS's actual Poisson-arrival storm sampling.

Concretely: it takes the per-storm response population that PROS's response-based (RB) machinery already computed from the full extratropical (XC) storm suite, treats those values as a Peaks-Over-Threshold sample, and fits a GPD tail to it — using a bootstrap over that same sample (resample-with-jitter, refit, repeat) to produce confidence limits. Nothing about the fit itself resamples or simulates new storms; the only resampling in play is the statistical bootstrap used to quantify the uncertainty of the fitted curve. This page corrects the framing used elsewhere on the wiki, which (before this page was written) described PST/SST as an external, resampling/simulation-style routine akin to a stochastic storm simulator — it is neither external to this codebase nor conceptually a storm simulator.

Theory: POT/GPD extreme value analysis

Classical block-maxima EVT (fitting a Generalized Extreme Value distribution to, say, one peak per year) throws away most of a record — only the single largest value per block is used. Peaks-Over-Threshold analysis instead keeps every storm peak that exceeds some chosen threshold, and a well-known EVT result (the Pickands–Balkema–de Haan theorem) says that, for a suitably high threshold, the distribution of excesses above that threshold converges to a Generalized Pareto Distribution (GPD). This is the standard basis for building a hazard curve's tail — the rare, high-magnitude end where a raw empirical record is thin — from a POT sample.

The threshold choice is the crux of a POT analysis, and it is a real trade-off, not a formality: set it too high and too few excesses remain to fit a stable GPD; set it too low and the excesses no longer behave like a GPD tail (the asymptotic justification breaks down), risking a biased fit. Mean Residual Life (MRL) analysis is the standard diagnostic for this: the mean of the excesses above a candidate threshold should be approximately linear in the threshold if the GPD approximation holds, so scanning candidate thresholds and looking for where that linearity holds (or where a fit-error measure is minimized) gives an objective way to pick a defensible threshold instead of eyeballing a plot.

Once a threshold and GPD fit are in hand, the resulting hazard curve is a hybrid: the GPD-fitted quantiles describe the tail above the threshold (where extrapolation into rare, unobserved severities is needed), while the empirical points describe the body below the threshold (where the observed record already provides a reasonable direct estimate) — the curve is not a single parametric fit over the whole range, but a merge of the two.

Because any single GPD fit to a finite sample carries sampling uncertainty, StormSim wraps the whole fit — threshold selection included in spirit, though in practice the threshold criterion itself is fixed — in a parametric bootstrap: resample the empirical peaks (with replacement, jittered — see the algorithm walkthrough below) many times, refit the whole pipeline per replicate, and read the spread across replicates as the confidence-limit band. This is the same general logic behind the confidence limits described on StormSim PROS's "Hazard curves: the central object" section — it just happens, for the XC/PST engine specifically, that the CL spread comes from a bootstrap resample of the observed peaks rather than from a closed-form/replicate-based uncertainty propagation like the "444 replicates" mechanism PROS uses for TC forcing uncertainty.

Algorithm walkthrough

1. Entry. stormsim_pros.m computes per-storm structure/forcing responses, then, whenever config.storm_sampling contains 'XC'/'CC', hands them to call_hazard_curve_builder(config, [], Resp, 'XC', use_aep, outPath).
2. Response packaging. call_hazard_curve_builder.m builds response_data with DataType hardcoded to 'POT' — the per-storm response values are treated as one discrete event-peak per XC storm already (declustering already happened by construction of the synthetic XC storm suite upstream), not as a continuous record needing inter-event-time declustering. Nyrs = config.Nyrs_XC*dcols.
3. Tuning knobs (eva_options) are set here, and every one of them is hardcoded, none config-driven: ind_Skew=0, use_AEP=config.pros_use_aep, prc (parsed from config.project_CLs), stat_print=0, tLag=0, GPD_TH_crit=2 (the WMSE/minimum-error MRL criterion — see step 6), apply_GPD_to_SS=1, bootstrap_sims=100.
4. Input handling inside StormSim_PST.m. Validates/normalizes inputs, strips flagged/NaN/Inf/≤0 values, and aborts if fewer than 4 unique values remain. For DataType='POT' (the only path exercised in production), the incoming data is used directly as the POT sample — StormSim_POT.m's inter-event-time declustering logic is present and complete but never actually invoked on this path (it would only run for DataType='Timeseries', which has no current caller in this codebase). Duplicate response values are given a tiny (1e-6) jitter to avoid downstream tie issues.
5. Fit step — StormSim_PST_Fit.m, the real core:
   • Sets up fixed AEP/AEF return-period grids and a dense log-spaced plotting grid.
   • Builds the empirical CDF via Weibull plotting position (P = m/(n+1)), then applies a Partial Duration Series (PDS) Lambda correction — Lambda_hist = Nstorms/Nyrs — that converts the rank-based exceedance probability into an annual rate, which is what a hazard curve's AEP/AEF axis actually needs.
   • Bootstrap (ecdf_boot.m). rng('default') fixes a deterministic seed, then for each of bootstrap_sims = 100 replicates the empirical peaks are resampled with replacement and given Gaussian jitter sized to the local spacing between adjacent order statistics — not a plain resample-and-refit. This bootstrap is what ultimately produces the confidence-limit spread.
   • GPD applicability gate. The nominal rule is "fit GPD only if the POT sample has ≥20 events and the record is ≥20 years, otherwise fall back to empirical-only," but apply_GPD_to_SS = 1 is hardcoded, so the gate is effectively always satisfied in production regardless of sample size.
   • Threshold selection via Mean Residual Life — StormSim_MRL.m. For each candidate threshold, computes the mean excess and a weighted linear regression of mean-excess-vs-threshold; the weighted mean-square error (WMSE) of that fit is denoised via Nadaraya-Watson kernel regression (KernReg_LocalMean.m) to locate its minimum. This denoised WMSE-minimum ("CritWMSE") is the criterion actually used in production, since GPD_TH_crit = 2 is hardcoded. A second criterion (GPD_TH_crit = 1, which instead picks the threshold giving a sample intensity closest to 2 events/year) exists in the code but is never selected on the live pipeline. If MRL fails to find a valid minimum for a given bootstrap replicate, the threshold falls back to 0.99 × min(that replicate's bootstrap sample).
   • GPD fit per bootstrap replicate. fitdist(exceedances, 'GeneralizedPareto', 'theta', threshold) is refit independently for each of the 100 bootstrap replicates — this per-replicate refit (not a single point-estimate fit with an assumed sampling distribution) is what makes the resulting spread a genuine resampling-based uncertainty quantification. The fitted shape parameter is clamped to [-0.5, 0.3] ("Limits determined by NCNC" per the code comment, referring to author Norberto C. Nadal-Caraballo).
   • Hazard curve construction. For each bootstrap replicate, GPD-fitted quantiles above the threshold are merged with empirical points below the threshold into one hybrid curve, then interpolated onto the fixed plotting grid in log-space.
   • Percentile summary. The mean across all 100 bootstrap-replicate curves becomes the "central"/labeled-"50" curve — this is a literal bootstrap mean, not an actual 50th-percentile/median computed via a percentile function. This is a real labeling nuance: don't assume the "50" curve is a median. prctile across the 100 replicates at the user-requested config.project_CLs values gives the CL bands.
   • Monotonic adjustment (Monotonic_adjustment.m). The raw bootstrap-mean/percentile curves can come out non-monotonic ("jumpy") when the MRL-selected threshold is too low. This function works in log-x space, finds any point where the response value increases as frequency decreases (physically backwards for a hazard curve), and replaces the offending point via linear extrapolation from the two prior slopes. It is a local patch heuristic, not a formal isotonic regression — describe it as such rather than as a rigorous statistical smoother.
   • Final interpolation onto fixed ARI/response-magnitude tables for downstream use, plus a soft (non-throwing) sanity-check flag if the mean curve's value at 0.1 AEP/AEF exceeds 1.75× the empirical value there.
6. Skew-tide augmentation path (the ind_Skew = 1 branch, present throughout StormSim_PST_Fit.m) is real, complete code that adds a predicted skew-tide component (via a Gaussian process metamodel) to the bootstrapped surge sample before building the curve — but it is never exercised in production, since call_hazard_curve_builder.m always hardcodes ind_Skew = 0. Present-but-dormant, same status as the POT-declustering note in step 4.

Config fields

Field	Effect	Documented in config_reference.md?
config.Nyrs_XC	Record length used for the PDS Lambda correction and the MRL sample-intensity rate	Yes — §10, "Set by call_input_parser.m"
config.pros_use_aep	AEP vs. AEF convention for the hazard curve axis	Yes — §9, PROS Module
config.project_CLs	Which percentiles get computed across the bootstrap (max 4, parsed from a string expression)	Yes — §1, Project & Output Settings
config.storm_sampling	Gates whether the XC/PST branch runs at all ('XC' or 'CC')	Yes — §1, Project & Output Settings

Not configurable: the algorithm's actual tuning knobs — bootstrap replicate count (bootstrap_sims = 100), the GPD threshold criterion (GPD_TH_crit = 2), the skew-tide toggle (ind_Skew = 0), and the inter-event lag (tLag = 0) — are all hardcoded in call_hazard_curve_builder.m's eva_options struct, not exposed via config or the Excel input template. Anyone wondering why PST behavior can't be tuned from the input schema: it currently can't be, short of editing call_hazard_curve_builder.m directly.

Implementation reference

This section is a pointer for readers who want to trace the code directly; it is intentionally not the focus of this page.

• Signature: [SST_output] = StormSim_PST(response_data, pst_options, plot_options). Full field contracts for all three input structs are documented in the file's own header comment block.
• Success-path output fields: staID, RL, MRL_output, HC_plt, HC_tbl, HC_tbl_rsp_x, HC_emp, HC_tbl_rsp_y, HC_plt_x, HC_tbl_x. Note: POT and Warning fields are preallocated on the struct but are only populated in the catch/error branch — on a successful run, callers must not assume they exist.
• Caller-side mapping. call_hazard_curve_builder.m maps SST_output onto the shared Output.x_plot/y_plot/x_table/y_table/x_table_ARI/tbl_rsp_x/tbl_rsp_y/CL struct (the same contract StormSim PROS describes hazard curves using), with Output.CL = [50, prc] — the "50" here is a label for the bootstrap mean, not a literal 50th percentile (see above). Output.POT in that struct is actually the caller's own pre-PST raw per-storm response data (data_to_process.(resp_var)), not SST_output.POT (which, again, doesn't exist on a successful run) — a coincidental field-name reuse between the caller's struct and PST's own internal (error-only) field, worth a one-line caveat if tracing this code.
• Entry chain: stormsim_pros.m (XC/CC branch) → call_hazard_curve_builder.m (storm_type = 'XC' case) → StormSim_PST.m → StormSim_PST_Fit.m (core fit) → StormSim_MRL.m + KernReg_LocalMean.m (threshold selection) + ecdf_boot.m (bootstrap) + Monotonic_adjustment.m (post-fit patch) → StormSim_PST_plot.m (plotting, gated by plot_options.create_plots, which call_hazard_curve_builder.m hardcodes to 0 on the production path).
• Two known-dormant latent issues (informational, not action items):
  • StormSim_PST.m references a bare path_out variable in a disp() call that is never assigned in that scope (it should read plot_options.path_out) — but this line only executes when pst_options.stat_print = 1, which is always hardcoded to 0, so it is currently unreachable/harmless.
  • StormSim_POT.m (the declustering routine) and the skew-tide branch inside StormSim_PST_Fit.m are both complete, real code that simply never executes on the current production path — treat them as dormant, not as "the" active step, when reading the source.

Naming and terminology notes

• PST vs. SST — a real, unresolved inconsistency, not a typo. StormSim_PST.m's own header calls this the "Probabilistic Simulation Technique (PST)." The very same file's runtime banner and save-name strings call it the "StormSim-SST Tool." StormSim_MRL.m's header calls it "StormSim-SST-Fit (Statistics)." StormSim_PST_plot.m's titles/filenames say "StormSim-SST." The caller, call_hazard_curve_builder.m, labels it the "Stochastic Simulation Technique (SST)" in its progress messages — this is likely the source of the "resampling/simulation" framing that appeared elsewhere on this wiki before this page was written. No file in the codebase reconciles these two different spelled-out names ("Probabilistic Simulation Technique" vs. "Stochastic Simulation Technique"); this is reported here as an open inconsistency, not something resolved by picking one as authoritative.
• Alpha status. Every file in the PST family is explicitly marked "ALPHA VERSION — FOR TESTING/INTERNAL TESTING ONLY" in its own header (StormSim_PST.m, StormSim_MRL.m, StormSim_POT.m, ecdf_boot.m). This is a real, current disclaimer in the source — worth surfacing plainly rather than omitting — even though this is, in practice, live production code sitting in the middle of the PROS XC hazard-curve pipeline.
• Author credits. StormSim_MRL.m, StormSim_POT.m, and ecdf_boot.m credit Norberto C. Nadal-Caraballo, PhD and Efrain Ramos-Santiago as authors. KernReg_LocalMean.m and Monotonic_adjustment.m credit E. Ramos-Santiago alone.

See also

• StormSim PROS — the workflow that calls this engine for extratropical (XC) storms and combines its output with JPM's for combined (CC) hazard curves.
• StormSim JPM — the sibling extreme-value engine PROS uses for tropical (TC) storms, based on the Joint Probability Method rather than POT/GPD.