Topological Conditioning via Natural Visibility Graphs for Monte Carlo Simulation of Power Prices

1. Introduction

1.1. Stylized Facts of Electricity Price Dynamics

Electricity prices exhibit a distinctive combination of statistical properties that collectively challenge standard modelling frameworks [1]. Among energy commodities, power markets are distinguished by the physical non-storability of electricity: supply and demand must be balanced instantaneously, so any transient mismatch between generation capacity and consumption produces price spikes that can reach values orders of magnitude above normal levels within minutes, only to revert once equilibrium is restored [2,3]. This mechanism generates time series with mean reversion, heavy-tailed, skewed marginal distributions and pronounced volatility clustering: periods of elevated price variability tend to persist over time, alternating with calmer regimes [4,5].

A distinctive feature of electricity prices is mean reversion: prices tend to revert toward a long-run equilibrium level at a speed that varies across markets and time periods. The AR(1) model

$$X_t=(1-\alpha )X_{t-1}+{\epsilon }_t,$$

(1)

where $X_t$ denotes the detrended log-price and $\alpha$ governs the speed of mean reversion, provides a natural decomposition of the price dynamics into a persistent component and a residual innovation sequence $\left\{{\epsilon }_t\right\}$. This decomposition motivates a two-stage modelling approach [6]: the mean-reversion parameter $\alpha$ is estimated in the first stage via a distribution-free neural estimator that requires no assumption on the innovation distribution; in the second stage, the statistical properties of the innovations — fat tails, asymmetry, and volatility clustering — are addressed. The two stages can be treated independently, avoiding the distributional misspecification that arises when both are estimated jointly under a single parametric model.

Three interrelated stylized facts define the statistical landscape of electricity innovations:

1.

Fat tails and asymmetry. The distribution of log-price increments departs substantially from normality, exhibiting positive excess kurtosis and non-zero skewness driven by price spikes.

2.

Volatility clustering. Large innovations tend to cluster in time. The autocorrelation function of squared innovations decays slowly from significantly positive values at short lags, confirming that the conditional variance is persistent.

3.

Non-stationarity. Long-term trends attributable to structural changes, seasonality, and the energy transition require explicit removal before statistical analysis of short-term dynamics.

Any simulation framework aiming to generate realistic electricity price scenarios must jointly reproduce these three features. Failure to do so compromises downstream applications including risk management, derivative pricing, and investment analysis under uncertainty [7,8,9].

1.2. Limitations of Parametric Approaches

The standard toolkit for modelling financial volatility dynamics consists of parametric models, each of which imposes a specific functional form on the conditional distribution of innovations. While tractable and widely used, these models share structural limitations when applied to electricity markets.

GARCH models [10,11] are the workhorse of financial volatility modelling. In their standard formulation they assume symmetric, zero-skewness innovations, which is a structural limitation for electricity markets where the empirical distribution of innovations is markedly asymmetric. Extensions of the basic GARCH framework can accommodate skewness, heavy tails and other non-Gaussian features, but at the cost of substantially more complex parameterisations that are harder to estimate and interpret, and that still impose a single, globally fixed distributional shape across all market regimes.

Jump-diffusion models [3,12] augment a continuous diffusion with a Poisson jump component, explicitly representing price spikes as discrete events. Although conceptually appealing, these models require estimation of jump intensity and jump size distribution from limited historical data. In markets where spikes are rare but extreme, maximum likelihood estimation of jump parameters is notoriously unstable, and out-of-sample performance degrades rapidly when the historical spike record is short or non-representative of future dynamics.

Regime-switching models [13] assume that the price process alternates among a finite number K of latent states, each characterised by its own conditional distribution. Transition probabilities are estimated from data via the EM algorithm or Bayesian inference. These models can capture structural changes in mean and variance, but they impose several constraints: the number of regimes must be pre-specified; the within-regime distribution remains fully parametric; and the state-space representation forces a discrete partition of a volatility landscape that is inherently continuous.

The common thread across all parametric approaches is deeper than the specific limitations noted above. Every parametric model must solve two problems simultaneously: specifying the marginal distribution of innovations and specifying the mechanism that generates temporal dependence. These two problems are structurally coupled in the parametric framework — the GARCH variance recursion determines both the clustering and, through the chosen innovation distribution, the tail behaviour — so that improving one aspect of the fit (e.g., adding skewness via a skewed-t distribution, or richer dynamics via an EGARCH specification) necessarily re-parameterises the other. Distributional misspecification at any level propagates systematically through Monte Carlo simulations.

These limitations call for an approach that imposes no global distributional assumption on the innovations and conditions the temporal dependence on the data structure itself, rather than on a pre-specified parametric model. For electricity markets — where the marginal distribution is heavy-tailed, asymmetric, and market-specific — such an approach eliminates the model-selection step that any parametric framework requires.

1.3. Visibility Graphs as a Non-Parametric Tool

Visibility graphs [14] offer precisely such a mechanism: a principled way to encode the temporal structure of a time series into a network representation without any distributional assumption. Given a time series ${\left\{{\epsilon }_t\right\}}_{t=1}^T$, two observations at times $i<j$ are connected by an edge in the Natural Visibility Graph (NVG) if and only if all intermediate observations ${\epsilon }_k$ ($i<k<j$) satisfy:

