0
$\begingroup$

I have used the TBATS model on my data and when I apply the forecast() function, it automatically forecasts two years in the future. I haven't specified any training set or testing set, so how do I know how to calculate RMSE, MASE, MAE and MAPE?

The data I'm dealing with is Uber travel times data from Jan 2016 to Jan 2020. I have daily data (sampling frequency = 1) for 18 cities and each city has a different sample size (they range from 1422 days to 1459 days).

Here is the Kaggle link to my datasets. I'd recommend using the "final_" datasets instead of the "formatted_" ones.

Here is a subsample:

Date                        MeanTravelTimeSeconds
1   2016-01-02              1441.684
2   2016-01-03              1597.952
3   2016-01-04              1326.496
4   2016-01-05              1355.872
5   2016-01-06              1451.740
6   2016-01-07              1525.991
7   2016-01-08              1657.299
8   2016-01-09              1727.375
9   2016-01-10              1703.772

Here is the example for data for London:

London

I have set the vector of travel times as an msts object, for it has multiple seasonality, which is used by the TBATS model.

Below is a graph with the fitted values and with the forecast values:

enter image description here

I mainly have 3 questions:

  1. How do I correctly calculate the KPIs for the forecast (RMSE, MASE, MAE, MSE)?
  2. What is the difference between finding out the accuracy of the red line to the black line, and the accuracy of the blue line?
  3. Is what I'm doing called an in-sample forecast or an out-of-sample forecast?

The values I'm getting for each are very low, here's my code:

data <- read.csv('C:/users/Datasets/Final Datasets/final_london.csv', TRUE, ",")
y <- msts(data$MeanTravelTimeSeconds, start=c(2016,1), seasonal.periods=c(7.009615384615385, 30.5, 91.3, 365.25))

fit <- tbats(y)
fc <- forecast(fit)

# RMSE
rmse <- sqrt(fit$variance)

# MAE
res <- residuals(fit)
mae <- mean(abs(res))

# MAPE
pt <- (res)/y
mape <- mean(abs(pt))

# MSE (Mean Squared Error)
mse <- mean(res^2)
$\endgroup$
5
  • 1
    $\begingroup$ please share your data. $\endgroup$ Dec 28, 2020 at 1:10
  • $\begingroup$ I made an edit above (to the post) where I include the link to the Kaggle dataset. There you have access to all cities and I'd recommend using the "final_" ones instead of the "formatted_" ones. $\endgroup$ Dec 28, 2020 at 11:21
  • $\begingroup$ use dput(data) and copy and paste the output from the console into your question please. $\endgroup$ Dec 28, 2020 at 11:46
  • 1
    $\begingroup$ Time series error can be a tricky thing. It’s very easy to confuse yourself Into thinking it’s a lot better than it really is. One of the ways that you can double check your error metric is To look at the ratio of absolute predicted difference to absolute lag difference. If your estimate is a lot larger than a one step lag as a predictor then it might not be that great of an estimator. $\endgroup$ Dec 28, 2020 at 12:40
  • $\begingroup$ @EngrStudent thanks! Bu my confusion is: when you say absolute difference, are you assuming I have split my data into train and test sets? Because what I am trying to figure out is how to evaluate the accuracy of my forecast for a future for which I don't even have data (2021 and 2022). If I did an out-of-sample forecast for 2020, I could compare the forecast for 2020 with real data for 2020. But I want to know how to measure the accuracy for a period in the not-yet-seen future. $\endgroup$ Dec 28, 2020 at 13:36

1 Answer 1

1
$\begingroup$

Strangely, I can not get the tsCV() function to not return NAs. At any rate, the accuracy measures of the forecast can be obtained using the accuracy() function. We can also split into a training and test set as follows:

 library(forecast)
# Read in the csv: 
final_london <- transform(read.csv('C:/users/Datasets/Final Datasets/final_london.csv', TRUE, ","),
                          Date = as.Date(Date, "%Y-%m-%d"))

# Expand the series to contain an entry for every day in the year: 
final_london_expanded <-
  transform(merge(transform(final_london, estimate = FALSE),
    data.frame(Date = with(final_london, seq.Date(
    from = as.Date(paste0(gsub("\\-.*", "", min(Date)), "-01-01")),
    to = max(Date), by = "days"))), by = "Date", all = TRUE),
    estimate = ifelse(is.na(estimate), TRUE, estimate))

# Interpolate the data for the missing days: 
final_london_expanded[] <- lapply(seq_len(ncol(final_london_expanded)), function(i){
    if(inherits(final_london_expanded[,i], "Date")){
      as.Date(final_london_expanded[,i])
    }else if(is.numeric(final_london_expanded[, i])){
      x <- final_london_expanded[,i, drop = TRUE]
      ifelse(is.na(x), approx(x, n = length(x))$y[which(is.na(x))], x)
    }else{
      final_london_expanded[,i]
    }
  }
)

# Isolate the year:
final_london_expanded$year <- as.integer(gsub("\\-.*", "", final_london_expanded$Date))


# Days in year counter: 
final_london_expanded$day_in_year_ctr <- suppressWarnings(
  with(final_london_expanded, ave(year, year, FUN = seq.int))
)

# Store the series start year: 
series_start_year <- min(final_london_expanded$year)

# Store the first day number of the minimum year: 
series_sy_start_day <- with(final_london_expanded, min(day_in_year_ctr[year == series_start_year]))

