Introduction

John Mount of Win-Vector LLC recently posted a a very interesting example of tide prediction using known theory generating frequencies to model tides vs. alternative methods of generated frequencies and prediction thereof. Please give the post a full read to familiarize yourself with the problem:

http://www.win-vector.com/blog/2019/08/lord-kelvin-data-scientist/

Predict Tides Using NNS

This is an excellent application of the time-series forecasting method in the NNS R-package.

Step 1: Install the Latest Version of NNS (>= 11.0)

#library(devtools); install_github('OVVO-Financial/NNS', ref = "NNS-Beta-Version")
library(NNS)

Step 2: Read & Load the Variables, Create Train and Test Sets

Download tides.RDS from the following URL:

https://github.com/WinVector/Examples/blob/master/Tides/tides.RDS

Create the training set dtrain and the test set dtest using the dates provided in the original post. In the final step, we determine the length of the training set to optimize on.

tides <- readRDS('tides.RDS')

base_date_time =  as.POSIXct('2001/01/01 00:00', tz = "UTC")
first_date_time =  as.POSIXct('2019/06/01 00:00', tz = "UTC")
cut_date_time = as.POSIXct('2019/07/15 00:00', tz = "UTC")

dtrain <- tides[tides$dt<cut_date_time, , drop = FALSE]
dtest <- tides[tides$dt>=cut_date_time, , drop = FALSE]

training_length <- dim(dtrain)[1] - dim(dtest)[1]

Step 3: Determine Frequencies, Optimize & Forecast

3.1: Determine Frequencies

In this step we will ascertain the most relevant periods via the NNS.seas() function. To isolate just the periods, we call nns_periods$periods.

nns_start.time <- Sys.time()
nns_periods <- NNS.seas(dtrain$tide_feet, modulo = 240, mod.only = FALSE)$periods

head(nns_periods)
## [1] 84719 91920 95267 84711 84718 84710
# Total relevant periods...
length(nns_periods)
## [1] 65591

3.2: Optimize the Combination of Frequencies

Now we will utilize the nns_periods in the NNS.ARMA.optim function, as well as the number of forecast periods h.

We are only going to use the first 200 relevant periods. There are thousands more relevant periods that can be included but would require additional computational resources.

nns_periods <- nns_periods[nns_periods<44000]

arma_parameters <- NNS.ARMA.optim(variable = dtrain$tide_feet, 
                                  h = nrow(dtest),
                                  training.set = length(dtrain$tide_feet) - nrow(dtest), 
                                  pred.int = .95,
                                  seasonal.factor = nns_periods[1:200],
                                  print.trace = FALSE)
## Time difference of 19.9947 mins

3.3: Forecast

Finally, we will extract the results using these optimum parameters, already provided in our arma_parameters object along with the other following parameters for the NNS.ARMA forecast.

arma_parameters[1:6]
## $periods
## [1] 39246 22979 39247
## 
## $weights
## NULL
## 
## $obj.fn
## [1] 0.01529544
## 
## $method
## [1] "lin"
## 
## $shrink
## [1] FALSE
## 
## $nns.regress
## [1] FALSE
nns_estimates <- arma_parameters$results

Step 4: Evaluate the Results

Let’s see how we did. The R-squared1 between predicted and actual is the presented metric.

\[ R^2 = \frac{[\sum_{i=1}^n (y_i - \bar{y})(\hat{y_i} - \bar{y})]^2}{\sum_{i=1}^n (y_i - \bar{y})^2\sum_{i=1}^n (\hat{y_i} - \bar{y})^2}\]

(sum((nns_estimates - mean(dtest$tide_feet)) * (dtest$tide_feet - mean(dtest$tide_feet))) ^ 2) / (sum((dtest$tide_feet - mean(dtest$tide_feet)) ^ 2) * sum((nns_estimates - mean(dtest$tide_feet)) ^ 2))
## [1] 0.9534088
library(ggplot2)
ggplot(aes(x=dt), data=dtest) +
  geom_line(aes(y=tide_feet), color='blue', alpha=1) + 
  geom_line(aes(y=nns_estimates), color='black', alpha=0.5) +
  geom_line(aes(y=arma_parameters$lower.pred.int), color='red', alpha=0.25) +
  geom_line(aes(y=arma_parameters$upper.pred.int), color='red', alpha=0.25) +
  ggtitle("prediction (blue) superimposed on actuals on test")