$${\epsilon }_k<{\epsilon }_i+({\epsilon }_j-{\epsilon }_i)·{\displaystyle \frac{k-i}{j-i}}.$$

(2)

The right-hand side is the linear interpolant between $(i,{\epsilon }_i)$ and $(j,{\epsilon }_j)$ evaluated at time k: two points are mutually visible if no intermediate observation rises above the straight line connecting them.

The NVG encodes rich local information about the time series: a high observation connects to many neighbours and spans long temporal distances, while a low observation is shielded by surrounding values and has few, short-range connections. This topological structure is sensitive to the local volatility regime — spikes, troughs, and periods of low variability each produce characteristic connectivity patterns — without requiring any parametric model for that regime.

1.4. Topological Conditioning for Monte Carlo Simulation

Previous work has shown that visibility graph topology carries statistically significant information about the underlying stochastic process [15,16] and has been applied to financial and energy time series to characterise distributional properties and detect structural changes [17,18,19,20,21]. While previous studies have used visibility graph topology as a descriptive or diagnostic tool, the present paper is the first to exploit it as a conditioning mechanism for Monte Carlo simulation: the topological structure of the innovation sequence directly governs the local distributional regime at each simulation step, without any parametric assumption on regime dynamics or transitions. Specifically, given the mean-reversion parameter $\widehat{\alpha }$ estimated in the first stage [6] and the resulting innovations

$${\epsilon }_t=X_t-(1-\widehat{\alpha })X_{t-1},$$

(3)

the NVG built on the innovation sequence $\left\{{\epsilon }_t\right\}$ is used to encode the local market state into a feature vector that conditions a locally parametric simulation engine.

Three backward-looking features are extracted from the NVG: the backward degree $d_t^{-}$ (number of edges connecting the current node to past nodes), the maximum backward edge length ${\ell }_t^{-}$ (temporal span of the longest such edge), and the squared current innovation ${\epsilon }_t^2$ (a dynamic proxy for the local volatility level). Together they form the conditioning vector ${\mathbf{f}}_t=({\epsilon }_t^2,d_t^{-},{\ell }_t^{-})$ that encodes the local volatility regime through the geometry of the time series itself, with no global assumption on the form of the innovation distribution.

The conditioning mechanism works as follows. At each point in the historical record, the feature vector ${\mathbf{f}}_t$ characterises the local market state: the squared innovation captures the current volatility level, while the two topological quantities describe the geometric structure of the recent price trajectory. States with similar feature vectors are interpreted as instances of the same volatility regime, regardless of when they occurred. The k historical states most similar to the current one are identified by nearest-neighbour search in the feature space, and the innovations that followed those states are used to fit a local parametric distribution. As the local working model, a Jones–Pewsey Sinh-ArcSinh (SAS) distribution [23] is adopted: a flexible family that simultaneously accommodates heavy tails, positive or negative asymmetry, and a wide range of kurtosis values. Its scale, skewness, and tail weight are estimated locally for each historical state, so that a single conditioning operation jointly determines all three distributional properties.

This mechanism extends naturally to Monte Carlo simulation: starting from the last observed state, the feature vector is updated at each step using the most recent simulated innovations, a new set of nearest neighbours is identified in the historical record, and the next innovation is drawn from the corresponding local SAS distribution. The result is a collection of synthetic price paths that inherit the distributional and temporal structure of the observed series — fat tails, asymmetry, and volatility clustering — without any global assumption on the innovation distribution. The simulated innovations are then fed back into the AR(1) structure to generate synthetic price paths directly applicable to risk quantification and derivative pricing.

The framework is validated on two electricity markets with contrasting distributional characteristics: the Italian Power Exchange (negative skewness, moderate kurtosis) and PJM West Hub (strong positive skewness, heavy tails). With a single neighbourhood size $k=10$ and no market-specific re-calibration, the framework achieves simultaneous coverage of $\widehat{\sigma }$, $\widehat{\gamma }$, $\widehat{\kappa }$, and ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ for both markets, assessed via a coverage test on $N=500$ Monte Carlo simulated paths. The neighbourhood size k is the sole tunable parameter of the framework and serves as a natural adaptation lever for markets with substantially different distributional characteristics.

The remainder of the paper is organised as follows. Section 2 describes the data and preprocessing pipeline. Section 3 details the methodology: NVG construction, feature extraction, and local SAS fit. Section 4 presents the empirical results: coverage analysis ($\widehat{\sigma }$, $\widehat{\gamma }$, $\widehat{\kappa }$, ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$) and ACF reproduction for both markets. Section 5 discusses the single-mechanism nature of the topological conditioning, the locally parametric and distribution-free character of the framework, and its market-specific interpretation. Section 6 concludes and outlines limitations and directions for future research. The overall pipeline is summarised in Figure 1.

2. Data and Preprocessing

2.1. Markets and Data

The empirical analysis is conducted on two electricity spot markets with contrasting geographic and structural characteristics.

PUN (Prezzo Unico Nazionale) is the reference price of the Italian day-ahead power exchange, operated by GME (Gestore dei Mercati Energetici SpA). Daily price observations are used, yielding $T=1825$ AR(1) innovations after the first lag is consumed. PUN prices reflect Italy’s dependence on imported gas and hydropower, and exhibit recurrent negative spikes associated with periods of excess renewable generation.

