Edge-Optimized Reservoir Computing for Forecasting Complex Ecological Systems: Predicting Fishery Collapse with Incomplete Data

1. Introduction

Modern socio-ecological system (SES) models encounter significant challenges due to high-dimensional interactions and non-linear dynamics that can hinder sustainable decision-making [1]. These complex dynamics generate feedbacks and critical regime thresholds [2], where variability in noisy time series can obscure biological reorganizations [3,4], making long-term monitoring and global sensitivity analysis essential to mitigate analytical bias [1,5]. A representative case is the tropical shrimp fishery, where the transition from artisanal methods to large-scale industrial trawling altered benthic habitats and food webs [6,7,8]. Overexploitation and environmental degradation have subsequently driven severe landing declines in tropical regions [9,10], often characterized by landing drops below $10\%$ of historical maximums and accompanying socio-ecological instability [11].

Recent studies have documented severe landing declines across tropical and subtropical shrimp fisheries driven by cumulative fishing pressure and management challenges, as observed in Brazil [12], Tanzania [13], and Bangladesh [14]. In the Colombian Caribbean, industrial shrimp trawling initiated during the 1960s and expanded during the 1970s and 1980s, supported by favorable market margins and international demand [15]. During its historical peak, this fleet constituted one of the primary industrial marine sectors in the region, where the deployment of Florida-type trawlers equipped with double-rigged nets increased sweep efficiency and operational catchability [15,16].

Limited data availability has historically hindered precise stock assessments and the implementation of effective management measures, preventing the application of robust evaluation models [17]. Consequently, the cumulative impacts of overexploitation were not detected in time, becoming evident in the late 1990s when reported landings began a sustained decline. Operational constraints, including exchange rate fluctuations, escalating fuel costs, and aquaculture market competition, led to the official closure of the fishery in 2005. The daily operational dynamics of commercial fleets yield highly noisy datasets with an elevated percentage of missing values [18,19], structured as complex, high-dimensional time series [20,21,22,23] that present a critical challenge for conventional predictive frameworks [24,25].

Traditional listwise deletion or linear interpolation methods are inadequate and introduce bias when applied to non-linear dynamic time series, as they ignore underlying missing data distributions and distort system trajectories by omitting temporal correlations [26]. In non-linear systems, the loss of temporal continuity compromises structural representation by failing to capture higher-order dependencies [19,21]. To address these limitations, the Reservoir Computing (RC) framework integrates an input layer, a non-linear dynamic reservoir with fixed randomized weights, and a linear readout layer [18,20]. Its primary computational advantage lies in training exclusively the output weights through linear optimization, such as Ridge regression or regularized least squares [21,25,27,28]. By avoiding backpropagation through time (BPTT) [29], the RC paradigm ensures rapid convergence without vanishing gradient phenomena and with low computational overhead, making it suitable for deployment on edge devices [26,30]. Furthermore, its capacity to manage data loss relies on iterative schemes grounded in fixed-point convergence, enabling the recovery of structural dynamics without exhaustive preprocessing [18,19,20,24].

Echo State Networks (ESN) represent a prominent Reservoir Computing framework for processing temporal data [23,30] by employing fixed, randomized weight matrices ($W_{\mathrm{in}},W$) and training exclusively the output readout layer $W_{\mathrm{out}}$ via linear optimization [18,20]. Dynamical stability relies on satisfying the Echo State Property (ESP) with spectral radius $\rho <1$ [20,26,30], enabling phase-space reconstruction and imputation in highly variable time series without BPTT, even under elevated missing data rates ($\ge 85\%$) or near non-linear transition thresholds [18,19,24,26,30,31]. To overcome Python execution constraints [32,33], functional differentiable programming via JAX [34] integrates automatic differentiation (grad) [35], explicit array vectorization (vmap) [36], and Just-In-Time (JIT) compilation powered by OpenXLA [37,38]. This implementation optimizes memory throughput and computational overhead through kernel operator fusion [39], achieving execution efficiency comparable to low-level C++/CUDA implementations [40,41]. Deploying a JAX-native ESN on an NVIDIA Jetson Orin Nano edge device enables real-time inference and dynamic variable imputation [33,39,41,42], supporting the historical landing series reconstruction of the Colombian Caribbean industrial shrimp fishery under a missing data rate of $37.10\%$ [43].

Despite these methodological advances, two gaps remain in the application of Reservoir Computing to data-scarce fisheries monitoring. First, comparative benchmarks of ESNs against conventional forecasting baselines (SARIMA, LightGBM, LSTM) have rarely been conducted under a strict edge-deployment constraint that jointly considers predictive accuracy, inference latency, memory footprint, and energy consumption. Second, while imputation-capable reservoir architectures have been proposed for missing data recovery [18,26], their use for long-range, bidirectional historical reconstruction of collapsed or discontinued fisheries—where monitoring gaps exceed one-third of the observed record—remains largely unexplored, particularly in tropical small-scale fisheries with fragmented statistical systems [43].

To address these gaps, this study proposes an edge-optimized ESN framework implemented in JAX and evaluates it through two distinct experimental tasks applied to the historical landing series of the Colombian Caribbean industrial shrimp fishery. The first task addresses the question of which model offers the best one-step-ahead forecasting performance under a strict out-of-sample Walk-Forward Validation (WFV) scheme, benchmarking the JAX-native ESN against SARIMA, LightGBM, and LSTM baselines in terms of predictive accuracy (MAE, RMSE, RMSSE) and edge-deployment cost (latency, memory, energy) on an NVIDIA Jetson Orin Nano. The second task addresses the distinct question of whether a bidirectional ESN architecture can retrospectively reconstruct large, structurally discontinuous historical gaps (missing data rate of $37.10\%$) in the same series, providing a continuous landing trajectory suitable for exploratory comparison against independent biological and stock-assessment evidence, without being used as a substitute for formal stock assessment.

The main contributions of this work are threefold: (i) a systematic, hardware-aware benchmark of Reservoir Computing against conventional statistical and deep learning baselines for fisheries time-series forecasting under small-data regimes; (ii) a bidirectional ESN reconstruction protocol for long-gap univariate marine landing series, evaluated independently from the forecasting task; and (iii) an open, energy- and latency-characterized deployment of a JAX-native ESN on embedded edge hardware, establishing a reproducible computational baseline for autonomous, low-resource fisheries monitoring systems.

2. Materials and Methods

2.1. Historical Landing Records (1993–2023)

Landings of Shallow-Water Shrimp (SWS) taken by the industrial trawl fleet along the Colombian Caribbean coast comprise several commercial taxa, including Penaeus notialis, Penaeus schmitti, Penaeus brasiliensis, Penaeus subtilis, Trachypenaeus similis, Trachypenaeus constrictus, Xiphopenaeus kroyeri, and Penaeus monodon [43,44]. A historical database containing monthly landing records from 1993 to 2023 was analyzed [43]. Data processing involved descriptive analyses at annual and monthly scales, including inter- and intra-annual comparisons of total landings and taxonomic composition.

