PROPHET model

We use Prophet models to:

Prophet can potentially produce superior results to more traditional time series models such as ARIMA by identifying structural breaks in a time series and making forecasts by taking change points as well as seasonality patterns into account.

Model Building

library(prophet)
library(lubridate)

prophet_data = annual_temp
prophet_data$ds <-  as.character(prophet_data$ds)
prophet_data$ds = paste(prophet_data$ds,'-01-01',sep = "")


train = prophet_data %>% 
    filter(ds<ymd("2012-01-01"))
test = prophet_data %>%
    filter(ds>=ymd("2012-01-01"))


model = prophet(train, weekly.seasonality=FALSE, daily.seasonality=FALSE)

The model is then used to predict 5 years forward, with the predictions compared to the test set (actual values).

future = make_future_dataframe(model,periods = 5, freq="year")
forecast = predict(model,future)

plot(model, forecast) +
ylab("Temp Anomalies")+xlab("Year")+theme_bw()

Visualizing the Components

Very similar to general time series decomposition

prophet_plot_components(model,forecast)

Plotting Changepoints

Most real time series frequently have abrupt changes in their trajectories. These are referred to as changepoints. Identifying them could be important for accurately modelling changes in trends. By default, prophet identifies the changepoints. However, we can also control this using arguments.

plot(model,forecast)+add_changepoints_to_plot(model)+theme_bw()+xlab("Year")+ylab("Temp Anomalies")

model = prophet(train, n.changepoints=50, weekly.seasonality=FALSE, daily.seasonality=FALSE)
forecast = predict(model,future)
plot(model,forecast)+add_changepoints_to_plot(model)+theme_bw()+xlab("Year")+ylab("Temp Anomalies")

Saturation Point

The model doesn’t need to consider the minimum/maximum point.

Seasonality

Prophet can handle multiple types of seasonality simultaneously

For example – day of week AND yearly seasonality

Yearly seasonality estimated as smooths, we can control “wiggliness”

Since we only have yearly datawe cant expect much seasonality

prophet_plot_components(model,forecast)

mul = prophet(train, ,seasonality.mode = 'multiplicative', weekly.seasonality=FALSE, daily.seasonality=FALSE)
mul_fcst = predict(mul,future)
plot(mul,mul_fcst)

prophet_plot_components(mul, mul_fcst)

Estimation of in-sample performance of the model based on RMSE and associated visualizations of model performance

forecast_metric_data = forecast %>% 
  as_tibble() %>% 
  mutate(ds = as.Date(ds)) %>% 
  filter(ds>=ymd("2012-01-01"))
RMSE = sqrt(mean((test$y - forecast_metric_data$yhat)^2))
MAE = mean(abs(test$y - forecast_metric_data$yhat))
MAPE = mean(abs((test$y - forecast_metric_data$yhat)/test$y))
print(paste("RMSE:",round(RMSE,2)))
## [1] "RMSE: 0.26"
print(paste("MAE:",round(MAE,2)))
## [1] "MAE: 0.22"
print(paste("MAPE:",round(MAPE,2)))
## [1] "MAPE: 0.26"
df.cv <- cross_validation(model, initial = 730, period = 30, horizon = 30, units = 'days')
head(df.cv)
## # A tibble: 6 × 6
##       y ds                     yhat yhat_lower yhat_upper cutoff             
##   <dbl> <dttm>                <dbl>      <dbl>      <dbl> <dttm>             
## 1 -0.21 1883-01-01 00:00:00 -0.0804    -0.0804    -0.0804 1882-12-02 00:00:00
## 2 -0.28 1884-01-01 00:00:00 -0.234     -0.266     -0.200  1883-12-02 00:00:00
## 3 -0.32 1885-01-01 00:00:00 -0.157     -0.191     -0.128  1884-12-02 00:00:00
## 4 -0.31 1886-01-01 00:00:00 -0.323     -0.380     -0.265  1885-12-02 00:00:00
## 5 -0.33 1887-01-01 00:00:00 -0.374     -0.428     -0.321  1886-12-02 00:00:00
## 6 -0.2  1888-01-01 00:00:00 -0.407     -0.455     -0.356  1887-12-02 00:00:00