PJM West Hub is the day-ahead locational marginal price at the PJM Western Hub, the largest organised electricity market in the United States. Business day price observations are used (weekends and US market holidays excluded), yielding $T=1224$ AR(1) innovations after the first lag is consumed. PJM prices are driven by temperature-sensitive demand, thermal generation costs, and transmission constraints, and exhibit recurrent upward spikes associated with cold or heat waves.

Both datasets cover the same five-year window, 1 January 2019 – 31 December 2023, providing a consistent temporal basis for comparison. This period encompasses major common shocks: the COVID-19 demand disruption (2020), the post-pandemic recovery, and the extreme price volatility triggered by the 2022 energy crisis following the Russia–Ukraine conflict. PUN is characterised by moderate kurtosis and negative skewness, while PJM exhibits higher kurtosis and pronounced positive asymmetry (Table 1).

2.2. Preprocessing Pipeline

Raw price series $P_t$ are transformed into innovations ${\epsilon }_t$ through three sequential steps, illustrated in Figure 1.

Step 1 — Logarithmic transformation. Log-prices ${\tilde{P}}_t=ln{P}_t$ are used, which stabilise variance and convert multiplicative price dynamics into additive form.

Step 2 — LOESS detrending. A nonparametric LOESS smoother [24,25] with bandwidth $h=0.10$, selected by leave-one-out cross-validation over the grid $\{0.05,0.10,0.20\}$, is fitted to ${\tilde{P}}_t$, yielding the smooth trend ${\widehat{m}}_t$. The detrended series is $X_t={\tilde{P}}_t-{\widehat{m}}_t$. LOESS is preferred over parametric detrending because it imposes no functional form on the long-term trend and is robust to outliers.

Step 3 — Innovation extraction. The price process $X_t$ is modelled as an AR(1) mean-reverting process:

$$X_t=(1-\alpha )X_{t-1}+{\epsilon }_t,\mathrm{Cov}({\epsilon }_t,{\epsilon }_s)=0(t\ne s),$$

(4)

where $\alpha$ is the mean-reversion speed; the condition $|1-\alpha |<1$ ensures stationarity of $\left\{X_t\right\}$. No further assumption is imposed on ${\epsilon }_t$: in particular, innovations may be heteroscedastic (as in GARCH-type processes), and the Yule–Walker identification of $\alpha$ remains valid under zero serial correlation alone. Given an estimate $\widehat{\alpha }$, innovations are recovered as ${\epsilon }_t=X_t-(1-\widehat{\alpha })X_{t-1}$.

The mean-reversion parameter $\alpha$ is estimated by a Temporal Convolutional Network (TCN) [22] trained on simulated AR(1) paths with diverse innovation distributions [6]. Unlike classical estimators (OLS, Yule–Walker, GARCH-based), which condition on a specific innovation model and therefore yield different estimates of $\alpha$ on the same sample depending on the assumed distribution, the TCN learns to read the autocorrelation structure of the series directly, without any distributional conditioning. The values $\widehat{\alpha }=0.2857$ (PUN) and $\widehat{\alpha }=0.3540$ (PJM) are taken directly from [6].

2.3. Descriptive Statistics

Table 1 reports the main descriptive statistics of the innovation sequences. Both markets exhibit substantial departures from normality: Pearson kurtosis of $5.48$ (PUN) and $9.90$ (PJM), well above the Gaussian benchmark of 3, and significant asymmetry. The autocorrelation of squared innovations at lag 1, ${\rho }_1\left({\epsilon }_t^2\right)$, is positive and statistically significant in both cases, confirming the presence of volatility clustering. Figure 2 shows the innovation time series and empirical densities alongside a Gaussian reference.

3. Methodology

3.1. Natural Visibility Graph

Given the innovation sequence ${\left\{{\epsilon }_t\right\}}_{t=1}^T$, the Natural Visibility Graph (NVG) is the undirected graph $G=(V,E)$ where $V=\{1,\dots ,T\}$ is the set of nodes, each representing one observation in the time series, and $E\subseteq V\times V$ is the set of edges defined by the visibility condition: $(i,j)\in E$ if and only if $i\ne j$ and condition (2) holds for the pair $(i,j)$ with $i<j$. Each edge $(i,j)\in E$ therefore represents a direct line-of-sight between observations ${\epsilon }_i$ and ${\epsilon }_j$, with no intermediate observation rising above the straight line connecting them. The NVG is constructed once on the full historical sequence using an $O\left(T^2\right)$ algorithm [14,26].

The topological properties of the NVG encode local dynamical information about the time series. A node i with a large value ${\epsilon }_i$ dominates its neighbourhood: it is visible from many distant nodes, yielding high degree and long edges. Conversely, a node in a low-volatility trough is partially occluded by surrounding observations and has short, local connections. This sensitivity to the local volatility regime is the key property exploited in the simulation framework.

Figure 3 illustrates the geometric structure of the NVG on a representative window of 80 innovations for each market. The left panels display the arc representation: each arc connects a pair of observations satisfying condition (2) (arcs are drawn as curves for visual clarity; the underlying criterion is based on straight lines, as defined in (2)), with the arc colour encoding the temporal distance between the two nodes. High-valued observations (spikes) spawn long-range arcs spanning the entire window, while low-valued observations generate only short, local connections. The right panels show the same graph in a Kamada–Kawai network layout: node size is proportional to degree and node colour encodes the innovation value, revealing how spike nodes (large, red) occupy structurally central positions in the network while trough nodes (small, blue) cluster at the periphery.