The historical trajectory of reported landings shows a peak during the late 1990s, reaching a maximum of $179.82\mathrm{t}$ in 1997, followed by a decline that led to an operational collapse between 2000 and 2005 (Figure 1). Figure 1 also shows extensive observational gaps resulting from administrative and operational disruptions. The empirical landing distribution exhibits a strong positive skewness (Figure 2), demonstrating that the fishery operated predominantly at low yield levels (between 0 and $30\mathrm{t}$), whereas high landing values ($>100\mathrm{t}$) corresponded to isolated events prior to the decline in catches.

Intra-annual variability evaluated across months (January to December) reflects a seasonal pattern during historical high-landing periods, characterized by catch peaks in February, April, and December (Figure 3). Over the last decade, this seasonal signal flattened near zero, indicating a suppression of intra-annual landing fluctuations.

The taxonomic breakdown highlights shifts in taxon contributions and temporal dominance (Figure A1), notably the prominence of Pleoticus robustus during the late 1990s. Disaggregating these trajectories into individual panels (Figure 4) reveals taxon-specific landing dynamics and reporting discontinuities, demonstrating that the decline did not affect all reported categories uniformly or simultaneously.

Comparative boxplot distributions (Figure 5) confirm that Pleoticus robustus exhibited the highest median landings and interquartile ranges across the time series, whereas taxa such as Penaeus schmitti and Penaeus notialis maintained consistently lower and narrower distributions.

Table 1 summarizes the descriptive statistics for global landings and individual target taxa. All analyzed variables exhibit strong positive skewness, with mean values exceeding their respective medians. At the taxonomic level, Pleoticus robustus recorded the highest median landings ($44.81\mathrm{t}$) and a historical maximum of $177.93\mathrm{t}$. Conversely, Penaeus species (P. notialis and P. schmitti) yielded lower values, with medians of $4.00\mathrm{t}$ and $2.86\mathrm{t}$, respectively. The dataset shows missing observations across taxa, ranging from $28.47\%$ in P. notialis to $61.67\%$ in the general Penaeidae category.

2.1.1. Reconstructed Time Series Baseline

To enable rigorous benchmarking against classical statistical and machine learning models that strictly require complete time series without missing values, a fully reconstructed historical landing dataset was generated for the 1993–2023 period. Prior to imputation, a domain-specific temporal feature matrix was engineered to allow the non-parametric missForest algorithm [45] to capture complex temporal dynamics without imposing parametric assumptions. The target landing variable was contextualized by incorporating a continuous time index ($t\in [1,N]$) to model secular long-term trends; an interannual component (year); a categorical factor variable representing calendar months ($\mathrm{month}\in \{1,\dots ,12\}$); and harmonic sinusoidal transformations ($sin(2\pi m/12)$ and $cos(2\pi m/12)$, where m denotes the calendar month) to eliminate artificial boundary discontinuities between December and January.

Data imputation in R (version 4.2.2) [46] using the missForest package (version 1.6.1) was executed with 300 trees ($\mathrm{ntree}=300$), 3 variables randomly sampled at each split ($\mathrm{mtry}=3$), and a maximum limit of 15 iterations ($\mathrm{maxiter}=15$). Driven by this explicit temporal and seasonal representation, the imputation scheme yielded low internal imputation error, achieving a Normalized Root Mean Squared Error (NRMSE) of $0.0113$ for continuous variables and a Proportion of Falsely Classified (PFC) entries of $0.0000$ for categorical variables.

As shown in Figure 6, the imputed series restores temporal continuity across historical administrative gaps while preserving non-linear trajectories, peak events, and operational collapse dynamics. The statistical summary of the reconstructed dataset is presented in Table 2, yielding a complete set of $N=372$ monthly observations ($0\%$ missing values) with a mean of $40.12\mathrm{t}$ and a median of $32.72\mathrm{t}$, providing a fully gap-filled baseline for traditional benchmark models.

2.2. Experimental Scenarios and Dataset Partitioning

To evaluate model performance under realistic operational conditions and standardized benchmark configurations, an experimental scenario was designed using the reconstructed historical landing time series (${\mathrm{Landing}}_{\mathrm{rec}}$). To ensure a temporally consistent forecasting evaluation, the monthly historical landing time series ($N=372$ observations spanning January 1993 to December 2023) was chronologically partitioned into training, validation, and test subsets (Table 3). The training split encompasses 25 years (1993–2017, $n=300$ months, $80.65\%$), serving as the historical learning window. The validation set covers 3 years (2018–2020, $n=36$ months, $9.68\%$) for hyperparameter tuning and early stopping. Finally, an independent out-of-sample test set spanning 3 years (2021–2023, $n=36$ months, $9.68\%$) was reserved strictly for final performance evaluation.

To quantify the influence of the reconstruction stage, the proportion of reconstructed observations was calculated for each data partition. Most reconstructed values were concentrated within the historical training period (127 of 300 observations; 42.33%), corresponding to long monitoring interruptions prior to 2018. In contrast, the independent test partition contained only 2 reconstructed observations (5.56%), indicating that 94.44% of the final out-of-sample evaluation was based on directly observed landing records. Therefore, the reported forecasting performance is predominantly supported by empirical observations rather than imputed values. The ${\mathrm{Landing}}_{\mathrm{rec}}$ scenario utilizes the fully reconstructed historical landing dataset generated via non-parametric imputation with missForest (Figure 6), enabling a fair and consistent comparison among forecasting architectures under identical data conditions.

As shown in Table 4, this scenario reflects a regime shift in recent years. While historical training records (1993–2017) exhibit higher catch levels (mean of $48.81\mathrm{t}$), the validation (2018–2020) and test (2021–2023) periods show lower landing values ($2.08\mathrm{t}$ and $4.54\mathrm{t}$, respectively). The non-parametric imputation procedure preserves the distributional properties, quantiles, and extreme bounds of each partition without inducing variance collapse or tail distortion.

To evaluate the stationarity properties and order of integration of the target landing time series prior to model training, Augmented Dickey-Fuller (ADF) unit root tests were performed using the tseries package (version 0.10-62) [47] on the raw ${\mathrm{Landing}}_{\mathrm{raw}}$ series (excluding missing values) and on the reconstructed ${\mathrm{Landing}}_{\mathrm{rec}}$ series across both the full 1993–2023 observation window ($N=372$) and the isolated training partition ($n=300$). Across the entire evaluation period, the raw observed series rejected the unit root null hypothesis (${\mathrm{Landing}}_{\mathrm{raw}}$, $DF=-3.9973$, $\mathrm{lag}=6$, $p<0.01$). Similarly, the reconstructed series exhibited mean-reverting behavior (${\mathrm{Landing}}_{\mathrm{rec}}$, $DF=-5.9306$, $\mathrm{lag}=7$, $p<0.01$). These test statistics indicate that the imputation procedure preserved the overall mean-reverting properties without introducing deterministic linear trend artifacts into the univariate series.

