I’ve spoken to several young researchers lately that are enamored with FB’s Prophet.
It’s a great model, but two things stand out to me:
Philosophically, it’s more of a seasonal decomposition filter akin to ARIMA-SEAT-X13, than it is a forecasting model
Despite it’s great power, it’s not a panacea.
Let’s play a little game by throwing Prophet at the task of forecasting daily stock returns. I’ll use the IVV. We can quibble about the sample period and ticker choice, etc… Acknowledged. But, for my young analysts, out there, take a look at what happens when you treat Prophet like a cure-all for all your forecasting ills.
Housekeeping
cat("\014") # clear console
rm(list=ls()) # Clear the workspacelibrary(quantmod)
Loading required package: xts
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Loading required package: TTR
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
library(prophet)
Warning: package 'prophet' was built under R version 4.5.3
Loading required package: Rcpp
Loading required package: rlang
library(dplyr)
######################### Warning from 'xts' package ##########################
# #
# The dplyr lag() function breaks how base R's lag() function is supposed to #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
# source() into this session won't work correctly. #
# #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
# dplyr from breaking base R's lag() function. #
# #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
# #
###############################################################################
Attaching package: 'dplyr'
The following objects are masked from 'package:xts':
first, last
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Task: Create train, validation, and test periods, in that temporal order. Save the last 5 periods for testing, the 10 periods prior to the test period to be validation and the remaining to be training.
Model RMSFE
1 Prophet Full 0.007567685
2 Mean 0.007498113
Task: Create a plot of the last 30 days of the full training period with the test period spliced onto that. Then create a plot of the actual data, prophet forecasts, and a simple full training period sample mean as a benchmark.
window_size <-30train_tail <-tail(train_full, window_size)cutoff_date <-min(test$ds)# ---- Extract fitted values (in-sample) ----fit_full <-tail(fc_full$yhat[1:nrow(train_full)], window_size)# ---- Mean (same in-sample and out-of-sample) ----mean_line_train <-rep(mean_ret, window_size)mean_line_test <-rep(mean_ret, test_h)# ---- Combine data ----plot_df <-data.frame(ds =c(train_tail$ds, test$ds),actual =c(train_tail$y, test$y),prophet_full =c(fit_full, pred_full),mean =c(mean_line_train, mean_line_test))plot_long <-pivot_longer(plot_df, -ds)ggplot(plot_long, aes(x = ds, y = value, color = name)) +# Shade test regionannotate("rect",xmin = cutoff_date,xmax =max(plot_df$ds),ymin =-Inf,ymax =Inf,alpha =0.1) +# Linesgeom_line() +# 🔻 Train/Test splitgeom_vline(xintercept = cutoff_date, linetype ="dashed") +ggtitle("In-Sample Fit vs Out-of-Sample Forecast") +labs(color ="Series", x ="", y ="Return (1% = .01)" ) +theme_minimal()
EDA
Q: What did we do wrong? What should we always do first?
A: EDA. We SHOULD have done EDA first.
Task: plot the IVV over time
ggplot(df, aes(x = ds, y = y)) +geom_line() +ggtitle("IVV Daily Log Returns") +labs(x ="", y ="Return") +theme_minimal()
Q: Comment on the behavior as it relates to our modeling efforts.
A:
Mean 0
No trend
Little persistence
ACF
Task: Test for independence of the returns
Box.test(train_full$y, lag =10, type ="Ljung-Box")
There is little visible seasonality and not a strong linear time trend
Prophet relies on trend and seasonality. Since both are weak in returns, the model effectively collapses toward predicting the mean.
Minor note: The gray bars indicate the scale (amplitude) of each component within its own panel. They are not directly comparable to the original data scale and do not measure importance.
Task: Compute strength of trend and strength of seasonal
Strength of Trend = max(0,1-variance(noise)/var(Trend + noise))
Strength of Seasonal = max(0,1-variance(noise)/var(Seasonal + noise))
Running model 1 of 32 \nRunning model 2 of 32 \nRunning model 3 of 32 \nRunning model 4 of 32 \nRunning model 5 of 32 \nRunning model 6 of 32 \nRunning model 7 of 32 \nRunning model 8 of 32 \nRunning model 9 of 32 \nRunning model 10 of 32 \nRunning model 11 of 32 \nRunning model 12 of 32 \nRunning model 13 of 32 \nRunning model 14 of 32 \nRunning model 15 of 32 \nRunning model 16 of 32 \nRunning model 17 of 32 \nRunning model 18 of 32 \nRunning model 19 of 32 \nRunning model 20 of 32 \nRunning model 21 of 32 \nRunning model 22 of 32 \nRunning model 23 of 32 \nRunning model 24 of 32 \nRunning model 25 of 32 \nRunning model 26 of 32 \nRunning model 27 of 32 \nRunning model 28 of 32 \nRunning model 29 of 32 \nRunning model 30 of 32 \nRunning model 31 of 32 \nRunning model 32 of 32 \n
Task: Forecast this “best model”. Specifically, use the make_future_dataframe to create an object to store our forecasts. Then generate all of the in sample / out of sample forecasts and store within fc_best. Finally, grab only the out of sample and store in pred_best.
future_best <-make_future_dataframe(m_best, periods = test_h)fc_best <-predict(m_best, future_best)pred_best <-tail(fc_best$yhat, test_h)
Task: Create a simple baseline mean of the target variable as the forecast