3.2. Backward Topological Features

At any simulation step t, the future is not yet observed: only ${\epsilon }_1,\dots ,{\epsilon }_t$ are available. Computing the NVG on the full history would require $O\left(T^2\right)$ operations and would yield features whose meaning changes as T grows. The graph is therefore restricted to a rolling window of fixed length $W=100$ observations ending at t:

$${\mathcal{W}}_t=\{{\epsilon }_{t-W+1},\dots ,{\epsilon }_t\},$$

(5)

and build a local NVG $G_t=(V_t,E_t)$ on ${\mathcal{W}}_t$ at each step (for $t<W$, the window is zero-padded on the left to maintain fixed length). This choice serves two purposes: (i) it enforces causality, since ${\mathcal{W}}_t$ uses only observations available at time t; (ii) it ensures temporal locality, so that the topological features reflect the current market regime rather than a long-run average over the entire series.

Attention is restricted to backward edges of $G_t$ — edges connecting the last node of the window to all earlier nodes it can see — and the following causal feature vector is defined:

$${\mathbf{f}}_t=\left({\epsilon }_t^2,d_t^{-},{\ell }_t^{-}\right)\in {\mathbb{R}}^3,$$

(6)

where

$$\begin{array}{cc}\hfill d_t^{-}& =\sum _{\begin{array}{c}j\in {\mathcal{W}}_t\\ j<t\end{array}}\mathbf{1}\left[(j,t)\in E_t\right]\hfill \end{array}$$

(7)

is the backward degree — the number of past observations in ${\mathcal{W}}_t$ that are visible from t — and

$$\begin{array}{cc}\hfill {\ell }_t^{-}& =max\{t-j:j\in {\mathcal{W}}_t,j<t,(j,t)\in E_t\}\hfill \end{array}$$

(8)

is the maximum backward edge length — the temporal distance to the farthest visible predecessor within the window.

The three components carry complementary information. It is worth noting that ${\epsilon }_t^2$ is not a topological property of the NVG: it is the observed value of the current innovation, included in the feature vector to anchor the conditioning to the absolute volatility level. The two genuinely topological components $d_t^{-}$ and ${\ell }_t^{-}$ capture the relative structure of the local window, which may not fully identify the volatility level on its own. ${\epsilon }_t^2$ is preferred over the signed innovation ${\epsilon }_t$: conditioning on the signed value would make the drawn distribution depend on the direction of the innovation, creating a systematic link between successive innovations and inducing $\mathrm{Cov}({\epsilon }_{t+1},{\epsilon }_t)\ne 0$, in violation of the AR(1) assumption. The squared value captures the volatility level symmetrically, avoiding this directional bias. $d_t^{-}$ measures how many past observations are visible from t: a high value indicates that ${\epsilon }_t$ lies above its recent neighbours (low-volatility regime, unobstructed line of sight to many predecessors), while a low value signals that ${\epsilon }_t$ is buried in a volatile cluster where nearby spikes block the view. ${\ell }_t^{-}$ captures the temporal horizon of visibility: it is large for isolated spikes (visible far back in time) and small for values embedded in a volatile cluster (visibility blocked by recent high-amplitude neighbours).

For computational efficiency, the backward features are computed via a vectorised algorithm that precomputes the NVG index arrays for a window of length W once, requiring $O\left(W^2\right)$ memory and reducing the per-step cost to $O\left(W\right)$ floating-point operations. Figure 4 illustrates the NVG structure and backward features on a representative window of PJM innovations: the current node (red bar) is connected to all past nodes visible from it, with backward edges highlighted in red. The annotation shows the resulting values of $d_t^{-}$ and ${\ell }_t^{-}$ for that node.

Before nearest-neighbour search, features are standardised:

$${\tilde{\mathbf{f}}}_t={\displaystyle \frac{{\mathbf{f}}_t-{\widehat{\mu }}_f}{{\widehat{\sigma }}_f}},$$

(9)

where ${\widehat{\mu }}_f=({\widehat{\mu }}_1,{\widehat{\mu }}_2,{\widehat{\mu }}_3)$ and ${\widehat{\sigma }}_f=({\widehat{\sigma }}_1,{\widehat{\sigma }}_2,{\widehat{\sigma }}_3)$ are computed component-wise over all $T-1$ historical nodes:

$${\widehat{\mu }}_j={\displaystyle \frac{1}{T-1}}\sum _{t=1}^{T-1}{f}_{t,j},{\widehat{\sigma }}_j=\sqrt{{\displaystyle \frac{1}{T-1}}{\sum }_{t=1}^{T-1}{\left(f_{t,j}-{\widehat{\mu }}_j\right)}^2}.$$

(10)

The division in (9) is performed component by component, so that each feature is mapped to approximate zero mean and unit variance before the Euclidean distance is computed. This step is necessary because the three components operate on very different numerical scales: ${\epsilon }_t^2$ is typically small (of the order of the innovation variance), while $d_t^{-}$ counts visible predecessors and ${\ell }_t^{-}$ measures temporal distances, both of which can be orders of magnitude larger. Without standardisation, the Euclidean distance in feature space would be dominated by ${\ell }_t^{-}$, making the nearest-neighbour search effectively blind to the volatility signal carried by ${\epsilon }_t^2$.