2.2.1. Walk-Forward Validation

To prevent data leakage and simulate a realistic operational deployment, model parameters and hyperparameters were evaluated using Walk-Forward Validation (WFV) with a rolling-origin evaluation scheme [48]. Unlike standard k-fold cross-validation—which violates temporal dependencies by randomizing observations—the WFV protocol preserves chronological directionality [49]. Under this framework, parameter estimation at any target time step t relies exclusively on observable historical data.

The sequence evaluation procedure is formalized by establishing an initial training window of length $T_{\mathrm{init}}$, defining the operational observations required for model initialization and stable feature scaling. At each iteration i, the model utilizes the historical landings subsequence $y_{1:t}$ to generate multi-step out-of-sample forecasts over an h-step horizon, defined as ${\widehat{y}}_{t+1:t+h}=f\left(y_{1:t}\right)$. Following forecast generation over h, the forecast origin advances by s discrete time steps ($t\leftarrow t+s$). An expanding window strategy is implemented, wherein subsequent training iterations cumulatively integrate historical observations from preceding steps, enabling the model to adapt to structural variations and non-stationary trends in reported landings as new data become available.

The paired sequence of predicted $\widehat{y}$ and observed landings y generated at each forecast origin is stored sequentially to construct the out-of-sample error matrix. Predictive accuracy is quantified by aggregating forecasting discrepancies across all evaluated temporal folds, providing an empirical assessment of model generalization capacity under operational deployment conditions.

2.3. Bidirectional Landing Reconstruction Experiment

To evaluate the capability of the ESN to capture the underlying non-linear dynamics and temporal dependencies of reported landings from fragmented observations, a bidirectional reconstruction experiment was designed using the raw landing time series. The historical record was partitioned into three continuous subseries corresponding to effective monitoring periods: $s_1$ (1993–2000), $s_2$ (2006–2009), and $s_3$ (2013–2023). Table 5 presents the statistical summary for each subseries. The reconstruction of data gaps between adjacent subseries was performed using a dual-estimation strategy combining an autonomous forward forecast (forecasting) from the antecedent subseries ($s_A$) and a reversed backward forecast (backcasting) from the subsequent subseries ($s_B$). This approach leverages the high-dimensional dynamics of the reservoir in both temporal directions, mitigating error accumulation inherent to long-term recursive forecasting over intervals with high missing data density. This procedure was applied exclusively to the ESN architecture to evaluate its intrinsic capacity for time-series reconstruction in data-scarce settings.

To ensure smooth continuity across the boundaries of observed data and eliminate structural discontinuities at the center of the unmonitored gaps, the bidirectional forecasts were integrated using a soft temporal weighting scheme. For each time step $t\in [1,T]$ within a missing interval of length T, the consolidated estimate $\widehat{y}\left(t\right)$ is defined as:

$$\widehat{y}\left(t\right)=w\left(t\right)·\widehat{y}*\mathrm{forward}\left(t\right)+\left(1-w\left(t\right)\right)·\widehat{y}*\mathrm{backward}\left(t\right)$$

(1)

where the linear decay weight $w\left(t\right)$ decreases monotonically from 1 to 0 according to:

$$w\left(t\right)=1-{\displaystyle \frac{t-1}{T-1}}$$

(2)

Through this formulation, the prediction assigns higher confidence to the forward forecast near block $s_A$ and transitions smoothly toward reliance on the backward forecast near the boundary of $s_B$, preserving physical continuity and local temporal dynamics at the edges of the observed data.

2.4. NVIDIA Jetson Orin Nano

To evaluate edge deployment feasibility and operational efficiency, inference and temporal validation were executed on an NVIDIA Jetson Orin Nano embedded module (8 GB variant). This platform features a 6-core Arm Cortex-A78AE v8.2 CPU running at a nominal frequency of 1.5 GHz (up to 1.73 GHz) and an NVIDIA Ampere architecture GPU with 1024 CUDA cores and 32 Tensor cores. Hardware acceleration for vectorized matrix operations and non-linear state updates is provided within this system, which is required for the real-time execution of recurrent reservoir computing architectures.

The platform integrates 8 GB of unified LPDDR5 memory shared across the CPU and GPU cores, alongside an NVMe solid-state drive (SSD) to deliver the throughput required for continuous data ingestion without volatile memory bottlenecks. Memory consumption metrics reported during execution reflect resource allocations within this shared physical RAM pool. Operating within a configurable power budget of 7–15 W (15 W mode enabled), the system provides an energy-efficient deployment platform suitable for integration into autonomous environmental monitoring stations and field-level edge acquisition systems.

The software execution environment was standardized on NVIDIA JetPack 6.2.3 (L4T R36.5.2) running Linux Kernel 6.8 (aarch64), accelerated via CUDA Toolkit 12.6. Model inference and reservoir dynamics were natively compiled and executed in Python using JAX (v0.5.2) configured with full CUDA GPU backend support (CudaDevice). The technical characterization of the hardware and software stack is detailed in Table 6.

2.5. Model Architecture and Edge Optimization

The ESN architecture comprises three distinct processing stages: a fixed input mapping, a dynamic non-linear reservoir, and a linear output readout layer. Mathematically, the continuous state space and output transformations are defined by an input weight matrix $W_{\mathrm{in}}\in {\mathbb{R}}^{N_x\times N_u}$, a recurrent internal adjacency matrix $W\in {\mathbb{R}}^{N_x\times N_x}$, and an output readout weight matrix $W_{\mathrm{out}}\in {\mathbb{R}}^{N_y\times (N_x+N_u+1)}$.

The discrete-time dynamics of the reservoir internal state $x(n+1)\in {\mathbb{R}}^{N_x}$ are updated according to the non-linear state transition equation:

$$x(n+1)=(1-\alpha )x\left(n\right)+\alpha tanh\left({\gamma }_{\mathrm{in}}{W}_{\mathrm{in}}u(n+1)+Wx\left(n\right)+b\right)$$

(3)

where:

• $u(n+1)\in {\mathbb{R}}^{N_u}$ represents the vector of reported landing features at time step $n+1$;
• $\alpha \in (0,1]$ is the leaking rate governing system memory dissipation across temporal scales;
• ${\gamma }_{\mathrm{in}}$ denotes the input scaling factor controlling the degree of non-linearity activated within the reservoir state space;
• $b\in {\mathbb{R}}^{N_x}$ is the bias vector, and $tanh(·)$ serves as the element-wise non-linear activation function mapping inputs into a high-dimensional state space (Equation (4)).

