| Title: | Hybrid EMD-SSA Forecasting with Artificial Neural Networks for Time Series |
|---|---|
| Description: | Implements a hybrid empirical mode decomposition, singular spectrum analysis, and artificial neural network forecasting framework for time series data. The framework extracts the trend from the first intrinsic mode function using singular spectrum analysis, forecasts the extracted trend and remaining decomposition components using artificial neural networks, and combines the component forecasts to obtain the final forecast. For methodological details, refer to Kumar et al. (2026) <doi:10.1007/s10791-026-10076-2>. |
| Authors: | Prabhat Kumar [aut, cre], Girish Kumar Jha [aut], Rajeev Ranjan Kumar [aut], Santosha Rathod [aut], R N Singh [aut] |
| Maintainer: | Prabhat Kumar <prabhatkkv@gmail.com> |
| License: | MIT + file LICENSE |
| Version: | 0.1.0 |
| Built: | 2026-10-15 20:34:39 UTC |
| Source: | https://github.com/cran/EMDSSAForecast |
Calculate common forecast accuracy metrics including RMSE, MAE, MAPE, and R-squared.
accuracy_metrics(actual, forecast)accuracy_metrics(actual, forecast)
actual |
Numeric vector containing the observed values. |
forecast |
Numeric vector containing the forecasted values. |
Missing values are removed pairwise before calculating the metrics.
Observations with an actual value equal to zero are excluded from
the MAPE calculation to avoid division by zero. If all actual
values are zero, MAPE is returned as NA.
R-squared is calculated as
A one-row data frame containing:
Root Mean Squared Error.
Mean Absolute Error.
Mean Absolute Percentage Error in percent.
Coefficient of determination (R-squared).
# Example using the built-in AirPassengers dataset actual <- as.numeric(AirPassengers) # Illustrative forecasts forecast <- actual + 1 accuracy_metrics( actual = actual, forecast = forecast )# Example using the built-in AirPassengers dataset actual <- as.numeric(AirPassengers) # Illustrative forecasts forecast <- actual + 1 accuracy_metrics( actual = actual, forecast = forecast )
Fit a feed-forward artificial neural network to lagged time-series observations and produce recursive multi-step forecasts.
annTSforecast( series, h = 1, lag = 5, hidden = c(70, 20), maxit = 2500, seed = 123 )annTSforecast( series, h = 1, lag = 5, hidden = c(70, 20), maxit = 2500, seed = 123 )
series |
Numeric univariate training time series. |
h |
Forecast horizon. Default is 1. |
lag |
Number of lagged observations used as ANN inputs. Default is 5. |
|
Hidden-layer sizes. The default configuration is
|
|
maxit |
Maximum number of neural-network training steps. Default is 2500. |
seed |
Random seed for reproducible model fitting. Default is 123. |
The function converts the univariate time series into a supervised learning dataset using lagged observations. The input variables and response are scaled to the interval from 0 to 1 before fitting the artificial neural network. The fitted network is then used to generate recursive multi-step forecasts, where each predicted value is appended to the historical series and used as an input for subsequent forecasts.
A list containing:
In-sample fitted values. The first lag observations
are returned as NA.
Recursive multi-step forecasts.
Fitted neural network object.
Number of lagged observations used as inputs.
Hidden-layer configuration.
Maximum number of training steps.
Random seed used for model fitting.
Minimum value used for min-max scaling.
Maximum value used for min-max scaling.
# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) fit <- annTSforecast( series = series, h = 3, lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) fit$forecast# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) fit <- annTSforecast( series = series, h = 3, lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) fit$forecast
Decompose a univariate time series into intrinsic mode functions (IMFs)
and a residual component using Rlibeemd::emd().
emd_decompose(train, S_number = 4, num_siftings = 50)emd_decompose(train, S_number = 4, num_siftings = 50)
train |
Numeric univariate time series. |
S_number |
S-number stopping criterion used by
|
num_siftings |
Maximum number of sifting iterations. Default is 50. |
Empirical Mode Decomposition (EMD) decomposes a time series into a set of intrinsic mode functions and a residual component. The resulting components can subsequently be used for further analysis or hybrid forecasting.
A numeric matrix whose columns contain IMF1, IMF2, ..., and Residual. Each row corresponds to an observation in the original time series.
# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) # Perform EMD decomposition result <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) # View the dimensions and component names dim(result) colnames(result)# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) # Perform EMD decomposition result <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) # View the dimensions and component names dim(result) colnames(result)
Run the complete hybrid forecasting workflow: EMD -> SSA trend extraction from IMF1 -> ANN component forecasting -> ensemble reconstruction.
emd_ssa_annTSforecast( train, h = 1, L = 12, trend_rate = 0.9, lag = 5, hidden = c(70, 20), maxit = 2500, S_number = 4, num_siftings = 50, seed = 123 )emd_ssa_annTSforecast( train, h = 1, L = 12, trend_rate = 0.9, lag = 5, hidden = c(70, 20), maxit = 2500, S_number = 4, num_siftings = 50, seed = 123 )
train |
Numeric univariate training time series. |
h |
Forecast horizon. Default is 1. |
L |
SSA window length. Default is 12. |
trend_rate |
Proportion of SSA variance used for trend extraction. Default is 0.90. |
lag |
Number of lagged observations used as ANN inputs. Default is 5. |
|
Hidden-layer sizes for the ANN. Default is |
|
maxit |
Maximum number of neural-network training steps. Default is 2500. |
S_number |
S-number stopping criterion for EMD. Default is 4. |
num_siftings |
Maximum number of EMD sifting iterations. Default is 50. |
seed |
Random seed for reproducible ANN fitting. Default is 123. |
The hybrid forecasting framework consists of four main stages. First, Empirical Mode Decomposition (EMD) decomposes the training series into intrinsic mode functions and a residual component. Second, Singular Spectrum Analysis (SSA) is applied to IMF1 to extract its trend component. Third, the reconstructed components are forecast using an artificial neural network based on lagged observations. Finally, the individual component forecasts are aggregated to obtain the overall forecast.
An object of class emd_ssa_annTSforecast containing:
The matched function call.
The training time series.
Forecast horizon.
EMD components including the residual.
SSA trend extraction results.
SSA-extracted trend component.
Overall in-sample fitted values.
Overall forecasts.
Fitted values for individual components.
Forecasts for individual components.
List containing all model parameters.
# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) fit <- emd_ssa_annTSforecast( train = series, h = 5, L = 12, trend_rate = 0.90, lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) # Display the forecasts fit$forecast # Display the model summary fit # Plot the training series and forecasts plot(fit)# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) fit <- emd_ssa_annTSforecast( train = series, h = 5, L = 12, trend_rate = 0.90, lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) # Display the forecasts fit$forecast # Display the model summary fit # Plot the training series and forecasts plot(fit)
Plot the intrinsic mode functions (IMFs) and residual obtained from empirical mode decomposition.
plot_imfs(imfs, main = "EMD Decomposition")plot_imfs(imfs, main = "EMD Decomposition")
imfs |
Numeric matrix containing EMD components. Columns represent individual IMFs and the residual. |
main |
Main title of the plot. Default is "EMD Decomposition". |
Each EMD component is displayed in a separate panel. Different colours are used to distinguish the components visually.
The function produces a multi-panel plot and invisibly returns NULL.
# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) imfs <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) plot_imfs(imfs)# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) imfs <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) plot_imfs(imfs)
Plot the original IMF together with the trend extracted using singular spectrum analysis (SSA).
plot_ssa_trend(imf, trend, main = "SSA Trend Extraction from IMF1")plot_ssa_trend(imf, trend, main = "SSA Trend Extraction from IMF1")
imf |
Numeric vector containing the original IMF used for SSA trend extraction. |
trend |
Numeric vector containing the SSA-extracted trend. |
main |
Main title of the plot. Default is "SSA Trend Extraction from IMF1". |
The original IMF is displayed in grey and the SSA-extracted trend is displayed in red. This visualization allows the extracted trend to be compared directly with the corresponding IMF.
The function produces a plot and invisibly returns NULL.
series <- as.numeric(AirPassengers) imfs <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) result <- ssa_trend( imf = imfs[, "IMF1"], L = 12, trend_rate = 0.90 ) plot_ssa_trend( imf = imfs[, "IMF1"], trend = result$trend )series <- as.numeric(AirPassengers) imfs <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) result <- ssa_trend( imf = imfs[, "IMF1"], L = 12, trend_rate = 0.90 ) plot_ssa_trend( imf = imfs[, "IMF1"], trend = result$trend )
Plot the training series together with forecasts obtained from the hybrid EMD-SSA-ANN forecasting model.
## S3 method for class 'emd_ssa_annTSforecast' plot(x, ...)## S3 method for class 'emd_ssa_annTSforecast' plot(x, ...)
x |
An object of class |
... |
Additional graphical arguments. |
Invisibly returns the supplied model object.
series <- as.numeric(AirPassengers) fit <- emd_ssa_annTSforecast( train = series, h = 5, L = 12, trend_rate = 0.90, lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) plot(fit)series <- as.numeric(AirPassengers) fit <- emd_ssa_annTSforecast( train = series, h = 5, L = 12, trend_rate = 0.90, lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) plot(fit)
Forecast the SSA-extracted IMF1 trend and each remaining EMD component using artificial neural networks, then combine their forecasts to obtain the final hybrid forecast.
reconstruct_forecast( imfs, trend, h, lag = 5, hidden = c(70, 20), maxit = 2500, seed = 123 )reconstruct_forecast( imfs, trend, h, lag = 5, hidden = c(70, 20), maxit = 2500, seed = 123 )
imfs |
Numeric matrix or data frame containing EMD components,
typically returned by |
trend |
Numeric vector containing the SSA-extracted trend from IMF1. |
h |
Forecast horizon. |
lag |
Number of lagged observations used as ANN inputs. Default is 5. |
|
Hidden-layer sizes for the ANN. Default is |
|
maxit |
Maximum number of neural-network training steps. Default is 2500. |
seed |
Random seed for reproducible ANN fitting. Default is 123. |
The SSA-extracted trend from IMF1 is forecast separately using an artificial neural network. Each remaining EMD component, including the residual component, is also forecast separately. The component forecasts are then summed to obtain the final hybrid forecast.
A list containing:
Final reconstructed in-sample fitted values.
Final reconstructed forecasts.
A list containing fitted values for the SSA-extracted trend and individual EMD components.
A list containing forecasts for the SSA-extracted trend and individual EMD components.
A list containing the fitted ANN model objects for each component.
# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) # Perform EMD decomposition imfs <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) # Extract the trend from IMF1 using SSA ssa <- ssa_trend( imfs[, "IMF1"], L = 12, trend_rate = 0.90 ) # Forecast and reconstruct the components result <- reconstruct_forecast( imfs = imfs, trend = ssa$trend, h = 5, lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) # Display the final forecasts result$forecast# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) # Perform EMD decomposition imfs <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) # Extract the trend from IMF1 using SSA ssa <- ssa_trend( imfs[, "IMF1"], L = 12, trend_rate = 0.90 ) # Forecast and reconstruct the components result <- reconstruct_forecast( imfs = imfs, trend = ssa$trend, h = 5, lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) # Display the final forecasts result$forecast
Apply Singular Spectrum Analysis (SSA) to an intrinsic mode function and reconstruct the minimum number of eigentriples required to reach a specified Trend Rate.
ssa_trend(imf, L = 12, trend_rate = 0.9)ssa_trend(imf, L = 12, trend_rate = 0.9)
imf |
Numeric IMF series, normally IMF1 obtained from
|
L |
SSA window length. Default is 12. |
trend_rate |
Proportion of the cumulative SSA contribution retained for trend reconstruction. Must be between 0 and 1. Default is 0.90. |
Singular Spectrum Analysis is applied to the supplied IMF series.
The relative contribution of each SSA component is calculated from
the singular values, and the minimum number of leading eigentriples
whose cumulative contribution reaches the specified trend_rate
is selected for reconstruction.
A trend_rate of 0 returns a zero-valued trend, corresponding
to complete removal of the supplied IMF. A trend_rate of 1
returns the original IMF unchanged.
A list containing:
Reconstructed SSA trend component.
Fitted SSA model.
SSA window length.
Specified Trend Rate.
Number of eigentriples retained for trend reconstruction.
SSA eigenvalues.
Relative contribution of the SSA components.
Cumulative relative contribution of the SSA components.
# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) # Perform EMD decomposition imfs <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) # Extract the trend from IMF1 using SSA result <- ssa_trend( imf = imfs[, "IMF1"], L = 12, trend_rate = 0.90 ) # Display the extracted trend result$trend # Display the number of retained eigentriples result$eigentriples# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) # Perform EMD decomposition imfs <- emd_decompose( train = series, S_number = 4, num_siftings = 50 ) # Extract the trend from IMF1 using SSA result <- ssa_trend( imf = imfs[, "IMF1"], L = 12, trend_rate = 0.90 ) # Display the extracted trend result$trend # Display the number of retained eigentriples result$eigentriples
Evaluate combinations of SSA window lengths and Trend Rate values using a supplied hold-out test set.
tune_emd_ssa_annTSforecast( train, test, L_values = c(12, 36, 48), trend_rates = c(0.85, 0.9, 0.95), lag = 5, hidden = c(70, 20), maxit = 2500, S_number = 4, num_siftings = 50, seed = 123 )tune_emd_ssa_annTSforecast( train, test, L_values = c(12, 36, 48), trend_rates = c(0.85, 0.9, 0.95), lag = 5, hidden = c(70, 20), maxit = 2500, S_number = 4, num_siftings = 50, seed = 123 )
train |
Numeric training time series. |
test |
Numeric hold-out observations used for evaluating forecasts. |
L_values |
Numeric vector of SSA window lengths to evaluate. |
trend_rates |
Numeric vector of Trend Rate values to evaluate. Each value must be between 0 and 1. |
lag |
Number of lagged observations used as ANN inputs. Default is 5. |
|
Hidden-layer sizes for the ANN. Default is |
|
maxit |
Maximum number of neural-network training steps. Default is 2500. |
S_number |
S-number stopping criterion for EMD. Default is 4. |
num_siftings |
Maximum number of EMD sifting iterations. Default is 50. |
seed |
Random seed for reproducible ANN fitting. Default is 123. |
The function first decomposes the training series using EMD. For each combination of SSA window length and Trend Rate, SSA is applied to IMF1 to extract the trend component. The resulting hybrid EMD-SSA-ANN model is then used to forecast the complete hold-out period. Forecast accuracy is evaluated using RMSE, MAE, MAPE, and R-squared.
A data frame containing the evaluated combinations of SSA window length and Trend Rate, together with the forecast horizon and accuracy measures RMSE, MAE, MAPE, and R-squared. The returned data frame is ordered by RMSE, MAE, and MAPE.
# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) # Use the first 132 observations for training and the final # 12 observations as the hold-out test set train <- series[1:132] test <- series[133:144] tuning <- tune_emd_ssa_annTSforecast( train = train, test = test, L_values = c(12, 24), trend_rates = c(0.90, 0.95), lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) tuning# Example using the built-in AirPassengers dataset series <- as.numeric(AirPassengers) # Use the first 132 observations for training and the final # 12 observations as the hold-out test set train <- series[1:132] test <- series[133:144] tuning <- tune_emd_ssa_annTSforecast( train = train, test = test, L_values = c(12, 24), trend_rates = c(0.90, 0.95), lag = 5, hidden = c(10, 5), maxit = 500, seed = 123 ) tuning