# Automated Fourier Term Selection in Batch Forecasting via Foreach

I'm trying to parallelize my batch forecasting via foreach, and I'm pretty comfortable doing this via auto.arima and xreg terms, including Fourier terms that are pre-determined/calculated when the batch forecasting is done sequentially via a regular for or a foreach() %do% loop.

However, I was inspired by this chapter in Prof. Hyndman's textbook for automated selection of optimal pair of Fourier terms, but I'm not sure how I would implement this in a parallelized foreach loop, as even if it's just one univariate series, output of the foreach loop would output to a matrix via one of the options in the .combine argument, or just a list without it, and determining the "best K" or fitted model wouldn't be exported outside of the loop, (or at least I think?).

This question focuses on a paralellized with Fourier terms, but they're already determined beforehand.

I haven't tried the fable package yet, mainly because I know parallelization is still slow comparable to sequential. Also an automated process of determining the optimal K value in fable hasn't been implemented, though has been proposed.

Here's a sample of how I would do it sequentially:

# Sequential version
library(forecast)
library(lubridate)

list_of_means <-  c(5189.18109,2241.73218,2947.44890,18.62958,32.49026,19.06078,18.75038,24.59941,25.63607,18.11636)

list_of_stdevs <- c(992.247149,474.268766,610.983238,4.776820,9.558881,4.081409,4.441914,7.110687,6.635868,5.090847)

h <- 52

no_series <- 10

series_length <- (h*3)

maxK <- 25

bestKs <- rep(NA,no_series)

sim_models <- vector("list",length=no_series)

sim_series <- matrix(NA,nrow=series_length,ncol=no_series)

set.seed(12345)
for(i in seq(no_series)) {
sim_series[,i] <- rnorm(series_length,list_of_means[i],list_of_stdevs[i])
}

sim_series <- ts(sim_series,frequency=(365.25/7),start=decimal_date(ymd("2016-01-31")))

for(i in seq(no_series)) {
bestfit <- list(aicc=Inf)
for(K in seq(maxK)) {
fit <- auto.arima(sim_series[,i],xreg=fourier(sim_series[,i],K=K),seasonal=FALSE)
if(fit[["aicc"]] < bestfit[["aicc"]]) {
bestfit[["aicc"]] <- fit[["aicc"]]
sim_models[[i]] <- fit
bestKs[i] <- K
}
}
}

sim_forecasts <- vector("list",length=no_series)

for(i in seq(no_series)) {
sim_forecasts[[i]] <- forecast(sim_models[[i]],xreg=fourier(sim_series[,i],K=bestKs[i],h=h),h=h)
}


But this implementation is rather slow.

Here's what I'm doing, without doing an automated process of determining the optimal K value in foreach:

# Parallelized Version without Fouriers
library(forecast)
library(lubridate)
library(doParallel)

list_of_means <-  c(5189.18109,2241.73218,2947.44890,18.62958,32.49026,19.06078,18.75038,24.59941,25.63607,18.11636)

list_of_stdevs <- c(992.247149,474.268766,610.983238,4.776820,9.558881,4.081409,4.441914,7.110687,6.635868,5.090847)

h <- 52

no_cores <- detectCores() - 1

no_series <- 10

series_length <- (h*3)

maxK <- 25

bestKs <- rep(NA,no_series)

sim_models <- vector("list",length=no_series)

sim_series <- matrix(NA,nrow=series_length,ncol=no_series)

set.seed(12345)
for(i in seq(no_series)) {
sim_series[,i] <- rnorm(series_length,list_of_means[i],list_of_stdevs[i])
}

sim_series <- ts(sim_series,frequency=(365.25/7),start=decimal_date(ymd("2016-01-31")))

cl <- makeCluster(no_cores)

registerDoParallel(cl)

sim_models <- foreach(i=seq(no_series),.packages=c("forecast")) %dopar% {
auto.arima(sim_series[,i],xreg=fourier(sim_series[,i],K=maxK),seasonal=FALSE)
}