$$tanh\left(x\right)={\displaystyle \frac{e^x-e^{-x}}{e^x+e^{-x}}}.$$

(4)

To ensure system stability and maintain the fading memory property, the reservoir state space must satisfy the ESP. This guarantees that the internal state representation evolves as a unique function of the historical input sequence while asymptotically dissipating initial state influences. Mathematically, a sufficient condition for ESP compliance is enforced by scaling the internal recurrent adjacency matrix W such that its spectral radius $\rho \left(W\right)$ satisfies $\rho \left(W\right)<1$.

The training phase is restricted to optimizing the readout weight matrix $W_{\mathrm{out}}\in {\mathbb{R}}^{N_y\times (N_x+N_u+1)}$. This step maps the concatenated reservoir state matrix $X\in {\mathbb{R}}^{(N_x+N_u+1)\times T}$ to the target reported landing matrix $Y\in {\mathbb{R}}^{N_y\times T}$ using Ridge regression ($L_2$ Tikhonov regularization) to prevent overfitting and ensure numerical stability during matrix inversion:

$$W_{\mathrm{out}}=Y{X}^T{(X{X}^T+\beta I)}^{-1},$$

(5)

where $\beta >0$ denotes the Tikhonov regularization parameter and I represents the identity matrix.

2.5.1. Hyperparameter Tuning and Experimental Setup

Input sequences were constructed using a lag depth of $N_u=5$ historical observations to project the univariate landing signal into a high-dimensional recurrent reservoir space ($N_x=500$). Hyperparameter selection and reservoir initialization were implemented following the guidelines of Lukoševičius (2012) [50].

To configure internal network dynamics and prevent temporal data leakage, a grid search was conducted using a forward-chaining cross-validation split. The search space included settings for reservoir capacity $N_x\in \{100,250,500,1000\}$, leaking rate $\alpha \in \{0.1,0.3,0.5,0.8,1.0\}$, spectral radius $\rho \in \{0.85,0.90,0.95,0.99\}$, input scaling ${\gamma }_{\mathrm{in}}\in \{0.01,0.1,0.5,1.0\}$, and Ridge regularization coefficient $\beta \in \{{10}^{-4},{10}^{-3},{10}^{-2},{10}^{-1},1.0\}$. The optimal combination was selected based on the lowest Mean Squared Error (MSE) of validation predictions, evaluated on non-imputed historical landing data.

To remove transient initialization artifacts, a washout period of $T_{\mathrm{wash}}=200$ steps was discarded before fitting the readout model. The readout parameters were estimated via closed-form Ridge regression (Tikhonov regularization) using direct matrix inversion, which provides low computational overhead suitable for edge devices. The complete structural hyperparameters, dynamic settings, and operational parameters configured for the JAX-accelerated ESN architecture are detailed in Table 7.

2.6. Baseline Models

To evaluate the predictive performance of the proposed ESN relative to established modeling approaches, three baseline architectures were implemented: a linear seasonal autoregressive baseline, a tree-based ensemble framework, and a deep recurrent neural network.

2.6.1. Seasonal Autoregressive Integrated Moving Average

To model the statistical properties of the time series, a Seasonal Autoregressive Integrated Moving Average (SARIMA) model was evaluated as a linear benchmark. Hyperparameter selection and order estimation for seasonal and non-seasonal autoregressive, differencing, and moving average components were performed using the pmdarima library (version 2.1.1) [51] through Akaike Information Criterion (AIC) minimization. Table 8 outlines the parameter estimates and diagnostic metrics for the fitted $\mathrm{SARIMA}(2,1,3)\times {(2,0,0)}_{12}$ model evaluated on the landing time series ($N=372$). Although most parameters achieved statistical significance ($p<0.01$) and residual autocorrelation was reduced (Ljung-Box $Q\left(1\right)=0.01$, $p=0.910$), the Jarque-Bera test statistic ($JB=1106.14$, $p<0.001$) indicates heavy-tailed residual behavior, reflecting the abrupt regime shifts present in historical landing data.

2.6.2. Light Gradient Boosting Machine

To evaluate non-linear tree-based ensemble methods, LightGBM was implemented using the lightgbm library (version 4.7.0) [52]. Because Gradient Boosted Decision Trees (GBDT) lack intrinsic sequential memory, temporal dependencies were explicitly incorporated by extracting autoregressive lags alongside calendar features. The predictor matrix includes historical target lags ($t-1$ through $t-5$) in conjunction with temporal covariates—specifically year, month, and sine/cosine cyclical encoding—to capture multi-scale seasonality and long-term landing trends.

Structural hyperparameters were tuned using randomized search with temporal forward-chaining cross-validation (TimeSeriesSplit) evaluating $N=50$ parameter combinations on the training set, optimizing for Root Mean Squared Error (RMSE). To mitigate overfitting, an early stopping criterion of 100 boosting rounds monitored against the validation set was enforced during hyperparameter tuning and final model training. The optimal hyperparameter configuration is summarized in Table 9. Regularization penalties ($L_1$ and $L_2$), feature fractioning, and subsampling were applied to constrain tree depth and prevent overfitting under data-scarce conditions.

2.6.3. Long Short-Term Memory

A Long Short-Term Memory (LSTM) recurrent neural network was evaluated as a deep learning benchmark for sequential modeling. To reformulate the time series into a supervised learning structure, a sliding-window approach was applied, constructing input sequences of $T=10$ consecutive historical timesteps ($X_t=[y_{t-T+1},\dots ,y_t]$) to predict the single-step-ahead landing value ($y_{t+1}$). The architecture was implemented using Keras (version 3.15.1) [53] with JAX as the computational backend. The JAX framework enabled accelerated network training and Backpropagation Through Time (BPTT) via JIT compilation and vectorized operations, ensuring computational efficiency during evaluation.

The network topology comprises an initial LSTM layer with ReLU activation, followed by a Dropout regularization layer and a dense linear output layer. Hyperparameter selection was conducted via random search across $N=10$ sampled configurations from a discrete search space spanning hidden unit capacities $u\in \{32,64,100,128\}$, learning rates $\eta \in \{{10}^{-2},5\times {10}^{-3},{10}^{-3},5\times {10}^{-4},{10}^{-4}\}$, batch sizes $B\in \{16,32,64\}$, training epochs $E\in \{20,50,100\}$, and dropout rates $p_{\mathrm{drop}}\in \{0.0,0.2,0.3,0.4,0.5\}$. Models were compiled using Mean Squared Error (MSE) loss and trained with the Adam optimizer. The optimal candidate was selected based on the minimal RMSE computed on denormalized validation predictions. Table 10 details the resulting optimal configuration, comprising $66,689$ trainable parameters.

2.7. Performance Metrics

