Arima
1 Arima
ARIMA book: https://www.otexts.org/fpp
- Generally, seasonal versions of ARIMA and ETS models are designed for shorter periods such as 12 for monthly data or 4 for quarterly data.
comparing results of applying transformation and seasonal configuration in arima
- forecast
fit <- Arima(WWWusage, order=c(3,1,0))
# plot the model against the data
plot(fit$x,col="black")
lines(fitted(fit),col="red")
# auto arima
require(forecast)
fit <- auto.arima(ts)
plot(forecast(fit, h=5)) # 5 stps ahead
# Model Diagnostics for a Fitted ARIMA Model
# gof.lag = maximum lag used in ACF and Ljung-Box tests for the residuals
tsdiag(fit,gof=50,omit.initial=FALSE)
# the 5th value in 5 stps ahead
forecast(fit, h=5)$mean[5]
# extract the (p, d, q) ; more info for arima.string() function
order <- c(fit$arma[1], fit$arma[6], fit$arma[2])
# or
order <- arimaorder(fit)2 regressor/covariate
hour <- as.POSIXlt(time(xts1))$hour # extract the hour of time stamp
xreg <- model.matrix(~as.factor(hour)) # hourly
week <- as.POSIXlt(time(xts1))$wday # find the day of the week
xreg <- model.matrix(~as.factor(week)) # week and weekend days
# Remove intercept
xreg <- xreg[,-1]
# Rename columns
# colnames(xreg) <- c("Mon","Tue","Wed","Thu","Fri","Sat")
fit2 <- auto.arima(... , xreg=xreg)
tdata <- s
tdata <- ts(tdata)
fit1 <- auto.arima(tdata)
fit1
#plot arima fitted model with the original series
plot(fit$x) # original
lines(fitted(fit),col="blue")
system.time(fit <- Arima(tdata, model=fit1 ))
fcast <- forecast(fit, h=5)
plot(fcast,include=50)
lines(tdata, col="red")
accuracy(fcast,tdata)
order <- c(fit1$arma[1], fit1$arma[6], fit1$arma[2])
system.time(fit <- arima(tdata, order=order) )
fcast <- forecast(fit, h=5)
plot(fcast,include=50)
lines(tdata, col="red")
accuracy(fcast,tdata)3 Arima config
Comparing different configurations with Arima forecasting
library(fpp) # for a10 dataset## Loading required package: forecast## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: timeDate## This is forecast 7.3## Loading required package: fma## Loading required package: expsmooth## Loading required package: lmtest## Loading required package: tseriesplot_fit <- function(fit){
fcast <- forecast(fit, h=50)
plot(fcast)
lines(a10, col="black")
accuracy(fcast, a10)
}
a10train <- window(a10, end=2004.99)
#
# par(mfrow=c(5,1),
# mar=c(2,2,2,2), oma=c(0,0,0,0)
# )
plot(a10train)
# ------ auto.arima only -------------
fit <- auto.arima(a10train)
fit1 <- fit
plot_fit(fit)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05685677 0.502642 0.3617925 0.2155619 4.196051 0.3730364
## Test set -0.91056683 1.998748 1.5569010 -6.4951797 9.028741 1.6052868
## ACF1 Theil's U
## Training set -0.02572485 NA
## Test set 0.14316285 0.6173202# ------ auto.arima + transformation ------------------
# check
if(Box.test(a10train)$p.value < 0.05){
print("transformation is unnecessary")
}## [1] "transformation is unnecessary"# force transformation
lam <- BoxCox.lambda(a10train) # BoxCox transformation
fit2 <- auto.arima(a10train, lambda=lam)
plot_fit(fit2)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01335119 0.4564804 0.3284141 -0.07159435 3.796341 0.3386206
## Test set 1.02892014 2.0869887 1.6719146 3.82741618 8.089647 1.7238748
## ACF1 Theil's U
## Training set -0.06675867 NA
## Test set 0.30267576 0.5640996# ------ Arima + transformation + order + seasonal ---------------------
fit <- auto.arima(a10train)
# fit <- auto.arima(a10train, lambda=lam) gives error
order <- c(fit$arma[1], fit$arma[6], fit$arma[2])
seasonal <- list(order=c(fit$arma[3], fit$arma[7], fit$arma[4]), period=fit$arma[5])
fit3 <- Arima(a10train, order=order, seasonal=seasonal, lambda=lam)
plot_fit(fit3)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02975795 0.4373383 0.3096405 0.1333841 3.540816 0.3192636
## Test set -0.37381024 1.7016682 1.4207328 -3.4298686 7.848789 1.4648867
## ACF1 Theil's U
## Training set -0.0559194 NA
## Test set 0.1187889 0.5142006# ------ Arima + order ---------------------
fit <- auto.arima(a10train)
fit4 <- Arima(a10train, order=arimaorder(fit)[1:3])
plot_fit(fit4)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3049229 1.403439 0.9517912 1.430843 11.45887 0.9813712
## Test set 3.1729505 5.183330 4.0732277 12.397250 19.11872 4.1998165
## ACF1 Theil's U
## Training set -0.04940304 NA
## Test set 0.57701894 1.345497# ------ Arima + order + transformation ---------------------
fit <- auto.arima(a10train, lambda=lam)
fit5 <- Arima(a10train, order=arimaorder(fit)[1:3], lambda=lam)
plot_fit(fit5)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3554298 1.410858 0.9846879 1.877586 11.76367 1.015290
## Test set 3.6324311 5.564654 4.4160264 14.727246 20.61579 4.553269
## ACF1 Theil's U
## Training set -0.06488317 NA
## Test set 0.60180135 1.442611# ------ Arima + order + xreg ---------------------
# build xreg
months <- array( seq(1:12), 162+50+6)
months <- months[7:length(months)]
xreg <- model.matrix(~as.factor(months))
xreg <- xreg[, -1]
# get order
fit <- auto.arima(a10train, xreg=xreg[1:162,])
fit6<- Arima(a10train, order=arimaorder(fit)[1:3], xreg=xreg[1:162,])
fcast6 <- forecast(fit6, h=50, xreg=xreg[163:(162+50),])
plot(fcast6)
lines(a10, col="black")
accuracy(fcast6, a10)## ME RMSE MAE MPE MAPE MASE
## Training set 0.171548 0.6565298 0.5026238 2.130656 6.697535 0.5182445
## Test set 2.742901 4.3559570 3.4210346 11.012524 16.019891 3.5273544
## ACF1 Theil's U
## Training set -0.08619346 NA
## Test set 0.67299382 1.141562# ------- Arima + drift ---------------------
fit <- auto.arima(a10train)
fit7 <- Arima(a10train, order=arimaorder(fit)[1:3], include.drift=T)
fcast7 <- forecast(fit7, h=50)
plot_fit(fit7)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.02179321 1.333284 0.9158729 -3.642429 11.65504 0.9443367
## Test set 2.39311903 4.287831 3.3245427 8.930366 15.71104 3.4278637
## ACF1 Theil's U
## Training set 0.002848839 NA
## Test set 0.435836493 1.107696# ------------------------------