3.3. Local SAS Fit

Let $k\ge 1$ be a fixed integer (the neighbourhood size) and, for each historical node $i\in \{1,\dots ,T-1\}$, let ${\mathcal{N}}_k\left(i\right)$ denote the set of the k indices $j\in \{1,\dots ,T-1\}\setminus \left\{i\right\}$ with smallest Euclidean distance $\parallel {\tilde{\mathbf{f}}}_i-{\tilde{\mathbf{f}}}_j\parallel$ to node i in the standardised feature space (9). Note that the feature vectors ${\tilde{\mathbf{f}}}_i$ are each computed from a local rolling window of length W ending at i (Section 3.2), so they capture the local dynamical state at each point in time. The nearest-neighbour search, however, ranges over the entire historical record $\{1,\dots ,T-1\}$: it identifies which past episodes — wherever they occurred in the series — most closely resembled the current local state. The rolling window defines the geometry of each node; the k-NN search compares these geometries globally across all available history.

The successor of any node j is the immediately following innovation ${\epsilon }_{j+1}$, i.e. the observation recorded one step after j in the historical sequence. Rather than drawing directly from the observed successors of the k nearest neighbours — which would confine the simulation to the historical range — a Sinh-ArcSinh (SAS) distribution [23] is fitted by maximum likelihood on those successors, allowing the simulation to generate innovations beyond the observed sample.

The SAS Family

Given the successors $\{{\epsilon }_{j+1}:j\in {\mathcal{N}}_k\left(i\right)\}$ of the k nearest neighbours of node i, a SAS distribution is fitted by maximum likelihood. The family is parametrised as:

$$X\stackrel{d}{=}\lambda sinh\left({\displaystyle \frac{\mathrm{arcsinh}\left(Z\right)+\eta }{\delta }}\right),Z\sim \mathcal{N}(0,1),$$

(11)

where $\lambda >0$ is the scale parameter, $\eta \in \mathbb{R}$ controls skewness (positive $\eta$ induces right-skewed distributions), and $\delta >0$ controls tail weight (smaller $\delta$ yields heavier tails). The SAS family is flexible enough to capture the full range of skewness and kurtosis observed in power markets and encompasses symmetric unimodal distributions as a special case ($\eta =0$).

Before any simulation path is generated (offline phase), parameters $({\widehat{\lambda }}_i,{\widehat{\eta }}_i,{\widehat{\delta }}_i)$ are estimated by maximum likelihood (Nelder-Mead) for every historical node $i\in \{1,\dots ,T-1\}$: the SAS is fitted on the set of successors $\{{\epsilon }_{j+1}:j\in {\mathcal{N}}_k\left(i\right)\}$. Since all nodes and their successors belong to the historical record, no future observations are required at this stage.

The Unifying Principle

The critical property of this construction is that all three parameters — scale ${\lambda }_i$, skewness ${\eta }_i$, and tail weight ${\delta }_i$ — are functions of the topological neighbourhood of node i and therefore vary continuously with the local dynamical state. In high-volatility regimes (large ${\epsilon }_t^2$, high backward degree $d_t^{-}$, long maximum edge ${\ell }_t^{-}$), the k-NN identifies historically similar spike-time nodes: their successors have large variance, pronounced skewness, and heavy tails, yielding a fitted SAS with large ${\widehat{\lambda }}_i$, non-zero ${\widehat{\eta }}_i$, and small ${\widehat{\delta }}_i$. In low-volatility regimes, the SAS contracts toward a near-Gaussian shape. Scale, skewness, and tail weight are thus jointly determined by a single conditioning operation — the topological k-NN search — without any additional parameters or regime-switching equations.

Locally Parametric and Distribution-Free

The framework is locally parametric in the sense of [27]: within each topological neighbourhood a parametric family (SAS) is fitted as a working approximation to the local distribution, while globally no parametric form is imposed on the innovation process and the total parameter count grows with the sample size. This is the defining property that distinguishes locally parametric methods from both fully parametric models (fixed parameter count, global distributional assumption) and classical non-parametric methods (no local approximating family). The SAS family serves as a flexible working model; any universal approximator (Gaussian mixture, normal-inverse-Gaussian, generalised hyperbolic) could be substituted without altering the conditioning mechanism.

3.4. Monte Carlo Simulation

In the online phase, at each simulation step t, the current feature vector ${\tilde{\mathbf{f}}}_t$ acts as a real-time regime indicator: it jointly encodes the current volatility level (through ${\epsilon }_t^2$) and the local topological structure (through $d_t^{-}$ and ${\ell }_t^{-}$), and is matched against all precomputed historical feature vectors to identify $j^{*}$, the historical node whose feature fingerprint most closely resembles the current one. Proximity in this hybrid feature space is interpreted as similarity of market regime. The next innovation is then drawn from the SAS already fitted for $j^{*}$, which encodes how the process behaved after historically similar regimes. The simulation never requires a future successor of t: the entire conditional distribution is inherited from the closest historical analogue. The per-step cost reduces to a single nearest-neighbour lookup and one scalar draw.