To evaluate predictive accuracy across scale-dependent and scale-independent dimensions, both point-error metrics and relative scaled forecasting indicators were implemented. Initial forecasting accuracy was quantified using the Mean Absolute Error (MAE) and the Root Mean Squared Error (RMSE), defined in Equations (6) and (7), respectively:

$$\mathrm{MAE}={\displaystyle \frac{1}{h}}{\displaystyle \sum _{t=n+1}^{n+h}}|y_t-{\widehat{y}}_t|$$

(6)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{h}}\sum _{t=n+1}^{n+h}{(y_t-{\widehat{y}}_t)}^2}$$

(7)

While RMSE imposes a quadratic penalty that heavily weights larger forecast errors, MAE provides a linear assessment that is less sensitive to extreme transient peaks. However, because the historical landing series exhibits pronounced positive skewness, intermittency, and near-zero values (Figure 2)—a statistical distribution analogous to the intermittent demand series analyzed in the M5 forecasting competition [54]—standard scale-dependent metrics alone are insufficient to evaluate model performance across varying catch dynamics without bias.

To address these constraints, performance was also evaluated using the Root Mean Squared Scaled Error (RMSSE) [55], defined as:

$$\mathrm{RMSSE}=\sqrt{{\displaystyle \frac{\frac{1}{h}{\sum }_{t=n+1}^{n+h}{(y_t-\widehat{y}*t)}^2}{\frac{1}{n-1}\sum *{t=2}^n{(y_t-y_{t-1})}^2}}}$$

(8)

The RMSSE scales the out-of-sample mean squared error of the forecast against the in-sample mean squared error of a first-order Naïve random-walk model [56]. This metric provides two main methodological advantages for fisheries time series: it is scale-independent, allowing comparison across different scenarios, and it avoids division-by-zero singularities that destabilize percentage-based metrics during zero-catch or low-landing periods.

In Equations (7)–(8), $y_t$ represents the observed landing value at time step t, ${\widehat{y}}_t$ denotes the forecast value, h indicates the out-of-sample forecast horizon ($h=36$ months), and n specifies the duration of the historical training sequence ($n=300$ months).

3. Results

3.1. Static Benchmark Performance

Table 11 presents the comparative evaluation of predictive performance across candidate architectures throughout the training, validation, and testing phases. Regarding global predictive accuracy, the ESN architecture demonstrates superior performance, achieving the lowest test error metrics with a MAE of $1.90$, RMSE of $2.67$, and RMSSE of $0.67$. This indicates the capacity of the dynamic reservoir to capture complex non-linear patterns without overfitting. LightGBM exhibits competitive test performance (MAE of $2.34$, RMSE of $2.93$), while recording the lowest training error (MAE of $9.47$, RMSE of $19.33$). Conversely, SARIMA (test MAE of $6.18$) and LSTM (test MAE of $6.23$) display substantially higher errors, highlighting the limitations of linear statistical formulations and the convergence challenges of deep recurrent architectures on this dataset scale.

When evaluating computational efficiency during real-time execution, LightGBM demonstrates the lowest processing latency, recording a mean sample latency of $0.07$ ms/sample across validation and test sets. This processing rate is approximately four times faster than ESN and SARIMA ($0.27\mathrm{ms}/\mathrm{sample}$) and over five times faster than LSTM ($0.38\mathrm{ms}/\mathrm{sample}$). All tested architectures maintain sub-millisecond execution latencies during the testing phase, supporting their preliminary technical feasibility for deployment on resource-constrained embedded monitoring systems.

Regarding memory consumption, LightGBM maintained the lowest baseline system RAM occupancy ($42.24\%$) across all evaluated execution splits. The ESN architecture preserved a stable baseline system memory load of $46.85\%$, exhibiting only a minor dynamic memory peak during training ($\Delta {\mathrm{RAM}}_{\mathrm{peak}}=4.34\mathrm{MB}$) associated with state matrix allocation, which was completely released during inference ($\Delta {\mathrm{RAM}}_{\mathrm{peak}}=0\mathrm{MB}$). In contrast, SARIMA exhibited the highest baseline system RAM load ($50.01\%$ in testing), whereas LSTM required the largest peak allocation during training ($\Delta {\mathrm{RAM}}_{\mathrm{peak}}=6.01\mathrm{MB}$). These empirical findings support the use of ESN when maximum predictive precision is required, while LightGBM represents an effective alternative for deployment scenarios with stringent latency and computational constraints.

3.2. Trajectory Analysis and Physical Domain Consistency

Figure 7 illustrates the temporal trajectory of predicted landing series across all evaluated models, while Table 12 summarizes their descriptive statistical properties over their respective effective evaluation sample sizes. A fundamental operational requirement for physical consistency in fishery landing forecasting is the non-negativity of catch predictions ($\mathrm{Landing}\ge 0$). Linear baseline models failed to strictly satisfy this domain constraint; SARIMA yielded unphysical negative values down to $-12.60\mathrm{t}$, and LightGBM also crossed the zero threshold with a minimum of $-1.78\mathrm{t}$. Conversely, both recurrent architectures—LSTM ($\mathrm{Min}=8.80\mathrm{t}$) and ESN ($\mathrm{Min}=0.77\mathrm{t}$)—preserved strictly positive predictions, demonstrating inherent stability in low-density regimes.

Regarding distributional behavior, all candidate architectures matched the right-skewed nature of the target landing series ($\mathrm{Mean}>\mathrm{Median}$), maintaining central values consistent with historical landings ($\mathrm{Mean}\approx 38.15--40.95\mathrm{t}$). However, discrepancies arose in the reconstruction of extreme peak events (Max). While LightGBM exhibited severe oversmoothing by truncating peak values at $103.00\mathrm{t}$, ESN ($\mathrm{Max}=140.11\mathrm{t}$) and LSTM ($\mathrm{Max}=132.10\mathrm{t}$) tracked high-amplitude nonlinear spikes without numerical explosion or dynamic range attenuation. Finally, the effective evaluation sample sizes for ESN ($N=357$) and LSTM ($N=342$) reflect the initial truncation induced by sequence lookback requirements ($T_{\mathrm{lookback}}$) and the initial transient state discard period ($T_{\mathrm{washout}}$).

Statistical significance tests revealed a trade-off between in-sample fitting and out-of-sample generalization. During the training period (1993–2017), LightGBM exhibited lower fitting errors than the ESN model, as confirmed by positive Diebold-Mariano statistics ($DM=+2.2966,p=0.0223$ for squared error; $DM=+6.0595,p=4.173\times {10}^{-9}$ for absolute error). Across the out-of-sample validation and test period (2018–2023), however, the performance trend reversed, and the JAX-based ESN achieved superior predictive accuracy compared to LightGBM. The Diebold-Mariano test yielded negative statistics for both squared-error loss ($DM=-3.8674,p=2.698\times {10}^{-4}$) and absolute-error loss ($DM=-4.4134,p=4.221\times {10}^{-5}$), rejecting the null hypothesis of equal forecast accuracy. The paired Wilcoxon signed-rank test further confirmed that the ESN model produced a significantly smaller median absolute error ($V=406,p=3.217\times {10}^{-5}$). These results show that while LightGBM overfits when trained on tabular lag structures, the ESN architecture maintains dynamic generalization under time-varying conditions, supporting its suitability for operational deployment.

