library (readxl)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(TTR)
library(TSA)
## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
library(graphics)
library(tseries)
dataawal<- read_excel("Data Time Series.xlsx")
dataawal
## # A tibble: 120 × 1
## `Sales (in thousands)`
## <dbl>
## 1 10618.
## 2 10538.
## 3 10209.
## 4 10553
## 5 9935.
## 6 10534.
## 7 10196.
## 8 10512.
## 9 10090.
## 10 10371.
## # ℹ 110 more rows
data1<-dataawal[-111:-120,]
data2<-dataawal[-1:-110,]
data.ts<-ts(data1)
head(data.ts)
## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## Sales (in thousands)
## [1,] 10618.1
## [2,] 10537.9
## [3,] 10209.3
## [4,] 10553.0
## [5,] 9934.9
## [6,] 10534.5
ts.plot(data.ts, xlab="Time Period ", ylab="Sales", main= "Time Series Plot Data Sales")
points(data.ts)
adf.test(data.ts)
## Warning in adf.test(data.ts): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: data.ts
## Dickey-Fuller = -5.9291, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
acf(data.ts)
pacf(data.ts, lag.max = 50)
eacf(data.ts)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o x o o o o o o o o o o o
## 1 o o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x o o o o o o o o o o o o o
## 4 x o o o o o o o o o o o o o
## 5 o o x o o o o o o o o o o o
## 6 x o o o o o o o o o o o o o
## 7 x x o o x o o o o o o o o o
#Penentuan Model Terbaik berdasar AIC
arima(data.ts, order=c(0,0,3), method="ML")
##
## Call:
## arima(x = data.ts, order = c(0, 0, 3), method = "ML")
##
## Coefficients:
## ma1 ma2 ma3 intercept
## 0.0141 -0.0518 -0.2326 10374.4514
## s.e. 0.1002 0.1171 0.0967 14.9337
##
## sigma^2 estimated as 45079: log likelihood = -745.56, aic = 1499.12
arima(data.ts, order=c(3,0,0), method="ML")
##
## Call:
## arima(x = data.ts, order = c(3, 0, 0), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 intercept
## 0.0904 0.0086 -0.2312 10374.9540
## s.e. 0.0927 0.0929 0.0926 17.9023
##
## sigma^2 estimated as 44745: log likelihood = -745.15, aic = 1498.3
arima(data.ts, order=c(1,0,1), method="ML")
##
## Call:
## arima(x = data.ts, order = c(1, 0, 1), method = "ML")
##
## Coefficients:
## ar1 ma1 intercept
## 0.0672 0.0246 10375.1957
## s.e. 0.4398 0.4314 22.7734
##
## sigma^2 estimated as 47350: log likelihood = -748.18, aic = 1502.36
arima(data.ts, order=c(1,0,1), method="ML")
##
## Call:
## arima(x = data.ts, order = c(1, 0, 1), method = "ML")
##
## Coefficients:
## ar1 ma1 intercept
## 0.0672 0.0246 10375.1957
## s.e. 0.4398 0.4314 22.7734
##
## sigma^2 estimated as 47350: log likelihood = -748.18, aic = 1502.36
arima(data.ts, order=c(3,0,1), method="ML")
##
## Call:
## arima(x = data.ts, order = c(3, 0, 1), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ma1 intercept
## 0.6498 -0.0437 -0.2197 -0.6094 10374.0433
## s.e. 0.2868 0.1149 0.1105 0.2950 12.8663
##
## sigma^2 estimated as 43484: log likelihood = -743.64, aic = 1497.28
ARIMA_data <- Arima(data.ts, order=c(3,0,1), method="ML")
hasil_forecastt <- forecast(ARIMA_data)
hasil_forecast <- as.data.frame(hasil_forecastt)
hasil_forecast
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 111 10371.02 10097.49 10644.55 9952.692 10789.35
## 112 10361.97 10088.22 10635.72 9943.305 10780.64
## 113 10362.54 10088.75 10636.33 9943.808 10781.27
## 114 10367.76 10086.66 10648.86 9937.849 10797.67
## 115 10373.12 10088.65 10657.58 9938.063 10808.17
## 116 10376.24 10090.72 10661.76 9939.579 10812.91
## 117 10376.89 10091.37 10662.41 9940.230 10813.56
## 118 10376.00 10090.29 10661.72 9939.037 10812.97
## 119 10374.71 10088.73 10660.69 9937.338 10812.08
## 120 10373.76 10087.68 10659.84 9936.242 10811.29
plot(hasil_forecastt)
data_forecast1 <- hasil_forecast$`Point Forecast`
data_forecast1 <- as.data.frame(data_forecast1)
head(data_forecast1)
## data_forecast1
## 1 10371.02
## 2 10361.97
## 3 10362.54
## 4 10367.76
## 5 10373.12
## 6 10376.24
dataaktual <- as.data.frame(data2)
dataaktual
## Sales (in thousands)
## 1 10210.1
## 2 10352.5
## 3 10423.8
## 4 10519.3
## 5 10596.7
## 6 10650.0
## 7 10741.6
## 8 10246.0
## 9 10354.4
## 10 10155.4
akurasi.arima <- accuracy(data_forecast1$data_forecast1,dataaktual$`Sales (in thousands)`)
akurasi.arima
## ME RMSE MAE MPE MAPE
## Test set 53.57804 193.8149 161.3917 0.4821229 1.538676
plot(data.ts,main="Plot semua data")
points(data.ts)
data1.ses<-HoltWinters(data.ts,alpha=0.1, beta=FALSE, gamma=FALSE)
data1.ses
## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = data.ts, alpha = 0.1, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.1
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 10393.24
SSE <- data1.ses$SSE
SSE
## [1] 5954375
MSE <- SSE/110
MSE
## [1] 54130.69
hasil_forecast.sess <- forecast(data1.ses)
hasil_forecast.ses <- as.data.frame(hasil_forecast.sess)
hasil_forecast.ses
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 111 10393.24 10093.50 10692.98 9934.830 10851.65
## 112 10393.24 10092.01 10694.48 9932.543 10853.94
## 113 10393.24 10090.52 10695.96 9930.268 10856.22
## 114 10393.24 10089.04 10697.44 9928.004 10858.48
## 115 10393.24 10087.57 10698.92 9925.751 10860.73
## 116 10393.24 10086.10 10700.38 9923.509 10862.97
## 117 10393.24 10084.64 10701.84 9921.278 10865.21
## 118 10393.24 10083.19 10703.29 9919.057 10867.43
## 119 10393.24 10081.74 10704.74 9916.846 10869.64
## 120 10393.24 10080.30 10706.18 9914.645 10871.84
plot(hasil_forecast.sess)
akurasi.arima
## ME RMSE MAE MPE MAPE
## Test set 53.57804 193.8149 161.3917 0.4821229 1.538676
akurasi.ses
## ME RMSE MAE MPE MAPE
## Test set 31.73833 189.8569 161.3 0.2723997 1.540696
fit <- Arima(data.ts, order=c(3,0,1), method="ML")
fit
## Series: data.ts
## ARIMA(3,0,1) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ma1 mean
## 0.6498 -0.0437 -0.2197 -0.6094 10374.0433
## s.e. 0.2868 0.1149 0.1105 0.2950 12.8663
##
## sigma^2 = 45555: log likelihood = -743.64
## AIC=1499.28 AICc=1500.1 BIC=1515.49
hasil_forecast
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 111 10371.02 10097.49 10644.55 9952.692 10789.35
## 112 10361.97 10088.22 10635.72 9943.305 10780.64
## 113 10362.54 10088.75 10636.33 9943.808 10781.27
## 114 10367.76 10086.66 10648.86 9937.849 10797.67
## 115 10373.12 10088.65 10657.58 9938.063 10808.17
## 116 10376.24 10090.72 10661.76 9939.579 10812.91
## 117 10376.89 10091.37 10662.41 9940.230 10813.56
## 118 10376.00 10090.29 10661.72 9939.037 10812.97
## 119 10374.71 10088.73 10660.69 9937.338 10812.08
## 120 10373.76 10087.68 10659.84 9936.242 10811.29
interval.prediction <- hasil_forecast[,-1]
interval.prediction
## Lo 80 Hi 80 Lo 95 Hi 95
## 111 10097.49 10644.55 9952.692 10789.35
## 112 10088.22 10635.72 9943.305 10780.64
## 113 10088.75 10636.33 9943.808 10781.27
## 114 10086.66 10648.86 9937.849 10797.67
## 115 10088.65 10657.58 9938.063 10808.17
## 116 10090.72 10661.76 9939.579 10812.91
## 117 10091.37 10662.41 9940.230 10813.56
## 118 10090.29 10661.72 9939.037 10812.97
## 119 10088.73 10660.69 9937.338 10812.08
## 120 10087.68 10659.84 9936.242 10811.29
plot(hasil_forecast$`Point Forecast`, type="n", ylim=range(hasil_forecast$`Lo 80`,hasil_forecast$`Hi 80`))
polygon(c(time(hasil_forecast$`Point Forecast`),rev(time(hasil_forecast$`Point Forecast`))), c(hasil_forecast$`Hi 80`,rev(hasil_forecast$`Lo 80`)),
col=rgb(0,0,0.6,0.2), border=FALSE)
lines(hasil_forecast$`Point Forecast`)
lines(fitted(fit),col='red')
out <- (hasil_forecast$`Point Forecast`< hasil_forecast$`Lo 80` | hasil_forecast$`Point Forecast` > hasil_forecast$`Hi 80`)
plot(hasil_forecast$`Point Forecast`, type="n", ylim=range(hasil_forecast$`Lo 95`,hasil_forecast$`Hi 95`))
polygon(c(time(hasil_forecast$`Point Forecast`),rev(time(hasil_forecast$`Point Forecast`))), c(hasil_forecast$`Hi 95`,rev(hasil_forecast$`Lo 95`)),
col=rgb(0,0,0.6,0.2), border=FALSE)
lines(hasil_forecast$`Point Forecast`)
lines(fitted(fit),col='red')
out <- (hasil_forecast$`Point Forecast`< hasil_forecast$`Lo 95` | hasil_forecast$`Point Forecast` > hasil_forecast$`Hi 95`)
out
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE