Hi all! In this section, we will explore some basic times series.
Special thanks to NUS ISS for teaching the codes.
Download TimesSeriesMen.csv here
Times Series - ts()
# read data and store to tsm
tsm <- read.csv("TimeSeriesMen.csv", stringsAsFactors = FALSE)
str(tsm)
## 'data.frame': 120 obs. of 12 variables:
## $ date : chr "1/1/1989" "2/1/1989" "3/1/1989" "4/1/1989" ...
## $ men : num 11358 10606 16999 6564 6608 ...
## $ women : num 16579 18236 43394 30908 28702 ...
## $ jewel : num 10776 10822 22846 11103 16067 ...
## $ mail : int 7978 8290 8029 7752 8685 7847 7881 8121 7811 8706 ...
## $ page : int 73 88 65 85 74 87 79 72 83 111 ...
## $ phone : int 34 29 24 20 17 30 28 27 35 25 ...
## $ print : num 22294 27426 27979 28950 22642 ...
## $ service: int 20 20 26 22 21 23 22 20 15 20 ...
## $ YEAR_ : int 1989 1989 1989 1989 1989 1989 1989 1989 1989 1989 ...
## $ MONTH_ : int 1 2 3 4 5 6 7 8 9 10 ...
## $ DATE_ : chr "Jan-89" "Feb-89" "Mar-89" "Apr-89" ...
# create a time series call men, freq of 12 months.
men <- ts(tsm$men, frequency = 12, start = c(1989,1))
# plot
plot(men, lwd = 2, col = "blue", main = "Men Series")
# similarly
women <- ts(tsm$women, frequency = 12, start = c(1989,1))
jewel <- ts(tsm$jewel, frequency = 12, start = c(1989,1))
par(mfrow = c(1,2))
plot(women, lwd = 2, col = "green3", main = "Women Series")
plot(jewel, lwd = 2, col = "red", main = "Jewel Series")
par(mfrow = c(1,1))
# Comparing
plot(men, lwd = 2, col = "blue")
lines(jewel, lwd = 2, col = "red")
plot(men, lwd = 2, col = "blue")
lines(0.75*(women-15000), lwd =2, col = "green")
Holt-Winters - HoltWinters()
# Holt Winters Addditive
fitHWA <- HoltWinters(men, seasonal = "additive")
fitHWA
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = men, seasonal = "additive")
##
## Smoothing parameters:
## alpha: 0.004152723
## beta : 0.4349592
## gamma: 0
##
## Coefficients:
## [,1]
## a 22287.9136
## b 132.5304
## s1 423.6983
## s2 -3722.2834
## s3 -2116.2092
## s4 -798.1288
## s5 -2058.7850
## s6 -2127.7738
## s7 -2523.7142
## s8 -1469.6505
## s9 184.1891
## s10 1409.7529
## s11 1878.7212
## s12 10920.1833
plot(fitHWA, lwd = 2, main = "Holt Winters Additive")
# Auto Cross Variance and Correlation Function - acf
acf(residuals(fitHWA), 24, lwd = 2, col = "blue")
# Box-Ljung test
Box.test(residuals(fitHWA), lag = 12, type = "Ljung")
##
## Box-Ljung test
##
## data: residuals(fitHWA)
## X-squared = 14.738, df = 12, p-value = 0.2561
# forecast
library(forecast)
f <- forecast(fitHWA, h = 12, level = c(80,95))
plot(f, lwd = 1, col = "green3")
# accuracy
accuracy(f)
## ME RMSE MAE MPE MAPE MASE
## Training set 325.7033 3534.413 2313.992 -5.742507 18.67779 0.642404
## ACF1
## Training set 0.05592869
# Holt Winters Multiplicative
fitHWM <- HoltWinters(men, seasonal = "multiplicative")
fitHWM
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = men, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha: 0.00361763
## beta : 0.2455624
## gamma: 0.07278723
##
## Coefficients:
## [,1]
## a 2.184028e+04
## b 1.190637e+02
## s1 1.017983e+00
## s2 7.876398e-01
## s3 8.425741e-01
## s4 9.055215e-01
## s5 8.703500e-01
## s6 8.560653e-01
## s7 8.317946e-01
## s8 9.497353e-01
## s9 9.717868e-01
## s10 1.133625e+00
## s11 1.129739e+00
## s12 1.849141e+00
plot(fitHWM, lwd = 2, main = "Holt Winters Multiplicative")
# Auto Cross Variance and Correlation Function - acf
acf(residuals(fitHWM), 24, lwd = 2, col = "blue")
# Box-Ljung test
Box.test(residuals(fitHWM), lag = 12, type = "Ljung")
##
## Box-Ljung test
##
## data: residuals(fitHWM)
## X-squared = 17.543, df = 12, p-value = 0.1303
# forecast
library(forecast)
f <- forecast(fitHWM, h = 12, level = c(80,95))
plot(f, lwd = 1, col = "green3")
# accuracy
accuracy(f)
## ME RMSE MAE MPE MAPE MASE
## Training set 298.1092 3537.055 2308.113 -4.958027 18.78873 0.6407718
## ACF1
## Training set 0.08713956
Error Trend Seasonality - ets()
# ETS
fitE <- ets(men)
fitE
## ETS(M,A,M)
##
## Call:
## ets(y = men)
##
## Smoothing parameters:
## alpha = 0.0394
## beta = 0.0037
## gamma = 1e-04
##
## Initial states:
## l = 11599.5462
## b = 35.3159
## s=1.6666 1.0477 1.082 1.0212 0.9791 0.852
## 0.918 0.8819 0.8468 0.8998 0.8671 0.9378
##
## sigma: 0.2158
##
## AIC AICc BIC
## 2553.219 2559.219 2600.607
plot(fitE, lwd = 1, col = "red")
acf(residuals(fitE), 24, lwd = 2, col = "blue")
Box.test(residuals(fitE), lag = 12, type = "Ljung")
##
## Box-Ljung test
##
## data: residuals(fitE)
## X-squared = 9.6485, df = 12, p-value = 0.6468
Arima with Regressors
# splitting data set into train & test
train <- subset(tsm, YEAR_ != 1998)
test <- subset(tsm, YEAR_ == 1998)
xregtrain <- train[c(5,6,7,8,9)]
xregtest <- test[c(5,6,7,8,9)]
men_train <- ts(train$men, frequency = 12, start = c(1989,1))
men_test <- ts(test$men, frequency = 12, start = c(1998,1))
TrainArimaReg <- auto.arima(men_train, xreg = xregtrain)
TrainArimaReg
## Series: men_train
## Regression with ARIMA(0,0,1)(1,0,0)[12] errors
##
## Coefficients:
## ma1 sar1 intercept mail page phone print
## 0.1309 0.3351 -20945.860 1.8392 19.7837 373.1267 0.1727
## s.e. 0.0935 0.1003 2995.366 0.2086 19.2173 41.0745 0.0713
## service
## -23.8064
## s.e. 39.5992
##
## sigma^2 estimated as 6826800: log likelihood=-999.58
## AIC=2017.15 AICc=2018.99 BIC=2041.29
fTestArimaReg <- forecast.Arima(TrainArimaReg, xreg = xregtrain)
Box.test(residuals(fTestArimaReg), lag = 12, type = "Ljung")
##
## Box-Ljung test
##
## data: residuals(fTestArimaReg)
## X-squared = 8.1743, df = 12, p-value = 0.7714
accuracy(fTestArimaReg, men)
## ME RMSE MAE MPE MAPE MASE
## Training set -29.0428 2514.182 1795.227 -2.46764 12.33812 0.5130855
## Test set 11738.4411 12171.368 11738.441 52.81488 52.81488 3.3549099
## ACF1 Theil's U
## Training set 0.01336445 NA
## Test set -0.25024289 2.043187
plot(fTestArimaReg, lwd = 2, col = "green3")
lines(men_test, lwd = 2, col = "red")