3.3. Walk-Forward Validation

To ensure a strictly equal baseline comparison, all models (SARIMA, LightGBM, LSTM, and ESN) were evaluated exclusively on the fully reconstructed dataset (${\mathrm{Landing}}_{\mathrm{rec}}$) generated via missForest.

Table 13 presents the evaluation metrics under the Walk-Forward Validation (WFV) scheme, which models sequential real-time deployment conditions through one-step-ahead ($h=1$) retraining and inference. The ESN yielded the lowest overall error metrics, attaining a MAE of $1.37$, RMSE of $1.96$, and RMSSE of $0.70$. This performance demonstrates the ability of the dynamic reservoir to update its readout weights via Ridge regression at each iterative step. SARIMA (MAE of $1.83$) and LSTM (MAE of $1.92$) showed lower errors than in static evaluation baselines, whereas LightGBM recorded a MAE of $2.01$ and RMSE of $2.47$.

For computational efficiency and potential edge-device retraining, LightGBM recorded the lowest overall execution time ($24.41\mathrm{s}$) and average latency per evaluation fold ($339.08\mathrm{ms}$). The ESN displayed comparable temporal efficiency, processing evaluation folds in $30.14\mathrm{s}$ ($424.44\mathrm{ms}/\mathrm{fold}$), which represents a minor computational trade-off given its forecast accuracy gains under the WFV protocol. Conversely, the LSTM required $1388.81\mathrm{s}$ ($19.29\mathrm{s}/\mathrm{fold}$) due to backpropagation through time and computation graph operations. SARIMA proved operationally unfeasible for continuous edge deployment, requiring $4348.07\mathrm{s}$ ($60.39\mathrm{s}/\mathrm{fold}$) due to iterative ARIMA parameter re-estimation at each step.

Regarding memory dynamics during WFV, LightGBM and LSTM maintained zero net memory leakage ($\Delta {\mathrm{RAM}}_{\mathrm{WFV}}=0.00\mathrm{MB}$), with LightGBM maintaining overall system usage at $70.01\%$. The ESN exhibited minor state dynamic variation ($\Delta {\mathrm{RAM}}_{\mathrm{WFV}}=111.94\mathrm{MB}$) due to reservoir state matrix updating, keeping overall system memory occupancy stable at $70.12\%$. Although SARIMA exhibited the lowest initial baseline memory footprint ($58.00\%$) due to the absence of deep learning framework dependencies, it showed significant state accumulation over continuous iterations ($\Delta {\mathrm{RAM}}_{\mathrm{WFV}}=631.91\mathrm{MB}$), indicating potential memory growth over extended deployment horizons.

Power profiling revealed distinct thermal and electrical behaviors across architectures. While average power consumption remained within a narrow band across all algorithms ($6.13\mathrm{W}$ to $6.60\mathrm{W}$), the LSTM recorded the highest transient power peak ($7.91\mathrm{W}$), driven by burst GPU compute cycles during recurrent backpropagation. Evaluating total energy expenditure ($E=\mathrm{Avg}.\mathrm{Power}\times \mathrm{Total}\mathrm{Time}$) highlighted the operational efficiency of reservoir computing: LightGBM and ESN consumed $161.11\mathrm{J}$ and $198.32\mathrm{J}$, respectively, across the WFV scheme. In contrast, the LSTM consumed $8513.41\mathrm{J}$ ($\approx 2.36\mathrm{Wh}$), and SARIMA required $27001.51\mathrm{J}$ ($\approx 7.50\mathrm{Wh}$). Consequently, the ESN provided an effective balance between predictive precision and energy efficiency for resource-constrained edge hardware.

3.4. Bidirectional Landing Reconstruction Results

The ESN-based bidirectional reconstruction strategy successfully imputed the unmonitored historical gaps (2000–2005 and 2010–2012) (Figure 1, Table 1). As illustrated in Figure 8, the reconstructed landing trajectory preserves the temporal variability and dynamics observed in adjacent subseries.

As summarized in Table 14, the reconstructed time series ($N=372$ monthly observations) exhibits a mean of $39.98\mathrm{t}$ ($\mathrm{Median}=38.13\mathrm{t}$), ranging from $0.10$ to $179.82\mathrm{t}$. For the 107 reconstructed observations, confidence intervals derived from the ensemble dispersion of multiple reservoir initializations indicate bounded uncertainty, with a mean lower boundary of $43.98\mathrm{t}$ for the $95\%$ confidence interval and a mean upper boundary of $61.39\mathrm{t}$.

3.4.1. Seasonal Patterns and Dynamic Transition

During the first reconstructed interval (2000–2005), the ESN architecture captured both the long-term memory and the seasonal fluctuations characteristic of the fishery landing records. Following the cessation of monitoring in February 2000, when a relative peak of $136.65\mathrm{t}$ was recorded, the reconstructed trajectory exhibited a progressive attenuation while preserving intra-annual oscillations. Systematic increases in reported landing levels occurred during the final quarter of each year and the initial months of the following year, reaching secondary peaks of $66.12\mathrm{t}$ in July 2002 alongside recurring recoveries in early 2003 and 2004 at $59.70\mathrm{t}$ and $45.53\mathrm{t}$, respectively. These upward trends were followed by moderate mid-year troughs, culminating in a decline to $26.20\mathrm{t}$ in November 2005 prior to the junction with subseries $s_2$.

This behavior indicates that the high-dimensional state space of the reservoir performs geometric interpolation between boundaries while projecting the dominant frequencies learned during training on subseries $s_1$. Soft Temporal Weighting ($w\left(t\right)$) eliminated sharp discontinuities at the center of the imputation blocks, ensuring a smooth transition into $s_2$ in 2006.

3.4.2. Uncertainty Behavior and Resource Status

The morphology of the confidence intervals (Figure 8) displays an hourglass pattern consistent with dynamic systems theory: uncertainty remains minimal near the boundaries adjacent to observed subseries ($s_A$ and $s_B$) and reaches its maximum amplitude at the center of the gap (around 2002–2003 and 2011), where the temporal distance to real observations is greatest.

Regarding fishery dynamics, the continuous reconstruction confirms that throughout the 2000–2012 period, relative landing levels remained above the operational landing-based collapse threshold of $17.98\mathrm{t}$ ($10\%$ of historical maximum landing). The lower boundary of the $95\%$ confidence interval ($\mathrm{Min}=14.32\mathrm{t}$) touched this critical threshold only during brief transition states in late 2005. This indicates that the subsequent decline in reported landings occurred rapidly after 2013 rather than as a continuous degradative process during the unmonitored years.

