Packages

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")

Data

Data berupa penjualan produk farmasi (dalam ribuan) yang diobservasi setiap minggu sebanyak 120 kali. Berikut eksplorasi data tersebut

## # A tibble: 6 x 2
##    Week  Sales
##   <dbl>  <dbl>
## 1     1 10618.
## 2     2 10538.
## 3     3 10209.
## 4     4 10553 
## 5     5  9935.
## 6     6 10534.

Pemeriksaan Kestasioneran Y

#Berdasarkan plot y, secara subjektif data terlihat stasioner
plot.ts(y, lty=1, xlab="Minggu ke-", ylab="Penjualan(Ribuan)", main="Plot Y")
points(y)

#Plot ACF sebagai kelanjutan pemeriksaan kestasioneran. Secara grafis data terlihat stasioner dalam mean
acf(y, lag.max=10, main="ACF data y")

#Uji ADF sebagai kelanjutan pemeriksaan kestasioneran. Secara hipotesis (p.value<0.05) data telah stasioner
adf.test(y,alternative=c("stationary"),k=trunc((length(y)-1)^(1/3)))
## Warning in adf.test(y, alternative = c("stationary"), k = trunc((length(y) - :
## p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y
## Dickey-Fuller = -6.2859, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
#Data tidak dilakukan pembedaan

Penyesuaian Model

#Kandidat model menurut plot ACF (tidak teridentifikasi cuts off)
acf(y, lag.max=10, main="ACF data y")

#Kandidat model menurut plot PACF (tidak teridentifikasi cuts off)
pacf(y, lag.max=10, main="PACF data y")

#Kandidat model menurut EACF adalah ARMA(1,1)
eacf(y)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o x x 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 x 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 o o o o o o o o o o  o  o  o
#Dapat disimpulkan model kandidat yakni ARMA(1,1)

Model ARMA (1,1)

Berikut dugaan 10 data terakhir

model<-arma(y,order=c(1,1))
## Warning in arma(y, order = c(1, 1)): singular Hessian
dugaan<-fitted.values(model)
cbind(y,dugaan)[111:120,]
##             y   dugaan
##  [1,] 10210.1 10378.24
##  [2,] 10352.5 10358.08
##  [3,] 10423.8 10372.56
##  [4,] 10519.3 10381.38
##  [5,] 10596.7 10392.52
##  [6,] 10650.0 10401.80
##  [7,] 10741.6 10408.30
##  [8,] 10246.0 10418.85
##  [9,] 10354.4 10364.65
## [10,] 10155.4 10373.19

dan plot fitted model

plot.ts(y,xlab="Week",ylab="Penjualan(Ribuan)")
points(y)
par(col="Red")
lines(dugaan)

par(col="black")

Simple Exponential Smoothing

Berikut plot SES dengan \(\lambda=0.1\)

modelses<-HoltWinters(y,alpha=0.1,beta=F,gamma=F)
plot(modelses,xlab="Week",ylab="Penjualan(Ribuan)",main = "SES 0.1")
points(y)

Dugaan data 10 terakhir

dugaant<-fitted(modelses)
ye<-y[110:120]
dugaany<-dugaant[109:119,-2]
cbind(ye,dugaany)
##            ye  dugaany
##  [1,] 10391.3 10393.46
##  [2,] 10210.1 10393.24
##  [3,] 10352.5 10374.93
##  [4,] 10423.8 10372.68
##  [5,] 10519.3 10377.80
##  [6,] 10596.7 10391.95
##  [7,] 10650.0 10412.42
##  [8,] 10741.6 10436.18
##  [9,] 10246.0 10466.72
## [10,] 10354.4 10444.65
## [11,] 10155.4 10435.62

Perbandingan ARMA (1,1) dengan SES lambda=0.1

Indikator perbandingan SSE bahwa model ARMA (1,1) lebih baik (SSE kecil) dari SES pada kasus data tersebut.

#SSE bagi SES
modelses$SSE
## [1] 6338093
#SSE bagi ARIMA (1,1)
sum((y[2:120]-dugaan[2:120])^2)
## [1] 5493730

Prediction Interval

futurVal<-forecast(y,h=10, level=c(95.0))
plot(futurVal)