sim_forecasts <- foreach(i=seq(no_series),.packages=c("forecast")) %dopar% {
forecast(sim_models[[i]],xreg=fourier(sim_series[,i],K=maxK,h=h),h=h)
}


Is there something I'm missing in terms of the implementation of foreach that would help me accomplish this?

• You might want to consider shifting this question over to stackoverflow given it is more of a programming question and you might find a more helpful audience over there... In saying that I don't think it is inappropriate here. – André.B Nov 14 '19 at 20:25
• Ahh, I didnt know. Should I delete the question and move it over to SO? – nmck160 Nov 14 '19 at 21:20

I'm an idiot; I could have just used a nested foreach loop via %:% to store each pair-wise combination of K-value and each of my series/models in a nested list.

By not using the .combine argument, things are left simpler (IMO) to have a nested list where each element of our models can be examined (i.e. checking coefficients, AICc/BIC values).

For those who want to know a solution, albeit a clunky, inelegant one, here's my solution:

library(forecast)
library(doParallel)
library(lubridate)

list_of_means <-  c(5189.18109,2241.73218,2947.44890,18.62958,32.49026,19.06078,18.75038,24.59941,25.63607,18.11636)

list_of_stdevs <- c(992.247149,474.268766,610.983238,4.776820,9.558881,4.081409,4.441914,7.110687,6.635868,5.090847)

h <- 52

no_cores <- detectCores() - 1

no_series <- 10

series_length <- (h*3)

maxK <- 25

sim_series <- matrix(NA,nrow=series_length,ncol=no_series)

#Reproducible rnorm()
set.seed(12345)
for(i in seq(no_series)) {
sim_series[,i] <- rnorm(series_length,list_of_means[i],list_of_stdevs[i])
}

#Converting simulated data to time-series
sim_series <- ts(sim_series,frequency=(365.25/7),start=decimal_date(ymd("2016-01-31")))

#Making clusters and registering our parallel backend
cl <- makeCluster(no_cores)
registerDoParallel(cl)

#Parallelized loop where the outer loop is each of the ARIMA models, each iteration of the inner loop is a K-value within [1,25] pertaining to our ARIMA model
sim_models <- foreach(i=seq(no_series)) %:%
foreach(K=seq(maxK),.packages=c("forecast")) %dopar% {
auto.arima(sim_series[,i],xreg=fourier(sim_series[,i],K=K),seasonal=FALSE)
}

#Making a matrix, where the rows are the K-values, the columns are each of our ARIMA models, the elements comprising of the AICc values of each pairwise combination of K and model.
aicc_matrix <-matrix(NA,nrow=maxK,ncol=no_series)
for(col in seq(no_series)) {
for(row in seq(maxK)) {
aicc_matrix[row,col] <- sim_models[[col]][[row]][["aicc"]]
}
}

#Creating an array to store the K-value which minimizes our AICc value for each of our models.
bestK <- rep(NA,length=no_series)
for(i in seq(no_series)) {
bestK[i] <- which(aicc_matrix[,i]==min(aicc_matrix[,i]))
}

#Creating a new list which will consist of only the ARIMA models where the K-value Fourier term has our AICc minimized.
sim_models_best <- vector("list",length=no_series)
for(i in seq(no_series)) {
sim_models_best[[i]] <- sim_models[[i]][[bestK[i]]]
}

#Creating forecasts based off the models we just created and their corresponding optimal K-value.
sim_forecasts <- foreach(i=seq(no_series),.packages=c("forecast")) %dopar% {
forecast(sim_models_best[[i]],xreg=fourier(sim_series[,i],K=bestK[i],h=h),h=h)
}


This is my first post, and I realize now (and based off of Andre B's suggestion) that this should have probably been posted on StackOverflow.

However I do not want to spam this question so I will leave it here (unless advised otherwise), as I believe it can still be relevant here and of some help to people wanting to speed up batch forecasts while not going through tedious testing to find the optimal K-value for pairs of Fourier terms for their weekly-based forecasts.