3.5. Model Stability and Seed Sensitivity Analysis

To evaluate whether the reconstruction performance is dependent upon favorable random state initializations, a comprehensive seed sensitivity analysis was conducted across ten independent random instances ($N_{\mathrm{seeds}}=10$). For each seed $s\in \{0,\dots ,9\}$, the internal weight matrices (W and $W_{\mathrm{in}}$) were stochastically generated utilizing JAX pseudo-random keys, while maintaining all structural hyperparameters strictly constant ($\rho =0.99,\alpha =0.3,N_x=500$).

As detailed in Table 15, the WFV metrics exhibited exceptional stability across all evaluated seeds and error formulations. The MAE yielded an average performance of $1.3802$ with a minimal standard deviation of $0.0008$ (ranging from a minimum of $1.3791$ to a maximum of $1.3814$). Similarly, the RMSE demonstrated high consistency, averaging $1.9729\pm 0.0008$ ($[1.9712,1.9737]$). Finally, the RMSSE confirmed these findings with an average value of $0.7117$ and a negligible standard deviation of $0.0003$ ($[0.7110,0.7120]$).

This marginal variance across all considered error norms indicates that the high-dimensional state space generated by the reservoir, in conjunction with Ridge regression ($\beta ={10}^{-1}$), achieves a robust and consistent representation of the underlying dynamics irrespective of the initial stochastic configuration. Consequently, the ESN architecture demonstrates substantial structural robustness, ensuring that the predictive performance of the reconstructed time series remains invariant to initialization fluctuations.

Furthermore, computational tracking demonstrates the efficiency and scalability of the underlying implementation. The validation procedure for the initial execution (Seed 0) required $33.61$ seconds, with subsequent evaluations maintaining an almost identical computational cost and stabilizing at an average execution time of $32.81$ seconds per complete validation sequence. This highly consistent processing time, devoid of significant initial compilation overheads, confirms the computational suitability and reliability of the proposed framework for deployment in resource-constrained edge architectures.

4. Discussion

The results obtained from both the static partition and the sequential WFV demonstrate that the ESN achieves competitive predictive accuracy, structural robustness, and computational efficiency when modeling non-stationary landing time series characterized by high variability and a high frequency of zero-value observations.

In the static evaluation scheme, the ESN achieved the lowest testing error among all evaluated architectures, recording an $\mathrm{MAE}=1.90\mathrm{t}$, an $\mathrm{RMSE}=2.67\mathrm{t}$, and an $\mathrm{RMSSE}=0.67$. These metrics reflect a substantial performance improvement over traditional linear statistical baselines (SARIMA: $\mathrm{MAE}=6.18\mathrm{t}$, $\mathrm{RMSE}=7.81\mathrm{t}$), standard tree-based ensemble methods (LightGBM: $\mathrm{MAE}=2.34\mathrm{t}$, $\mathrm{RMSE}=2.93\mathrm{t}$), and deep recurrent baselines (LSTM: $\mathrm{MAE}=6.23\mathrm{t}$, $\mathrm{RMSE}=6.69\mathrm{t}$). While linear and autoregressive baselines exhibited pronounced structural degradation during periods of rapid landing fluctuations, the recurrent dynamic structure of the ESN captured complex non-linear interactions without suffering from overfitting or phase-lag artifacts.

Under the sequential walk-forward validation scheme ($h=1$), where models were continuously updated with new temporal observations, the performance advantage of the ESN framework became more pronounced. The ESN attained an $\mathrm{MAE}=1.37\mathrm{t}$, an $\mathrm{RMSE}=1.96\mathrm{t}$, and an $\mathrm{RMSSE}=0.70$. An RMSSE below $1.0$ indicates that the ESN outperforms the baseline Naïve benchmark by approximately $30\%$, effectively handling data sparsity and volatile landing yields. This performance stems from the online adaptation capability of the linear readout layer via Ridge regression, which allows the network to track local temporal regime shifts without retraining the fixed internal reservoir weights.

Deep learning baselines such as the LSTM ($\mathrm{MAE}=1.92\mathrm{t}$, $\mathrm{RMSSE}=0.86$ under WFV) demonstrated competitive pattern recognition capabilities but required higher computational complexity and memory allocation ($19.29\mathrm{s}/\mathrm{fold}$). Conversely, while LightGBM executed rapidly during training, it failed to generate realistic continuous trajectories over multi-step or out-of-sample prediction horizons, frequently predicting negative landing values ($\mathrm{Landing}<0\mathrm{t}$). In contrast, the internal state dynamics of the ESN constrained predicted values within physically plausible domain boundaries ($\mathrm{Min}=0.77\mathrm{t}$), preserving the operational logic of the fishery landing series.