ggplot_data <- data.frame(cbind(nns_estimates,dtest, arma_parameters$lower.pred.int, arma_parameters$upper.pred.int))

ggplot(aes(x=nns_estimates,y=tide_feet), data = ggplot_data) +
  geom_point(alpha=0.1) + 
  ggtitle("prediction versus actual on test")

SCUM Forecasts and Timings

In addition to the NNS forecasts, the SCUM method was employed to generate forecasts and compare timings.

library(forecast)
library(smooth)
## Loading required package: greybox
## Package "greybox", v2.0.3 loaded.
## This is package "smooth", v4.1.0
# Start timing for SCUM
scum_start.time <- Sys.time()

# Generate SCUM forecasts
h <- length(dtest$tide_feet)

# Helper function for safe forecasting
safe_forecast <- function(forecast_function, data, h) {
  tryCatch({
    forecast_function(data, h = h)$mean
  }, error = function(e) {
    message(paste("Error in", deparse(substitute(forecast_function)), ":", e$message))
    return(rep(NA, h)) # Return NA for each horizon in case of an error
  })
}

# Perform individual forecasts safely
ets_forecast <- safe_forecast(function(data, h) forecast(ets(data), h = h), dtrain$tide_feet, h)
arima_forecast <- safe_forecast(function(data, h) forecast(auto.arima(data), h = h), dtrain$tide_feet, h)
theta_forecast <- safe_forecast(function(data, h) forecast(thetaf(data), h = h), dtrain$tide_feet, h)
## Error in forecast_function : Please select a longer horizon when the forecasts are first computed
ces_forecast <- safe_forecast(function(data, h) forecast(ces(data), h = h), dtrain$tide_feet, h)

# Combine forecasts, excluding errors (NA)
forecasts_matrix <- rbind(ets_forecast, arima_forecast, theta_forecast, ces_forecast)
forecasts_matrix <- forecasts_matrix[complete.cases(forecasts_matrix), ] # Remove rows with NA
scum_forecast <- apply(forecasts_matrix, 2, median)


# End timing for SCUM
scum_end.time <- Sys.time()

# SCUM execution time
scum_time <- scum_end.time - scum_start.time
scum_time
## Time difference of 6.835971 mins
# SCUM evaluation
(sum((scum_forecast - mean(dtest$tide_feet)) * (dtest$tide_feet - mean(dtest$tide_feet))) ^ 2) / (sum((dtest$tide_feet - mean(dtest$tide_feet)) ^ 2) * sum((scum_forecast - mean(dtest$tide_feet)) ^ 2))
## [1] 2.294788e-07
# Plot SCUM forecasts vs actual
ggplot() +
  geom_line(aes(x = dtest$dt, y = dtest$tide_feet), color = 'blue', alpha = 1) +
  geom_line(aes(x = dtest$dt, y = scum_forecast), color = 'green', alpha = 0.7) +
  ggtitle("SCUM Forecasts vs Actual")

Comments

The NNS method produces a very good result without any knowledge of or access to the tide machine and significantly better than the R-squared of 0.81 obtained with FFT frequencies used in the Elastic Net Regularized Linear Regression in the original blog post. Could NNS perform better? Yes. NNS.ARMA and NNS.ARMA.optim is an active area of development.

Furthermore, a larger cloud based instance would be able to handle more relevant periods and more granular data than the 6-minute intervals would also help. How much better could NNS perform is an open question…

NNS is not a one-trick pony, as it has been demonstrated to excel in time-series forecasting, nonlinear continuous regressions, and provide solutions for econometric applications. See the following examples:

I look forward to further discussions and collaboration with those equally as passionate about these issues, and open to embracing alternative solutions. If you found this presentation interesting or useful, please feel free to reach out via e-mail:

Thanks for your interest!

  1. See Hayfield and Racine https://cran.r-project.org/web/packages/np/vignettes/np.pdf↩︎