Formally, the procedure at each step t is as follows:

1.

Compute ${\tilde{\mathbf{f}}}_t=({\tilde{\epsilon }}_t^2,{\tilde{d}}_t^{-},{\tilde{\ell }}_t^{-})$ from the current rolling window.

2.

Find the closest historical analogue

$$j^{*}=\underset{i\in \{1,\dots ,T-1\}}{arg\; min}\parallel {\tilde{\mathbf{f}}}_t-{\tilde{\mathbf{f}}}_i\parallel .$$

3.

Draw ${\epsilon }_{t+1}\sim \mathrm{SAS}({\widehat{\lambda }}_{j^{*}},{\widehat{\eta }}_{j^{*}},{\widehat{\delta }}_{j^{*}})$ using transformation (11).

4.

Update the detrended log-price: $X_{t+1}=(1-\widehat{\alpha })X_t+{\epsilon }_{t+1}$.

Iterating steps 1–4 over $t=T,T+1,\dots ,T+H$ produces a Monte Carlo trajectory $\left\{X_t^{\left(s\right)}\right\}$ of detrended log-prices of length H, where each $X_{t+1}$ is generated from ${\epsilon }_{t+1}$ via the AR(1) recursion in step 4. Repeating the procedure for $s=1,\dots ,N$ independent seeds yields an ensemble of N synthetic price paths, each initialised at the last observed state $X_T$. The loop implements a continuous regime-matching mechanism: at every step, the current topological fingerprint selects the historical market state that most resembles the present one, and the simulation inherits the distributional properties that followed that state in the past. Regime transitions are not governed by a parametric equation but emerge naturally from the geometry of the feature space.

Internal consistency. The model assumption $\mathrm{Cov}({\epsilon }_t,{\epsilon }_s)=0$ for $t\ne s$, required for Yule–Walker identification of $\alpha$, is preserved by the simulation. The SAS parameters $({\widehat{\lambda }}_{j^{*}},{\widehat{\eta }}_{j^{*}},{\widehat{\delta }}_{j^{*}})$ are estimated in the offline phase on historical innovations that satisfy this assumption. In the online phase, each draw is made from a fixed distribution indexed by $j^{*}$: the topological conditioning determines which distribution to use, but does not introduce a linear dependence between ${\epsilon }_{t+1}$ and ${\epsilon }_t$. Serial dependence in the simulated sequence is therefore confined to the second moment (volatility clustering through ${\epsilon }_t^2$), which is precisely the structure targeted by the conditioning mechanism. The use of ${\epsilon }_t^2$ rather than the signed innovation ${\epsilon }_t$ in the feature vector is essential to this argument: conditioning on the signed value would make the drawn distribution depend on the direction of the innovation, creating a systematic link between successive innovations and inducing $\mathrm{Cov}({\epsilon }_{t+1},{\epsilon }_t)\ne 0$, in violation of $\mathrm{Cov}({\epsilon }_t,{\epsilon }_s)=0$.

4. Results

4.1. Coverage Analysis

Throughout this section, the following statistics are used to characterise the innovation distribution: $\widehat{\sigma }$ (standard deviation), $\widehat{\gamma }$ (skewness), $\widehat{\kappa }$ (Pearson kurtosis, Gaussian $=3$), and ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ (first-order autocorrelation of squared innovations). The last quantity captures volatility clustering and is not a distributional moment, but an autocorrelation measure of the squared process.

Table 2 reports the main coverage results for both markets. For each statistic, the empirically observed value is compared to the median and $[5\mathrm{th},95\mathrm{th}]$ percentile band across $N=500$ simulated paths of length T. All three statistics and ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ are covered simultaneously for both markets with a single choice of $k=10$ and no market-specific re-calibration.

The coverage results warrant several observations. First, the observed values generally lie near the lower portion of the $[p_5,p_{95}]$ band for skewness and the median for kurtosis, indicating that the simulated distribution of these moments is wide. The width reflects genuine sampling variability: for a series of length $T\le 1825$, a single observed path is one realisation of a stochastic process whose moment estimates carry substantial uncertainty. The topological SAS framework correctly reproduces this uncertainty rather than artificially concentrating simulated moments around a point estimate.

Second, the coverage of ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ — the first-order autocorrelation of squared innovations, the direct measure of volatility clustering — is achieved without any explicit model for the conditional variance. The clustering emerges endogenously from the topological conditioning: in high-volatility states (large ${\epsilon }_t^2$, high $d_t^{-}$), the k-NN selects historically similar spike-time nodes whose successors remain in the elevated-volatility tail, perpetuating the cluster.

4.2. Sensitivity to Neighbourhood Size k

Table 3 reports, for each $k\in \{10,20,30,50\}$, the median and 90% prediction interval $[p_5,p_{95}]$ of $\widehat{\sigma }$, $\widehat{\gamma }$, $\widehat{\kappa }$, and ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ across $N=500$ simulated paths. All other settings are unchanged.

For PUN, all three statistics and ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ are covered across the entire range $k\in \{10,20,30,50\}$: the simulation band consistently contains the observed values regardless of neighbourhood size. The median estimates remain close to the observed values, while the interval width narrows modestly as k increases, reflecting the smoothing effect of larger neighbourhoods.