Regarding the pronounced discrepancy observed in tree-based ensemble baselines, such as the training-to-validation error divergence in LightGBM (training MAE of $9.47$ versus validation MAE of $1.96$, as detailed in Table 11, this behavior highlights a fundamental limitation of static machine learning models under structural breaks. Because the evaluation protocol establishes a strict temporal partition to test out-of-distribution generalization, the transition from the training period to validation and test horizons encompasses a shift in the underlying dynamical regime. Although standard regularization constraints (specified in Table 9) prevented excessive tree growth during training, LightGBM relies on axis-aligned split decisions over static feature spaces, rendering it incapable of extrapolating temporal trajectories when the statistical properties of the data change post-training. Consequently, the apparent overfitting is not merely a failure of hyperparameter tuning, but a structural vulnerability of tabular learning algorithms when confronted with non-stationary time series transitions. In contrast, the recurrent reservoir dynamics of the ESN absorb these regime shifts by projecting temporal histories into a high-dimensional continuous manifold.

Statistical significance testing confirmed the quantitative advantage of the proposed ESN architecture. Both the Diebold–Mariano test and the non-parametric Wilcoxon signed-rank test indicated statistically significant differences ($p<0.05$) in predictive error distributions when comparing the ESN against SARIMA and LightGBM on out-of-sample data. While the error difference between the ESN and the LSTM was narrower in certain folds, the ESN achieved competitive accuracy at a fraction of the computational and time cost during training and inference.

Analysis of the ESN bidirectional temporal reconstruction reveals a clear distinction between operational decline and catch-based variations within the fishery timeline. Specifically, the operational decline observed around 2005—characterized by the cessation of monitoring and a sharp reduction in industrial fleet activity—represented an economic and operational contraction driven by exogenous macroeconomic pressures, including fuel price fluctuations and fleet dynamics. As evidenced by the bidirectional ESN reconstruction, where the mean reconstructed landing remained at $39.98\mathrm{t}$, well above the landing-based operational threshold of $17.98\mathrm{t}$, positive landing levels were nevertheless sustained throughout the unmonitored window of 2000–2012. In contrast, the post-2013 trajectory constitutes a sustained catch-based decline. Governed by cumulative exploitation patterns and compounding environmental variability, landings persistently decreased toward near-zero values—falling below $4.54\mathrm{t}$ in the test set—thereby indicating the firm persistence of a low-yield regime within the reported records.

Notably, the slight upward inflection observed at the terminal end of the reconstructed series is consistent with recent independent biological evidence for the same stock. Paramo et al. [57] reported low exploitation rates ($E=0.19$–$0.50$) for pink shrimp (Penaeus notialis) in the Colombian Caribbean, attributing this to the drastic reduction of the industrial trawling fleet to approximately five vessels, and suggested that the stock could be undergoing a recovery phase. This convergence between a purely data-driven reconstruction and an independent stock assessment strengthens confidence that the terminal signal reflects a genuine biological trend rather than a modeling artifact.

Beyond numerical error metrics, the contrast in predictive performance across historical periods can be interpreted through the principles of non-linear dynamical systems. Framing the initial monitoring phase ($s_1$) as a volatile dynamic regime suggests a system governed by sensitivity to initial conditions, characterized by a positive maximal Lyapunov exponent (${\lambda }_{\max }>0$). Under this theoretical framing, the Lyapunov prediction horizon (${\tau }_L\approx 1/{\lambda }_{\max }$) is restricted, rendering autonomous forward forecasting across extensive multi-year gaps challenging without bidirectional soft constraints. Conversely, following the catch-based decline post-2013 ($s_3$), the system transitioned into a low-variance, stable state. This dynamical stabilization corresponds to a reduced divergence rate and an extended predictability horizon, explaining why predictive metrics stabilize during this era.

From an architectural standpoint, the Echo State Network mirrors these dynamical constraints. The reservoir weight matrix satisfies the ESP, which mathematically guarantees the Fading Memory Property by ensuring that the internal conditional Lyapunov exponents of the network remain negative. This intrinsic property forces the reservoir to exponentially decay the influence of distant past states, preventing distant noise from corrupting current state estimations while maintaining the short- to medium-term non-linear dynamics required for reliable bidirectional reconstruction.

Furthermore, the stability of the evaluation metrics across multiple random initializations (Table 15) provides empirical evidence of the structural robustness of the proposed ESN framework. In Reservoir Computing, stochastic weight generation (W and $W_{\mathrm{in}}$) can introduce performance variability across different random seeds. However, the consistent response observed in these sensitivity trials demonstrates that the combination of a high-dimensional state space ($N_x=500$), an appropriate spectral radius ($\rho =0.99$), and Ridge regularization ($\beta ={10}^{-1}$) projects the non-stationary landing dynamics into a stable manifold where the readout optimization converges reliably. This deterministic behavior indicates that the reconstructed subseries are governed by the intrinsic non-linear memory of the historical data, supporting the reliability of the model for data recovery.

Finally, despite the performance demonstrated under univariate fishery conditions, several inherent methodological limitations exist. First, the benchmark dataset was reconstructed using missForest prior to model comparison. Although this approach was necessary to provide a common complete dataset for SARIMA, LightGBM, and LSTM baselines, a potential source of bias may arise from the imputation stage itself. Nevertheless, only 5.6% of the observations in the independent test partition required reconstruction, indicating that the final out-of-sample evaluation was overwhelmingly supported by directly observed landing records.

Second, because the analytical framework relies exclusively on a univariate historical time series without exogenous covariates (such as environmental, climatic, or economic indicators), the capacity of the model to anticipate structural regime shifts driven entirely by external factors remains constrained. Second, regarding scalability and extreme non-stationarity, while a fixed reservoir size ($N_x=500$) proved adequate for capturing the temporal dependencies of this specific single-variable domain, scaling the architecture to handle high-dimensional multivariate inputs would require expanding the input weight matrix dimensions and re-tuning the spectral radius to prevent chaotic state saturation under severe distribution shifts. Future work will focus on integrating multi-source exogenous streams and optimizing adaptive reservoir scaling algorithms for multi-variable edge deployments.

5. Conclusions

Under small-data regimes, RC implemented via an ESN demonstrated competitive forecasting accuracy compared with an autoregressive baseline (SARIMA), a gradient-boosted decision tree ensemble (LightGBM), and a conventional recurrent neural network (LSTM). Evaluated under a strictly out-of-sample WFV protocol ($h=1$), the ESN configuration yielded a MAE of $1.37$, a RMSE of $1.96$, and a RMSSE of $0.70$, reducing forecasting error by $30\%$ relative to the non-seasonal naive baseline ($\mathrm{RMSSE}=1.00$). The consistency of predictive metrics across temporal validation partitions indicates that the non-linear dynamic reservoir captures complex temporal patterns without exhibiting catastrophic overfitting.

From the perspective of edge intelligence deployed on embedded hardware (NVIDIA Jetson Orin Nano), the experimental benchmark established a Pareto frontier among predictive accuracy, latency, and memory footprint. Although LightGBM proved to be the least resource-intensive model, with an inference latency of $0.07\mathrm{ms}/\mathrm{sample}$ and a system RAM usage of $70.01\%$, the ESN provided a competitive trade-off by recording sub-millisecond inference latency ($0.27\mathrm{ms}$) and a training time of $424.44\mathrm{ms}$ per fold. Compared with the high latency of the LSTM ($19.29\mathrm{s}/\mathrm{fold}$) caused by BPTT and graph recompilation, as well as the high operational cost of SARIMA ($60.39\mathrm{s}/\mathrm{fold}$ with a RAM accumulation of $631.91\mathrm{MB}$), the ESN updates its weights through an exact linear ridge regression.

Although average instantaneous power draw remains comparable across model families, the LSTM consumes approximately 43 times more total energy than the ESN to complete the WFV task. This energy footprint indicates that the ESN is a more energy-efficient solution for autonomous edge monitoring systems powered by batteries or photovoltaic panels.

The high-dimensional state space generated by the fixed reservoir allows the network to map complex non-linear temporal dynamics—such as industrial shrimp landings—even in the presence of continuous missing data gaps. Operating on univariate time series, the bidirectional ESN framework provides historical reconstruction and enables identification of the post-2013 landing decline threshold, with the terminal recovery signal consistent with independent stock-level evidence for the same fishery [57].

The ESN implementation in JAX using 64-bit precision demonstrates that over-parameterized deep learning models or linear autoregressive methods are unnecessary for edge time-series forecasting. Projecting input dynamics into a high-dimensional fixed reservoir while updating only the linear readout layer constitutes an efficient strategy for continuous monitoring under hardware and data constraints.