library(fpp);
###################################
#[Link]()
plot(usconsumption[, 1], xlab = "Year", ylab = "Quarterly percentage change", main = "US consumption")
[Link](usconsumption[,1], seasonal=FALSE)
fit <- [Link](usconsumption[,1], seasonal=FALSE)
plot(forecast(fit, h=5))
#################################
#acf and pacf
par(mfrow=c(1,2))
Acf(usconsumption[,1],main="")
Pacf(usconsumption[,1],main="")
#################################
#applying the procedure
#plot the data
eeadj <- seasadj(stl(elecequip, [Link]="periodic"))
plot(eeadj, ylab = "Seasonally adjusted new orders index", xlab = "Year")
#no need for Box-Cox transformation to stabilize the data
#data are not stationary, so we try differencing
tsdisplay(diff(eeadj), main="")
#use unit-root test to make sure that the data are stationary
[Link](diff(eeadj),alternative = "stationary")
[Link](diff(eeadj))
#find the appropriate number of fist differences required
#ATTENSION, the following code only applies to non-seasonal time series
#nd will hold the degree of first differencing
#xstar will hold the appropriately differenced time series
ns <- nsdiffs(eeadj)
if (ns>0) {
xstar <- diff(eeadj, lag=frequency(eeadj),differences=ns)
}else{
xstar <- eeadj
}
nd <- ndiffs(xstar)
if(nd > 0){
xstar <- diff(xstar,differences = nd)
}
#after understanding that fist difference of the data is all we need
#we analyze again the ACF and PACF data
tsdisplay(diff(eeadj), main="")
#the ARIMA(3,1,0) is indicated from the PACF
#trying different variations
Arima(eeadj, c(3, 1, 0))
Arima(eeadj, c(3, 1, 1))
Arima(eeadj, c(4, 1, 0))
Arima(eeadj, c(2, 1, 0))
#the ARIMA(3,1,1) has slightly smaller AIC_c
fit <- Arima(eeadj, order=c(3,1,1))
summary(fit)
#checking the residuals for autocorrelation and doing a portmanteau test
Acf(residuals(fit))
[Link](residuals(fit), lag=24, fitdf=4, type="Ljung")
#ploting the forecast
plot(forecast(fit))
#using the [Link]()
[Link](seasadj(stl(elecequip, [Link]="periodic")))
###################################################
# Applying the procedure to a seasonal time series
#plotting the data for the quarterly European retail trade from 1996 to 2011
plot(euretail, ylab="Retail index", xlab="Year")
#step1
tsdisplay(diff(euretail,4))
#step2
tsdisplay(diff(diff(euretail,4)))
#step3
fit <- Arima(euretail, order=c(0,1,1), seasonal=c(0,1,1))
#step4
tsdisplay(residuals(fit))
#step5
fit2 <- Arima(euretail, order = c(0, 1, 2), seasonal = c(0, 1, 1))
fit3 <- Arima(euretail, order = c(0, 1, 3), seasonal = c(0, 1, 1))
fit4 <- Arima(euretail, order = c(1, 1, 1), seasonal = c(0, 1, 1))
#step6
# fit3 the best with the smallest AIC_c
res <- residuals(fit3)
tsdisplay(res, main = "")
[Link](res, lag=16, fitdf=4, type="Ljung")
#step7
plot(forecast(fit3, h=12))
#step8
[Link](euretail)
#step9
[Link](euretail, stepwise=FALSE, approximation=FALSE) #turnign the short-cuts off
##########################################################
#take a look at it too
lh02 <- log(h02)
par(mfrow=c(2,1))
plot(h02, ylab="H02 sales (million scripts)", xlab="Year")
plot(lh02, ylab="Log H02 sales", xlab="Year")
tsdisplay(diff(lh02,12),
main="Seasonally differenced H02 scripts", xlab="Year")
[Link](h02, lambda = 0)
Arima(h02, order = c(3, 0, 0), seasonal = c(2, 1, 0), lambda = 0)
Arima(h02, order = c(3, 0, 1), seasonal = c(2, 1, 0), lambda = 0)
Arima(h02, order = c(3, 0, 2), seasonal = c(2, 1, 0), lambda = 0)
Arima(h02, order = c(3, 0, 1), seasonal = c(1, 1, 0), lambda = 0)
Arima(h02, order = c(3, 0, 1), seasonal = c(0, 1, 1), lambda = 0)
Arima(h02, order = c(3, 0, 1), seasonal = c(0, 1, 2), lambda = 0)
Arima(h02, order = c(3, 0, 1), seasonal = c(1, 1, 1), lambda = 0)
fit <- Arima(h02, order=c(3,0,1), seasonal=c(0,1,2), lambda=0)
tsdisplay(residuals(fit))
[Link](residuals(fit), lag=36, fitdf=6, type="Ljung")
fit <- [Link](h02, lambda=0, d=0, D=1, [Link]=9,
stepwise=FALSE, approximation=FALSE)
tsdisplay(residuals(fit))
[Link](residuals(fit), lag=36, fitdf=8, type="Ljung")
getrmse <- function(x,h,...)
{
[Link] <- time(x)[length(x)-h]
[Link] <- time(x)[length(x)-h+1]
train <- window(x,end=[Link])
test <- window(x,start=[Link])
fit <- Arima(train,...)
fc <- forecast(fit,h=h)
return(accuracy(fc,test)[2,"RMSE"])
}
getrmse(h02,h=24,order=c(3,0,0),seasonal=c(2,1,0),lambda=0)
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(2,1,0),lambda=0)
getrmse(h02,h=24,order=c(3,0,2),seasonal=c(2,1,0),lambda=0)
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(1,1,0),lambda=0)
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(0,1,2),lambda=0) #best
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(1,1,1),lambda=0)
getrmse(h02,h=24,order=c(4,0,3),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(3,0,3),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(4,0,2),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(3,0,2),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(2,1,3),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(2,1,4),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(2,1,5),seasonal=c(0,1,1),lambda=0)
fit <- Arima(h02, order=c(3,0,1), seasonal=c(0,1,2), lambda=0)
plot(forecast(fit), ylab="H02 sales (million scripts)", xlab="Year")