title: “hw1forecasting” author: “Rufeng Liu” date: “7/22/2018” output: html_document
###6.9 Exercises #7
#Check if we need transformation
library(fma)## Loading required package: forecastplot(writing)
hist(writing)
#Using stlf
fcast.wrt.naive <- stlf(writing, method='naive')
fcast.wrt.drift <- stlf(writing, method='rwdrift')
res.naive<-as.numeric(residuals(fcast.wrt.naive))
res.drift<-as.numeric(residuals(fcast.wrt.drift))
res.naive<-na.omit(res.naive)
res.drift<-na.omit(res.drift)
#Compare RMSE
RMSE.naive<-sqrt(sum((res.naive)^2)/120)
RMSE.drift<-sqrt(sum((res.drift)^2)/120)
RMSE.naive## [1] 41.47904RMSE.drift## [1] 41.34016#Based on the results of the RMSE, the numbers are very close.
###6.9 Exercises #8
#Check if we need transformation
fancyts=ts(fancy,frequency=12)
plot(fancyts)
hist(fancyts)
#Using stlf
forecast.fancy.naive <- stlf(writing, method='naive')
forecast.fancy.drift <- stlf(writing, method='rwdrift')
#the histogram suggests we need a log transformation
forecast.fancy.log<- stlf(writing, method='rwdrift', lambda=0)
resnaive<-as.numeric(residuals(forecast.fancy.naive))
resdrift<-as.numeric(residuals(forecast.fancy.drift))
reslog<-as.numeric(residuals(forecast.fancy.log, type="response"))
resnaive<-na.omit(resnaive)
resdrift<-na.omit(resdrift)
reslog<-na.omit(reslog)
#Compare RMSE
RMSEnaive<-sqrt(sum((resnaive)^2)/119)
RMSEdrift<-sqrt(sum((resdrift)^2)/119)
RMSElog<-sqrt(sum((reslog)^2)/119)
RMSEnaive## [1] 41.65295RMSEdrift## [1] 41.5135RMSElog## [1] 40.89991#Model log has the the lowest RMSE.
###7.8 Exercise #13
#Compare ets(), snaive() and stlf() on the following six time series
library(fpp)## Loading required package: expsmooth## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: tseriesplot(ausbeer)
plot(dole)
plot(usmelec)
#t=train set,te=test set
aus.t<-window(ausbeer, end=c(2006,4))
aus.te<-window(ausbeer,start=c(2007,1))
dole.t<-window(dole, end=c(1988,12))
dole.te<-window(dole,start=c(1989,1))
usm.t<-window(usmelec,end=c(2009,12))
usm.te<-window(usmelec,start=c(2010,1))
#AUSBEER
aus.ets<-ets(aus.t, model="ZZZ")
fcast.aus.ets<-forecast(aus.ets,h=14)
aus.naive<-snaive(aus.t,h=14)
aus.stlf<-stlf(aus.t)
fcast.aus.stlf<-forecast(aus.stlf)
accuracy(fcast.aus.ets,aus.te)## ME RMSE MAE MPE MAPE MASE
## Training set -0.3389651 15.854046 12.074192 -0.06408683 2.886788 0.7584292
## Test set -4.7986612 7.869055 6.450648 -1.13019666 1.538190 0.4051914
## ACF1 Theil's U
## Training set -0.1956448 NA
## Test set 0.4708623 0.1938648accuracy(aus.naive,aus.te)## ME RMSE MAE MPE MAPE MASE
## Training set 3.440000 19.75171 15.92 0.9302419 3.792287 1.0000000
## Test set -7.428571 11.58817 10.00 -1.7139055 2.355378 0.6281407
## ACF1 Theil's U
## Training set 0.0007894127 NA
## Test set 0.2512162760 0.2546494accuracy(fcast.aus.stlf,aus.te)## ME RMSE MAE MPE MAPE MASE
## Training set -0.04281057 12.628969 9.858354 -0.01528831 2.364963 0.6192433
## Test set -8.07198572 9.528903 8.071986 -1.96356230 1.963562 0.5070343
## ACF1 Theil's U
## Training set -0.1473240 NA
## Test set 0.2183406 0.2281952#DOLE
dole.ets<-ets(dole.t, model="ZZZ")
fcast.dole.ets<-forecast(dole.ets,h=45)
dole.naive<-snaive(dole.t,h=45)
dole.stlf<-stlf(dole.t)
fcast.dole.stlf<-forecast(dole.stlf)
accuracy(fcast.dole.ets,dole.te)## ME RMSE MAE MPE MAPE
## Training set 114.4834 17359.03 9995.905 0.5432515 6.388437
## Test set 132352.4062 219319.58 163712.826 16.7937055 24.761948
## MASE ACF1 Theil's U
## Training set 0.3272434 0.4760381 NA
## Test set 5.3595897 0.9450812 8.872678accuracy(dole.naive,dole.te)## ME RMSE MAE MPE MAPE MASE
## Training set 14337.50 55786.85 30545.78 3.782368 27.87116 1.000000
## Test set 80642.42 188752.67 150555.30 6.793096 24.54455 4.928842
## ACF1 Theil's U
## Training set 0.9780301 NA
## Test set 0.9416940 8.055767accuracy(fcast.dole.stlf,dole.te)## ME RMSE MAE MPE MAPE MASE
## Training set 200.2779 5423.679 3243.546 0.4990998 4.072099 0.1061864
## Test set 28235.1679 64982.702 42682.864 5.9036207 9.470565 1.3973410
## ACF1 Theil's U
## Training set 0.005627286 NA
## Test set 0.815065718 4.039644#USMELEC
usm.ets<-ets(usm.t, model="ZZZ")
fcast.usm.ets<-forecast(usm.ets,h=42)
usm.naive<-snaive(usm.t,h=42)
usm.stlf<-stlf(usm.t)
fcast.usm.stlf<-forecast(usm.stlf)
accuracy(fcast.usm.ets,usm.te)## ME RMSE MAE MPE MAPE MASE
## Training set 0.6379805 7.417697 5.531263 0.1950478 2.149658 0.6176162
## Test set 10.0168896 14.238482 12.485876 2.6440021 3.474077 1.3941624
## ACF1 Theil's U
## Training set 0.2024459 NA
## Test set 0.3581895 0.3560936accuracy(usm.naive,usm.te)## ME RMSE MAE MPE MAPE MASE
## Training set 4.827868 11.52154 8.955826 1.981702 3.511752 1.000000
## Test set 15.079600 19.86604 15.521000 4.057444 4.211093 1.733062
## ACF1 Theil's U
## Training set 0.4899769 NA
## Test set 0.4406002 0.5318733accuracy(fcast.usm.stlf,usm.te)## ME RMSE MAE MPE MAPE
## Training set -0.06374684 5.557921 4.228934 -0.05122201 1.638225
## Test set 6.84787629 12.009666 10.978768 1.73808338 3.111847
## MASE ACF1 Theil's U
## Training set 0.4721992 0.08586016 NA
## Test set 1.2258799 0.30190797 0.3048635#Based on the accuracy, the STLF model has the best results for all three dataset.
###8.11 Exercise #18
#I choose the dataset passenger-miles-mil-flown-domest
library(fpp)
library(forecast)
library(fma)
library(expsmooth)
library(lmtest)
library(zoo)
#Loading Dataset
passenger.ts <- ts(rdatamarket::dmseries("http://bit.ly/209Fr0e")[,1],
start=c(1962,7), frequency=12)
plot(passenger.ts,xlab="Month",ylab="miles flown(millions)",main= "Passenger miles flown domestic U.K.1962-1972")
#ARIMA
Arima1=auto.arima(passenger.ts)
accuracy(Arima1)## ME RMSE MAE MPE MAPE MASE
## Training set -0.3408052 4.470344 3.420576 -0.738475 3.959842 0.5041345
## ACF1
## Training set -0.02803551Arima2=Arima(passenger.ts,order=c(1,1,1), seasonal=c(0,1,0))
accuracy(Arima2)## ME RMSE MAE MPE MAPE MASE
## Training set -0.1657602 4.87375 3.633692 -0.43681 4.182477 0.5355441
## ACF1
## Training set -0.01291758#Check for White Noise
tsdisplay(residuals(Arima1))
tsdisplay(residuals(Arima2))
#Forecast
forecast1<-forecast(Arima1, h=48)
plot(forecast1)
#ETS
ets1<-ets(passenger.ts,model="ZZZ")
ets1## ETS(A,Ad,A)
##
## Call:
## ets(y = passenger.ts, model = "ZZZ")
##
## Smoothing parameters:
## alpha = 0.7281
## beta = 1e-04
## gamma = 1e-04
## phi = 0.9788
##
## Initial states:
## l = 55.651
## b = 1.1837
## s=26.4599 10 -5.4469 -18.8031 -34.9832 -33.0782
## -30.3075 -31.1193 -7.665 27.3536 47.4661 50.1235
##
## sigma: 4.2322
##
## AIC AICc BIC
## 948.0833 954.9233 998.1075tsdisplay(residuals(ets1))
#ETS Forecast
forecastets<-forecast(ets1, h=48)
plot(forecastets)
#Check for accuracy
accuracy(Arima1)## ME RMSE MAE MPE MAPE MASE
## Training set -0.3408052 4.470344 3.420576 -0.738475 3.959842 0.5041345
## ACF1
## Training set -0.02803551accuracy(ets1)## ME RMSE MAE MPE MAPE MASE
## Training set 0.1104657 4.232196 3.313952 0.09477059 4.026201 0.4884199
## ACF1
## Training set 0.01556275#ETS model outperfrom based on the accuracy