For PJM, coverage holds for $k\le 30$ but deteriorates at $k=50$: the observed standard deviation ($0.180$), skewness ($1.178$), and kurtosis ($9.90$) fall above the respective 95th percentiles ($0.175$, $0.910$, $9.26$), while ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ remains covered. The failure is systematic: at $k=50$ the median skewness drops to $0.419$ and the median kurtosis to $5.75$, far below the observed values, indicating that the large neighbourhood averages over regimes with substantially different skewness and scale. This compresses the local SAS parameters toward a global mean and erases the heavy-tailed, right-skewed character of the PJM distribution. The first-order ACF of squared innovations is more robust to this averaging because volatility clustering is a coarser property than the precise shape of the marginal distribution.

These results confirm that $k=10$ is a conservative and well-motivated choice, and that the framework remains valid across $k\le 30$ for both markets. More generally, k is the natural adaptation parameter of the framework: markets with more heterogeneous dynamics may benefit from a smaller k to preserve local specificity, while larger k provides stability at the cost of smoothing over distinct regimes.

4.3. Volatility Clustering: Full ACF Profile

Figure 5 displays the median ACF(${\epsilon }_t^2$) profile over lags 1–30, together with the $[5\mathrm{th},95\mathrm{th}]$ percentile simulation band and the empirical profile.

For PUN (MAE-ACF $=0.080$), the empirical ACF oscillates between positive and negative values, consistent with a weak clustering signal (${\widehat{\rho }}_1=0.168$). The observed profile exits the $[p_5,p_{95}]$ band at a few lags; with 30 lags tested simultaneously at the 90% level, approximately three exceedances are expected even under a correctly specified model. The higher MAE-ACF relative to PJM reflects the oscillatory character of the observed profile: the median of 500 simulated paths is inherently smooth and cannot track the individual oscillations of a single realisation with a weak dependence signal. The coverage of ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ in Table 2 confirms that the overall clustering level is correctly reproduced.

For PJM (MAE-ACF $=0.031$), the empirical ACF exhibits a clear positive structure at short lags (${\widehat{\rho }}_1=0.210$), consistent with pronounced spike clustering. The simulation band covers the empirical profile across all 30 lags, and the median simulated ACF tracks the observed decay profile closely.

4.4. Simulated Price Paths

Figure 6 provides a qualitative assessment of the simulation framework by displaying $N=30$ simulated detrended log-price paths together with the observed series. The single path whose maximum and minimum most closely match the observed extremes is highlighted in blue; all other paths are shown as thin grey lines.

The fan of grey paths spans a plausible range of scenarios anchored to the empirical volatility structure, without any path collapsing to zero or diverging to unrealistic levels. The highlighted path (blue) reproduces the overall amplitude of the observed series, confirming that the local SAS extrapolation generates tail realisations consistent with the observed extremes. The clustering of large deviations visible in the observed series (particularly the 2022 energy crisis spike in PUN and the cold-wave spikes in PJM) reappears across the simulated fan.

5. Discussion

5.1. A Single Mechanism for All Distributional Properties

The central conceptual contribution of this paper is the identification of a single conditioning mechanism — topological k-nearest neighbours in the NVG backward feature space — that jointly determines the scale, skewness, and tail weight of the local innovation distribution. This is not merely an empirical observation but a structural consequence of the framework design.

All three SAS parameters $({\widehat{\lambda }}_i,{\widehat{\eta }}_i,{\widehat{\delta }}_i)$ are estimated jointly on the k neighbours of node i: they share the same conditioning variable, the same neighbourhood, and the same estimation procedure. No separate equations govern their dynamics. The mechanism is parsimonious at the global level — one conditioning variable, one neighbourhood size k — while being flexible at the local level: each node carries its own distributional triplet, adapted to the local volatility regime.

5.2. Locally Parametric and Distribution-Free

A natural question is whether the introduction of the SAS family compromises the locally parametric character of the framework. The answer depends on the level at which parametric assumptions are evaluated.

At the global level, no parametric assumptions are imposed: the total number of SAS parameters is $3n$ (one triplet per node), so it grows with the sample size — a defining property of non-parametric methods. No fixed distributional form is assumed for the innovation distribution across all observations.

At the local level, within each topological neighbourhood, a three-parameter SAS distribution is fitted to k successor values. This is a form of local maximum-likelihood estimation [27]: the parametric model is a working approximation to the local distribution, not a global structural assumption.

The SAS was selected for its parsimony and flexibility: three parameters suffice to represent a wide range of asymmetric, heavy-tailed distributions encountered in power markets. However, the conditioning mechanism is independent of this choice: any sufficiently rich family (Gaussian mixture, normal-inverse-Gaussian, generalised hyperbolic) would serve the same role, producing a different working approximation of the same locally varying distribution.

5.3. Structural Comparison with GARCH-Type Models

In GARCH-type models, the conditional variance is driven by a parametric recursion on past squared innovations and past variances, so scale adapts over time. Skewness and tail weight, however, remain global constants: a symmetric innovation distribution (e.g. Student-t) forces zero skewness regardless of the market state, and the tail parameter is fixed across all regimes. Any attempt to improve one component — adding skewness via a skewed-t innovation distribution [28], or richer scale dynamics via an EGARCH specification — necessarily re-parameterises the others. Further extensions with time-varying higher moments [29] add dedicated parametric equations for skewness and kurtosis dynamics, substantially increasing model complexity without providing a unified interpretable mechanism.

