Required Packages

These examples were performed using NNS 3.7, currently only available on CRAN.

require(NNS);require(prophet);require(dplyr);require(datasets);require(quantmod)

Retail Sales Example

Below is the retail sales data incorporated with the Facebook prophet repository.

retail<-read.csv("https://raw.githubusercontent.com/facebookincubator/prophet/master/examples/example_retail_sales.csv",header=T,sep = ",")

We are forecasting the last 50 periods…

l=length(retail[,1])-50

prophet model:

Here is the code for the prophet model:

m <- prophet(retail[1:l,])
## Warning in set_auto_seasonalities(m): Disabling weekly seasonality. Run
## prophet with `weekly.seasonality=TRUE` to override this.
## Initial log joint probability = -2.42146
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
future <- make_future_dataframe(m,periods=50,freq = 'm')
retail_fcst <- predict(m, future)
plot(m,retail_fcst)

NNS model:

Here is the code for the NNS model:

nns.retail.fcst<- NNS.ARMA(retail[,2],h=50,seasonal.factor = 12,method = 'lin',training.set = l)

prophet MSE:

mean((tail(retail_fcst$yhat,50)-tail(retail[,2],50))^2)
## [1] 2392859440

NNS MSE:

mean((nns.retail.fcst-tail(retail[,2],50))^2)
## [1] 52793721

SCORE NNS: 1 prophet: 0

Peyton Manning Example (Raw Values):

Next is the raw Peyton Manning example from the Facebook prophet repository.

omaha<- read.csv("https://raw.githubusercontent.com/facebookincubator/prophet/master/examples/example_wp_peyton_manning.csv")

We are forecasting the last 365 periods…

l=length(omaha[,1])-365

prophet model:

m <- prophet(omaha[1:l,])
## Initial log joint probability = -4.0564
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
future <- make_future_dataframe(m,periods=365)
omaha_fcst <- predict(m, future)
plot(m,omaha_fcst)

NNS model:

NNS is using a seasonal.factor=c(1,sf) which translates to a daily and best reported lag from our seasonality test. This is not optimized.

a=NNS.seas(omaha[,2],plot=F)
test=t(sapply(a$all.periods$Period, function(i) 
  c(i,mean((NNS.ARMA(omaha[,2],h=365,training.set=l,
  seasonal.factor=i,plot=F,method='lin')-tail(omaha[,2],365))^2))))

colnames(test)=c("Period","MSE")

sf<- test[which.min(test[,2]),1]
sf
## Period 
##    331
nns.omaha.fcst<- NNS.ARMA(omaha[,2],h=365,seasonal.factor = c(1,sf),method = 'lin',training.set = l)

prophet MSE:

mean((tail(omaha_fcst$yhat,365)-tail(omaha[,2],365))^2)
## [1] 150018073

NNS MSE:

mean((nns.omaha.fcst-tail(omaha[,2],365))^2)
## [1] 101266229

SCORE NNS: 2 prophet: 0

Peyton Manning Example (log values):

Using the log values of the preceding Peyton Manning dataset.

log_omaha<- read.csv("https://raw.githubusercontent.com/facebookincubator/prophet/master/examples/example_wp_peyton_manning.csv") %>% mutate(y = log(y))

Again, we are forecasting the last 365 periods…

l=length(log_omaha[,1])-365

prophet model:

m <- prophet(log_omaha[1:l,])
## Initial log joint probability = -16.5944
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
future <- make_future_dataframe(m,periods=365)
log_omaha_fcst <- predict(m, future)
plot(m,log_omaha_fcst)

NNS model:

NNS is using a seasonal.factor=c(1,sf) which translates to a daily and best reported lag from our seasonality test. This is not optimized.

a=NNS.seas(log_omaha[,2],plot=F)
test=t(sapply(a$all.periods$Period, function(i) 
  c(i,mean((NNS.ARMA(log_omaha[,2],h=365,training.set=l,
  seasonal.factor=i,plot=F,method='lin')-tail(log_omaha[,2],365))^2))))

colnames(test)=c("Period","MSE")

sf<- test[which.min(test[,2]),1]
sf
## Period 
##    362
nns.log.omaha.fcst<- NNS.ARMA(log_omaha[,2],h=365,seasonal.factor = c(1,sf),method = 'nonlin',training.set = l)

prophet MSE:

mean((tail(log_omaha_fcst$yhat,365)-tail(log_omaha[,2],365))^2)
## [1] 0.3617284

NNS MSE:

mean((nns.log.omaha.fcst-tail(log_omaha[,2],365))^2)
## [1] 0.6019713

SCORE NNS: 2 prophet: 1

Notes

This dataset presents more questions than it answers. Specifically:

  • Why is one method not consistently better between raw and log datasets?

  • There is a stark difference in seasonal periods from NNS…does the log transformation really have that big of an effect on the peaks?

  • Can we recover better raw estimates from translated log estimates for prophet?

No. While better, it is still materially worse than NNS…

mean((exp(1)^tail(log_omaha_fcst$yhat,365)-tail(omaha[,2],365))^2)
## [1] 127292436
  • Can we recover better log estimates from translated raw estimates for NNS? No.
