Sunspots: NNS vs. KERAS
Intro
We are comparing RMSE of eleven 10 year predictions (120 monthly
observations) for the sunspot.month dataset in R using
KERAS LSTM and NNS.ARMA.
Long Short-Term Memory (LSTM) Models
LSTMs are quite useful in time series prediction tasks involving autocorrelation, the presence of correlation between the time series and lagged versions of itself, because of their ability to maintain state and recognize patterns over the length of the time series.
In normal (or stateless) mode, KERAS
shuffles the samples, and the dependencies between the time series and
the lagged version of itself are lost. However, when run in
stateful mode, we can often get high accuracy results by
leveraging the autocorrelations present in the time series.
KERAS LSTM description and example with accompanying
code (lots of it!) available via the following link: https://www.business-science.io/timeseries-analysis/2018/04/18/keras-lstm-sunspots-time-series-prediction.html
NNS.ARMA
NNS in its forecasting routine NNS.ARMA
also maintains the dependencies between the time series and lagged
values of itself, and then uses these dependencies in a regression
(either linear or nonlinear).
See here for a thorough description: https://www.researchgate.net/publication/327495856_Forecasting
Install the Latest Version of NNS (>= 11.6.3)
#library(devtools); install_github('OVVO-Financial/NNS', ref = "NNS-Beta-Version")
library(NNS)
library(Metrics)
library(xts)
library(zoo)
library(doParallel)Step 1: Create the Subsets of sunspot.month
We need to align dates to the KERAS example, so we
create an xts object to manipulate the dates and verify
they match.
sunspots_xts <- as.xts(sunspot.month)
tmp <- tempfile()
write.zoo(sunspots_xts, sep = ",", file = tmp)
sun <- read.zoo(tmp, sep = ",", FUN = as.yearmon)
sun_xts <- as.xts(sun)
# Create the same end date for all slices
dates=c("/18081202", "/18290102", "/18490202", "/18690302", "/18890402",
"/19090502","/19290602","/19490702","/19690802","/19890902","/20091002")
Slice=list()
for(i in 1:11){
Slice[[paste0("Slice_",i)]] = sun_xts[dates[i]]
print(tail(Slice[[i]],1))
}## [,1]
## Dec 1808 12.3
## [,1]
## Jan 1829 43
## [,1]
## Feb 1849 131.8
## [,1]
## Mar 1869 52.7
## [,1]
## Apr 1889 4.3
## [,1]
## May 1909 36
## [,1]
## Jun 1929 71.9
## [,1]
## Jul 1949 125.8
## [,1]
## Aug 1969 98
## [,1]
## Sep 1989 176.7
## [,1]
## Oct 2009 4.8
Step 2: Determine the seasonality
Here’s the general procedure on a single Slice.
We use a modulo 12 (%%12) call on the generated seasonal
periods from NNS.seas to find the nearest logical annual
cycle data point.
No test set leakage into the detected seasonal periods as
we are eliminating the last 120 observations from the variable and store
it in the new variable training.
training = Slice[[6]][1:(length(Slice[[6]])-120)]
periods = NNS.seas(training, modulo = 12)$periods
head(periods)## [1] 684 876 564 804 420 588Step 3: Determine the Optimal Parameters
The NNS.ARMA.optim routine checks various parameter
combinations in order to best fit the training set according to any
objective function specified. We store this under variable
b. NNS.ARMA.optim returns the optimal seasonal
periods, the objective function result, and which NNS
regression method was used: linear, nonlinear, or a combination of
both.
We also specify the objective function via
obj.fn = expression(Metrics::rmse(actual, predicted)) and
the objective to minimize it. Any objective function can be used,
calling the specific terms actual and
predicted within the expression(...) call.
start.time = Sys.time()
nns.predict = NNS.ARMA.optim(variable = training,
h = 120,
seasonal.factor = periods,
obj.fn = expression(Metrics::rmse(actual, predicted)),
objective = "min",
print.trace = TRUE)$results## [1] "CURRNET METHOD: lin"
## [1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
## [1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 132 ) ...)"
## [1] "CURRENT lin OBJECTIVE FUNCTION = 27.162508733888"
## [1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 132, 276 ) ...)"
## [1] "CURRENT lin OBJECTIVE FUNCTION = 24.6431883515214"
## [1] "BEST method = 'lin', seasonal.factor = c( 132, 276 )"
## [1] "BEST lin OBJECTIVE FUNCTION = 24.6431883515214"
## [1] "CURRNET METHOD: nonlin"
## [1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
## [1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 132, 276 ) ...)"
## [1] "CURRENT nonlin OBJECTIVE FUNCTION = 35.3808838776509"
## [1] "BEST method = 'nonlin' PATH MEMBER = c( 132, 276 )"
## [1] "BEST nonlin OBJECTIVE FUNCTION = 35.3808838776509"
## [1] "CURRNET METHOD: both"
## [1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
## [1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 132, 276 ) ...)"
## [1] "CURRENT both OBJECTIVE FUNCTION = 27.9810674881252"
## [1] "BEST method = 'both' PATH MEMBER = c( 132, 276 )"
## [1] "BEST both OBJECTIVE FUNCTION = 27.9810674881252"Sys.time() - start.time## Time difference of 40.21502 secsStep 4: Calculate NNS RMSE
The performance measure used is the root mean squared error (RMSE)
against the last 120 observations of that particular
Slice.
Metrics::rmse(predicted = nns.predict, actual = tail(Slice[[6]], 120))## [1] 24.93657Step 5: Compare NNS RMSE with KERAS LSTM
RMSE
Below is an image of the last 10 years prediction of
Slice 11. NNS is significantly more
accurate…
Step 6: Let’s try all the slices…
We compare KERAS implementation of its rolling
forecasts, to NNS forecasts per Slice.
KERAS generates \(\mu_{RMSE}=34.4\) and \(\sigma_{RMSE} = 13.0\).
start.time=Sys.time()
NNS.RMSE = list()
# Run in parallel
cl <- makeCluster(detectCores()-1)
registerDoParallel(cl)
NNS.RMSE <- foreach(i = 1:11,.packages=c("NNS", "Metrics") ) %dopar% {
training = as.vector(Slice[[i]][1:(length(Slice[[i]])-120)])
# Seasonality per slice
periods = NNS.seas(training, modulo = 12, plot = FALSE)$periods
# Determine optimal parameters
nns.predictions = NNS.ARMA.optim(variable = training,
h = 120,
seasonal.factor = periods,
obj.fn = expression(Metrics::rmse(actual, predicted)),
objective = "min",
ncores = 1,
print.trace = FALSE)$results
# RMSE call
Metrics::rmse(predicted = nns.predictions, actual = tail(Slice[[i]], 120))
}
stopCluster(cl)
registerDoSEQ()
print(Sys.time()-start.time)## Time difference of 4.545993 minsResults & Comments
Per Slice results:
"NNS.RMSE"=as.vector(unlist(NNS.RMSE))
"KERAS.RMSE"=c(48.2,17.4,41,26.6,22.2,49,18.1,54.9,28,38.4,34.2)
print(as.data.frame(cbind(NNS.RMSE, KERAS.RMSE, "NNS % Improvement"=1-NNS.RMSE/KERAS.RMSE),
row.names = names(Slice)))## NNS.RMSE KERAS.RMSE NNS % Improvement
## Slice_1 23.20313 48.2 0.51860729
## Slice_2 26.36671 17.4 -0.51532841
## Slice_3 31.80686 41.0 0.22422288
## Slice_4 30.10533 26.6 -0.13177926
## Slice_5 34.32120 22.2 -0.54599999
## Slice_6 24.93657 49.0 0.49109032
## Slice_7 26.22320 18.1 -0.44879550
## Slice_8 67.10959 54.9 -0.22239690
## Slice_9 26.13460 28.0 0.06662145
## Slice_10 29.38657 38.4 0.23472487
## Slice_11 33.57341 34.2 0.01832141Aggregate results:
"Mean NNS.RMSE"=mean(NNS.RMSE)
"SD NNS.RMSE"=sd(NNS.RMSE)
rbind(`Mean NNS.RMSE`,
`Mean KERAS.RMSE`=mean(KERAS.RMSE),
`SD NNS.RMSE`,
`SD KERAS.RMSE`=sd(KERAS.RMSE))## [,1]
## Mean NNS.RMSE 32.10611
## Mean KERAS.RMSE 34.36364
## SD NNS.RMSE 12.15593
## SD KERAS.RMSE 12.99525NNS demonstrates both a significant reduction in RMSE
and a significant reduction in variance vs. KERAS LSTM over
11 samples of the sunspots data.
The parsimonious code and interpretability of parameter settings
clearly favors NNS as well. Run-times however, favor
KERAS significantly and that is solely a function of not
having comparable resources to Google! NNS is not optimized
/ parallelized which would enable even more combinatorial parameter
testing over the objective function space…this is a work in progress. So
to arguments that KERAS is not optimized…neither is
NNS!
If you have any related questions / comments, feel free to e-mail: ovvo.financial.systems@gmail.com
Thanks for your interest!