In the topological framework, scale ${\widehat{\lambda }}_i$, skewness ${\widehat{\eta }}_i$, and tail weight ${\widehat{\delta }}_i$ are all functions of the local topological state, estimated jointly on the k neighbours of node i. No global parameters govern the distributional shape; all three properties vary endogenously with the market regime through the geometry of the NVG, without any additional equations or parameters in the global model.

5.4. Market-Specific Interpretation

The two markets present complementary profiles that illuminate the behaviour of the proposed approach.

PUN (Italian power exchange) has moderate kurtosis ($\widehat{\kappa }=5.48$) and negative skewness ($\widehat{\gamma }=-0.351$), driven by the 2022 energy crisis which produced large but relatively symmetric price movements. The volatility clustering signal is moderate (${\widehat{\rho }}_1\left({\epsilon }_t^2\right)=0.168$). The topological conditioning captures the alternation between high- and low-volatility regimes; coverage of all three statistics and ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ is achieved with a wide $[p_5,p_{95}]$ band for skewness and kurtosis, reflecting the genuine sampling variability of these statistics at $T=1825$.

PJM West Hub exhibits stronger clustering (${\widehat{\rho }}_1\left({\epsilon }_t^2\right)=0.210$), higher kurtosis ($\widehat{\kappa }=9.90$), and pronounced positive skewness ($\widehat{\gamma }=+1.178$) from cold- and heat-wave spikes. The larger empirical skewness requires the local SAS fit to assign systematically positive ${\widehat{\eta }}_i$ in high-volatility regimes. The fact that coverage is achieved for both skewness and kurtosis simultaneously confirms that the topological neighbourhood correctly identifies the spike-clustering regimes in which both properties are most extreme.

6. Conclusions

This paper has introduced a framework for Monte Carlo simulation of power price innovations whose core principle is topological conditioning: at each simulation step, the local distribution of the next innovation is determined by the k-nearest neighbours of the current node in the Natural Visibility Graph feature space, rather than by a parametric recursion for the conditional variance or a global distributional assumption.

The main finding is that a single conditioning mechanism — k-nearest neighbours in the space $({\epsilon }_t^2,d_t^{-},{\ell }_t^{-})$ — simultaneously determines the scale, skewness, and tail weight of the local innovation distribution, through the Jones–Pewsey Sinh-ArcSinh family fitted on the topological neighbourhood. These three distributional properties are not modelled separately: they emerge jointly from one geometric operation, with no additional parameters or regime-switching equations. The framework is locally parametric — the total parameter count grows with the sample size — and the SAS family serves as a working model, not a structural assumption; any universal approximator could substitute it without altering the conditioning mechanism.

The empirical results confirm this principle. Applied to two power markets with contrasting distributional characteristics (PUN: negative skewness, moderate kurtosis; PJM: strong positive skewness, heavy tails), the framework achieves simultaneous coverage of $\widehat{\sigma }$, $\widehat{\gamma }$, $\widehat{\kappa }$, and ${\widehat{\rho }}_1\left({\epsilon }_t^2\right)$ — for both markets with a single neighbourhood size $k=10$ and no market-specific re-calibration. The observed values fall within the $[5\mathrm{th},95\mathrm{th}]$ percentile band of $N=500$ simulated paths of length T for all eight moment-market combinations.

The framework contrasts structurally with GARCH-type models, not merely in accuracy. In GARCH-skewed-t, scale varies through a parametric recursion while skewness and tail weight remain global constants. Extending GARCH to time-varying higher moments requires dedicated equations for each property, increasing parametric complexity without providing a unified mechanism. The topological approach subsumes all three into a single operation, with complexity emerging from the geometry of the NVG rather than from parametric specification.

The framework is portable by design: the same algorithm, with the same hyperparameters, was applied without modification to markets with opposite skewness signatures and different tail weights. The neighbourhood size k is the sole tunable parameter and serves as a natural adaptation lever for markets with substantially different distributional characteristics. The synthetic paths produced by the framework are directly applicable to operational tasks such as derivative pricing, Value-at-Risk estimation, and scenario generation for investment analysis under uncertainty.

Several limitations of the present study point to directions for future research. The analysis is restricted to two electricity markets; gas and oil markets, where kurtosis is dominated by rare extreme events [6], represent a natural next step. For such markets, a smaller k would better preserve the local specificity of spike regimes, but the SAS fit on very few points becomes unreliable: hybrid approaches combining topological conditioning with explicit extreme-value modelling for the far tails are a promising extension. The feature vector ${\mathbf{f}}_t=({\epsilon }_t^2,d_t^{-},{\ell }_t^{-})$ was chosen for parsimony; enriching it with additional backward topological features — mean backward edge length, lagged squared innovations, or rolling variance estimates — may sharpen the conditioning signal for markets with complex multi-step clustering patterns, though the trade-off between feature richness and overfitting in the k-NN search warrants careful evaluation. Data-driven selection of k via cross-validation or information-theoretic criteria would eliminate the sole remaining hyperparameter. Finally, the $O\left(T^2\right)$ NVG construction cost, feasible for the daily series considered here ($T\le 1825$), would require approximate nearest-neighbour structures for application to high-frequency intraday data.