# Coerce to multiple seasonality time series object: 
y <-
  msts(
    final_london_expanded$MeanTravelTimeSeconds,
    start = c(series_start_year, series_sy_start_day),
    seasonal.periods = c(7.009615384615385, 30.5, 91.3, 365.25)
  )
# Chart the decomposition of the series: 
plot(decompose(y))

# Isolate the year in which approximately 80 % of the series is represented;
# but still leaves at least a year for the test set: 
year_eighty_perc <- with(final_london_expanded, year[round(length(y) *  .8, 0)])

# Determine the cut-off point of the test set: 
train_max_year <- ifelse(year_eighty_perc >= (max(final_london_expanded$year) - 1), 
                         year_eighty_perc - 1, year_eighty_perc)

# Find the final day in the year in which the test set is cut off: 
series_my_end_day <- with(final_london_expanded, max(day_in_year_ctr[year == train_max_year]))

# Determine the starting year and period for the test set: 
test_start_year <- train_max_year + 1

# Determine the day in year that the test set should begin: 
test_sy_start_day <- with(final_london_expanded, min(day_in_year_ctr[year == test_start_year]))

# Partition the data into a test and training set;
# a training set:  
train <- window(y, start = c(series_start_year, series_sy_start_day), 
                end = c(train_max_year, series_my_end_day))

# Testing set:  
test <- window(y, start = c(test_start_year, test_sy_start_day))

# Chart the decomposition of the series: 
plot(decompose(train))

# Apply function to fit the tbats model: 
fit <- tbats(train)

# Generate the forecast (to measure against the test set): 
fc <- forecast(fit, h = length(test))

# Create a lookup table in order to extract the ts date values:
date_lkp <- data.frame(ts_rep = time(y), Date = as.Date(final_london_expanded$Date),
  mean_tt = final_london_expanded$MeanTravelTimeSeconds)

# Store a data.frame denoting the interpolated values: 
date_values_interped <- transform(subset(final_london_expanded, estimate == TRUE, 
                               select = c("Date", "MeanTravelTimeSeconds")), 
                        Date = time(y)[match(Date, date_lkp$Date)])


# Chart the forecast, actuals and interpolated values: 
autoplot(fc, xlab = "Date", ylab = "Number of Trips") +
  autolayer(test, series = "Actuals (test) data set") + 
  autolayer(fc$mean, series = "Forecasts") + 
  geom_point(data = date_values_interped, 
    aes(x = Date, y = MeanTravelTimeSeconds, color = "Interpolated Values"))

# Calculate the accuracy metrics for the model on the trainin and test sets;
# https://stackoverflow.com/questions/11092536/forecast-accuracy-no-mase-with-two-vectors-as-arguments
# http://freerangestats.info/blog/2019/07/20/time-series-cv
# NOTE:  the accuracy measures for the test set are the ones that you want to report: 
accuracy(fc, test)

# Regenerate the forecast model using the full data, but the previously selected model:
fcast_model <- tbats(y, model = fit)

# Generate the forecast: 
fcast <- forecast(fcast_model)

# Chart the forecast: 
autoplot(fcast, xlab = "Date", ylab = "Number of Trips")

# Chart the forecast / data partitions: 
autoplot(fcast, xlab = "Date", ylab = "Number of Trips") + 
  autolayer(train, series = "Training Data") + 
  autolayer(test, series = "Test Data") + 
  autolayer(fcast, series = "Forecasts")  + 
  geom_point(data = date_values_interped, 
             aes(x = Date, y = MeanTravelTimeSeconds, color = "Interpolated Values"))

# Calculate the accuracy metrics of the forecast model (on the training (entire data) set): 
accuracy(fcast)
$\endgroup$
5
  • $\begingroup$ Thanks, but when I do forecast(fc), I get the following error message: Error in ce(actual, predicted) : argument "predicted" is missing, with no default. Meaning, I need to tell it what data I'm testing it against, which I don't have because I'm not dividing my data into train and test. Simply forecasting the future (to which I don't have data: 2021 and 2022). $\endgroup$ Dec 28, 2020 at 16:28
  • 1
    $\begingroup$ @LuizScheuer please see my revised answer above where accuracy metrics are calculated on the test set. $\endgroup$ Dec 29, 2020 at 1:17
  • $\begingroup$ The visualization is a way better one then the one I have because of the inclusion of the test set. Three problems I ran into: (1) you use 365 to specify the last day of the year, but my datasets are not always 365 days for a year. is there a way to specify that the end is the last data point with year 2018? (2) Number 1 results in an error when I calculate the accuracy(fc, test) because the fc used training data with e.g. n = 364 but the test data has n = 363. (3) when I run accuracy(fcast) it tells me: Error in ce(actual, predicted): argument "predicted" is missing, with no default. $\endgroup$ Dec 29, 2020 at 20:47
  • $\begingroup$ @LuizScheuer The errors you are seeing you are introducing when translating my code back into your environment. If you want to ensure every year is full and complete please check my code above which expands the series to have values for every day in the series (using linear interpolation). $\endgroup$ Dec 30, 2020 at 1:18
  • 1
    $\begingroup$ @LuizScheuer Please upvote and accept my answer if you have found it useful / it does what you require. $\endgroup$ Dec 30, 2020 at 8:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.