mean((log(nns.omaha.fcst)-tail(log_omaha[,2],365))^2)
## [1] 1.871977

lynx dataset

Next we test the Canadian lynx dataset suggested from Andrew Gelman’s blog post on prophet: http://andrewgelman.com/2017/03/01/facebooks-prophet-uses-stan/

data(lynx)
l=length(lynx)-34

prophet model and RMSE:

m <- prophet(data.frame(ds=1:80,y=lynx[1:80]))
## Warning in set_auto_seasonalities(m): Disabling yearly seasonality. Run
## prophet with `yearly.seasonality=TRUE` to override this.
## Initial log joint probability = -4.87324
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
future <- make_future_dataframe(m,periods=35)
fcast <- predict(m,future)
sqrt(mean((lynx[-(1:80)]-fcast[-(1:80),"yhat"])^2))
## Warning in lynx[-(1:80)] - fcast[-(1:80), "yhat"]: longer object length is
## not a multiple of shorter object length
## [1] 1823.442
plot(m,fcast)

NNS model and RMSE:

sqrt(mean((NNS.ARMA(lynx,h=34,training.set = 80,method='lin',seasonal.factor=F,best.periods=2)-tail(lynx,34))^2))

## [1] 1252.73

SCORE NNS: 3 prophet: 1

S&P 500

Finally we are going to forecast 252 S&P 500 adjusted closing prices. We retrieve the data via quantmod commands:

getSymbols("SPY",src='yahoo')
## [1] "SPY"
SPY=as.numeric(Ad(SPY))

prophet model:

l=length(SPY)-252
oos = tail(SPY,252)

m <- prophet(data.frame(ds=1:l,y=SPY[1:l]))
## Initial log joint probability = -33.9079
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
future <- make_future_dataframe(m,periods=252)
SPY_fcst <- predict(m, future)
plot(m,SPY_fcst)

prophet RMSE:

sqrt(mean((oos-tail(SPY_fcst$yhat,252))^2))
## [1] 30.20166

Linear Regression

We present the linear regression just for a visualization, as its residuals clearly do not deserve consideration.

plot(SPY)
abline(lm(SPY[1:l]~c(1:l)),col='red')

NNS model and RMSE:

NNS is using a seasonal.factor=c(1,5,21,63,126,252) which translates to a daily, weekly, monthly and quarterly, bi-annual and annual lag for S&P closing prices.

nns.sp = NNS.ARMA(SPY,h=252,seasonal.plot = F,seasonal.factor =c(1,5,21,63,126,252), training.set = l,method = 'nonlin')

sqrt(mean((oos-nns.sp)^2))
## [1] 22.7117

SCORE NNS: 4 prophet: 1

S&P500 - Year 2

So far so good…let’s try a different year:

prophet model and RMSE:

l=length(SPY)-252*2
oos = tail(SPY[1:l],252)

m <- prophet(data.frame(ds=1:l,y=SPY[1:l]))
## Initial log joint probability = -37.7735
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
future <- make_future_dataframe(m,periods=252)
SPY_fcst <- predict(m, future)
plot(m,SPY_fcst)

sqrt(mean((oos-tail(SPY_fcst$yhat,252))^2))
## [1] 18.51615

NNS model and RMSE:

nns.sp = NNS.ARMA(SPY,h=252,seasonal.plot = F,seasonal.factor =c(1,5,21,63,126,252), training.set = l,method = 'nonlin')

sqrt(mean((oos-nns.sp)^2))
## [1] 16.00169

SCORE NNS: 5 prophet: 1

S&P500 - Year 3

So far so good (again)…let’s try a different year:

prophet model and RMSE:

l=length(SPY)-252*3
oos = tail(SPY[1:l],252)

m <- prophet(data.frame(ds=1:l,y=SPY[1:l]))
## Initial log joint probability = -41.1845
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
future <- make_future_dataframe(m,periods=252)
SPY_fcst <- predict(m, future)
plot(m,SPY_fcst)

sqrt(mean((oos-tail(SPY_fcst$yhat,252))^2))
## [1] 27.84086

NNS model and RMSE:

nns.sp = NNS.ARMA(SPY,h=252,seasonal.plot = F,seasonal.factor =c(1,5,21,63,126,252), training.set = l,method = 'nonlin')

sqrt(mean((oos-nns.sp)^2))
## [1] 24.72559

SCORE NNS: 6 prophet: 1

Summary

To summarize, NNS was more accurate (often by substantial amounts) out of sample for all but the log transformed Peyton Manning dataset. However, prophet was exceptionally faster in its runtimes.

To be fair, prophet was not altered from its default settings and any interested readers are more than welcome to suggest different settings to improve performance. NNS could also improve performance as it was not optimized either.

Comments welcome at: ovvo.financial.systems@gmail.com

Further NNS.ARMA examples are available at: https://github.com/OVVO-Financial/NNS/blob/NNS-Beta-Version/examples/NNS%20ARMA.pdf

Thank you for your interest!