library(forecast)
tumblr<-[Link]([Link](), header=TRUE, sep=",")
tumblr_ts<- ts(tumblr$People,start=c(2010,4), frequency=12)
#plot various decompositions into error/noise, trend and seasonality
fit <- decompose(tumblr_ts, type="multiplicative") #decompose using "classical"
method, multiplicative form
plot(fit)
fit <- decompose(tumblr_ts, type="additive") #decompose using "classical" method,
additive form
plot(fit)
fit <- stl(tumblr_ts, [Link]=12, [Link]="periodic") #decompose using STL
(Season and trend using Loess)
plot(fit)
plot(tumblr_ts)
# Create exponential smoothing models: additive vs multiplicative noise (first A vs
M), additive vs multiplicative trend (second A vs M) and no seasonality vs
automatic detection (third N vs Z) trend and no seasonlity (AAN), multiplicative
(MMN)
tumblr_AAN <- ets(tumblr_ts, model="AAN")
tumblr_AAZ <- ets(tumblr_ts, model="AAZ", damped=TRUE)
tumblr_MMN <- ets(tumblr_ts, model="MMN", damped=TRUE)
tumblr_MMZ <- ets(tumblr_ts, model="MMZ", damped=TRUE)
# Create their prediction "cones" for 115 months into the future with quintile
confidence intervals
tumblr_AAN_pred <- forecast(tumblr_AAN, h=115, level=c(0.8, 0.95))#80% price in
between 95% confidence
tumblr_AAZ_pred <- forecast(tumblr_AAZ, h=115, level=c(0.8, 0.95))
tumblr_MMN_pred <- forecast(tumblr_MMN, h=115, level=c(0.8, 0.95))
tumblr_MMZ_pred <- forecast(tumblr_MMZ, h=115, level=c(0.8, 0.95))
# Compare the prediction "cones" visually
par(mfrow=c(1,4)) # This command sets the plot window to show 1 row of 4 plots
plot(tumblr_AAN_pred, xlab="Year", ylab="Predicted People")
plot(tumblr_MMN_pred, xlab="Year", ylab="Predicted People")
plot(tumblr_AAZ_pred, xlab="Year", ylab="Predicted People")
plot(tumblr_MMZ_pred, xlab="Year", ylab="Predicted People")
# Lets look at what our models actually are -- ETS
tumblr_AAZ
tumblr_MMZ
#Create a trigonometric box-cox autoregressive trend seasonality (TBATS) model
tumblr_tbats <- tbats(tumblr_ts)
tumblr_tbats_pred <-forecast(tumblr_tbats, h=115, level=c(0.8, 0.95))
par(mfrow=c(1,1))
plot(tumblr_tbats_pred, xlab="Year", ylab="Predicted People")
par(mfrow=c(1,3)) # Lets look at the three models with seasonality on one graph on
the same scale
plot(tumblr_AAZ_pred, xlab="Year", ylab="Predicted People",
ylim=c(100000000,200000000))
plot(tumblr_MMZ_pred, xlab="Year", ylab="Predicted People",
ylim=c(100000000,200000000))
plot(tumblr_tbats_pred, xlab="Year", ylab="Predicted People",
ylim=c(100000000,200000000))
# Lets look at what our models actually are -- TBATS
tumblr_tbats
###
### Comparing models -- Time series Cross Validation (Rolling Horizon Holdout)
###
f_AAN <- function(y, h) forecast(ets(y, model="AAN"), h = h)
errors_AAN <- tsCV(tumblr_ts, f_AAN, h=1)
f_MMN <- function(y, h) forecast(ets(y, model="MMN"), h = h)
errors_MMN <- tsCV(tumblr_ts, f_MMN, h=1)
f_AAA <- function(y, h) forecast(ets(y, model="AAA"), h = h)
errors_AAA <- tsCV(tumblr_ts, f_AAA, h=1)
f_MMM <- function(y, h) forecast(ets(y, model="MMM"), h = h)
errors_MMM <- tsCV(tumblr_ts, f_MMM, h=1)
par(mfrow=c(1,1))
plot(errors_AAN, ylab='tsCV errors')
abline(0,0)
lines(errors_MMN, col="red")
lines(errors_AAA, col="green")
lines(errors_MMM, col="blue")
legend("left", legend=c("CV_error_AAN",
"CV_error_MMN","CV_error_AAA","CV_error_MMM"), col=c("black", "red", "green",
"blue"), lty=1:4)
mean(abs(errors_AAN/tumblr_ts), [Link]=TRUE)*100
mean(abs(errors_MMN/tumblr_ts), [Link]=TRUE)*100
mean(abs(errors_AAA/tumblr_ts), [Link]=TRUE)*100
mean(abs(errors_MMM/tumblr_ts), [Link]=TRUE)*100
f_TBATS <- function(y, h) forecast(tbats(y), h = h)
errors_TBATS <- tsCV(tumblr_ts, f_TBATS, h=1, window=37)
plot(errors_AAA, ylab='tsCV errors', col="green")
abline(0,0)
lines(errors_MMM, col="blue")
lines(errors_TBATS, col="gray")
legend("left", legend=c("CV_error_AAA", "CV_error_MMM","CV_error_TBATS"),
col=c("green", "blue", "gray"), lty=1:4)
mean(abs(errors_TBATS/tumblr_ts), [Link]=TRUE)*100
#ARIMA
par(mfrow=c(1,3))
plot(tumblr_ts, xlab="Year",
ylab="People")
plot(log(tumblr_ts), xlab="Year",
ylab="log People")
plot(diff(log(tumblr_ts),12), xlab="Year",
ylab="Annual change in monthly log People")
# decomposition
fit <- stl(log(tumblr_ts), [Link]=12, [Link]="periodic", robust=TRUE)
plot(fit)
# fit the heck
par(mfrow=c(1,2))
Acf(diff(log(tumblr_ts),12)) # auto-correlation function
Pacf(diff(log(tumblr_ts),12)) # partial auto-correlation function
# fit ARIMA models
# non-seasonal first
fit <- [Link](tumblr_ts,seasonal=FALSE)
fit
arima_nons <- fit # create an object for non seasonal arima
plot(forecast(arima_nons,115))
#check residuals for autocorrelation
par(mfrow=c(1,1))
Acf(residuals(fit))
# and now seasonal
fit <- [Link](tumblr_ts,seasonal=TRUE)
fit
arima_seas <- fit # create an object for seasonal arima
arima_ffff <-forecast(arima_seas,115)
plot(forecast(arima_seas,115))
#errors
mean(abs(residuals(fit)/tumblr_ts), [Link]=TRUE)*100
# Print the mean and confidence intervals for the MMZ model
tumblr_AAN_pred
# Export the results out
[Link](tumblr_AAN_pred, file = "Predicted [Link]") # export the selected
model's predictions into a CSV file