library(ggplot2)
library(tsibble)
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(MASS)
library(readxl)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tseries)
library(dplyr)
library(ggplot2)
library(tsibble)
library(tseries)
library(MASS)
library(forecast)
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(TTR)
library(aTSA)
## 
## Attaching package: 'aTSA'
## The following object is masked from 'package:forecast':
## 
##     forecast
## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test
## The following object is masked from 'package:graphics':
## 
##     identify
library(graphics)
data = read_excel("C:\\Users\\MUTHI'AH IFFA\\Downloads\\data_penjualan_ujian.xlsx")
data
## # A tibble: 500 Ɨ 2
##     Hari Penjualan
##    <dbl>     <dbl>
##  1     1       190
##  2     2       176
##  3     3       182
##  4     4       184
##  5     5       182
##  6     6       178
##  7     7       189
##  8     8       176
##  9     9       191
## 10    10       176
## # ā„¹ 490 more rows
str(data)
## tibble [500 Ɨ 2] (S3: tbl_df/tbl/data.frame)
##  $ Hari     : num [1:500] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Penjualan: num [1:500] 190 176 182 184 182 178 189 176 191 176 ...
data.baru <- data %>%
  mutate(Penjualan = Penjualan + 91)

data.baru
## # A tibble: 500 Ɨ 2
##     Hari Penjualan
##    <dbl>     <dbl>
##  1     1       281
##  2     2       267
##  3     3       273
##  4     4       275
##  5     5       273
##  6     6       269
##  7     7       280
##  8     8       267
##  9     9       282
## 10    10       267
## # ā„¹ 490 more rows
str(data.baru)
## tibble [500 Ɨ 2] (S3: tbl_df/tbl/data.frame)
##  $ Hari     : num [1:500] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Penjualan: num [1:500] 281 267 273 275 273 269 280 267 282 267 ...

Buat Data Time Series

ts.data <- as.ts(data.baru$Penjualan)
head(ts.data)
## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1] 281 267 273 275 273 269

Eksplorasi Data

# plot deret waktu 
plot(ts.data, 
     main = "Total Sales (Time Series)",
     ylab = "Penjualan", 
     xlab = "Observation ")

#ACF dan PACF
acf(data$Penjualan, main = "ACF Penjualan")

pacf(data$Penjualan, main = "PACF Penjualan")

# differencing  
diff.ts <- diff(ts.data)
head(diff.ts)
## Time Series:
## Start = 2 
## End = 7 
## Frequency = 1 
## [1] -14   6   2  -2  -4  11
plot(diff.ts, type = "l",
     main = "Penjualan Setelah Differencing",
     xlab = "Hari", 
     ylab = "Ī” Penjualan")

adf.test(diff.ts)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -38.1    0.01
## [2,]   1 -26.0    0.01
## [3,]   2 -21.9    0.01
## [4,]   3 -19.6    0.01
## [5,]   4 -14.5    0.01
## [6,]   5 -12.4    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -38.1    0.01
## [2,]   1 -26.1    0.01
## [3,]   2 -22.0    0.01
## [4,]   3 -19.7    0.01
## [5,]   4 -14.7    0.01
## [6,]   5 -12.5    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -38.1    0.01
## [2,]   1 -26.1    0.01
## [3,]   2 -22.1    0.01
## [4,]   3 -19.8    0.01
## [5,]   4 -14.8    0.01
## [6,]   5 -12.7    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Cek Autokorelasi

library(lmtest) 
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:tsibble':
## 
##     index
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(orcutt) 
library(HoRM)   
model = lm(Penjualan~Hari, data = data.baru)
summary(model)
## 
## Call:
## lm(formula = Penjualan ~ Hari, data = data.baru)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.281 -18.579  -0.433  17.052  69.022 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.353e+02  1.997e+00  117.83   <2e-16 ***
## Hari        2.585e-01  6.908e-03   37.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.29 on 498 degrees of freedom
## Multiple R-squared:  0.7377, Adjusted R-squared:  0.7372 
## F-statistic:  1401 on 1 and 498 DF,  p-value: < 2.2e-16
# Ambil sisaan dari model
sisaan.awal <- residuals(model)

plot(data$Hari, sisaan.awal, type = "o", pch = 20, col = "red",
     main = "Plot Sisaan vs Waktu",
     xlab = "Hari", ylab = "Sisaan")

dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.28812, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

model terindikasi autokorelasi dengan p-value<0.05

penanganannya adalah dengan metode:

  1. Cochrane-Orcutt = melakukan iteratif untuk menduga nilai rho dan transformasi data untuk hilangkan autokorelasi
  2. Hidreth-Lu =Ia akan mencoba berbagai kemungkinan nilai Ļ dan memilih satu yang menghasilkan Sum of Squared Errors (SSE) terkecil.

Peramalan

# Data latih 
train = data.baru$Penjualan[1:400]
train.ts = ts(train)
plot.ts(train, 
        lty = 1, 
        xlab = "Waktu", 
        ylab = "Penjualan", 
        main = "Plot Penjualan Train")

#data uji
test = data.baru$Penjualan[401:500]
test.ts = ts(test)
plot.ts(test, 
        lty = 1, 
        xlab = "Waktu",
        ylab = "Penjualan",
        main = "Plot Penjualan Test")

# acf
acf(test)

# Plot ACF 
acf(train)

#PACF Test
pacf(test)

#PACF 
pacf(train)

tseries::adf.test(train)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train
## Dickey-Fuller = -1.0344, Lag order = 7, p-value = 0.9325
## alternative hypothesis: stationary
eacf(train)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 x o o o o o o x o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x x x o o o o o o o o  o  o  o 
## 4 x x x x o o o o o o o  o  o  o 
## 5 x o o x x o o o o o o  o  o  o 
## 6 x o x x x x o o o o o  o  o  o 
## 7 x x x x o x x o o o o  o  o  o

identifikasi

AR

MA(2) karena cut off pada lag ke-2

ARMA (1,1,1) karena differencing sebanyak 1 kali

# ARIMA (1,1,1)
model1.da = Arima(train, order = c(1,1,1), method = "ML")
summary(model1.da)
## Series: train 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.0490  -0.7304
## s.e.   0.0601   0.0345
## 
## sigma^2 = 69.14:  log likelihood = -1410.69
## AIC=2827.37   AICc=2827.43   BIC=2839.34
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.639026 8.283944 6.511328 0.1330468 2.310579 0.7709258
##                      ACF1
## Training set -0.009635835
lmtest::coeftest(model1.da)
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1 -0.048993   0.060069  -0.8156   0.4147    
## ma1 -0.730376   0.034516 -21.1604   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
prediksi_30 <- forecast::forecast(model1.da, h = 30)
print(prediksi_30)
##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 401       341.8420 331.1857 352.4983 325.5445 358.1395
## 402       342.0947 331.1821 353.0073 325.4053 358.7841
## 403       342.0823 330.8266 353.3381 324.8681 359.2965
## 404       342.0829 330.4989 353.6669 324.3667 359.7991
## 405       342.0829 330.1795 353.9863 323.8782 360.2876
## 406       342.0829 329.8684 354.2974 323.4025 360.7633
## 407       342.0829 329.5651 354.6007 322.9385 361.2273
## 408       342.0829 329.2689 354.8969 322.4856 361.6802
## 409       342.0829 328.9795 355.1863 322.0429 362.1229
## 410       342.0829 328.6962 355.4695 321.6098 362.5560
## 411       342.0829 328.4189 355.7469 321.1856 362.9802
## 412       342.0829 328.1471 356.0187 320.7699 363.3959
## 413       342.0829 327.8805 356.2853 320.3622 363.8036
## 414       342.0829 327.6188 356.5470 319.9619 364.2039
## 415       342.0829 327.3617 356.8041 319.5688 364.5970
## 416       342.0829 327.1091 357.0567 319.1824 364.9834
## 417       342.0829 326.8606 357.3052 318.8024 365.3634
## 418       342.0829 326.6162 357.5496 318.4286 365.7372
## 419       342.0829 326.3755 357.7903 318.0605 366.1053
## 420       342.0829 326.1385 358.0273 317.6980 366.4678
## 421       342.0829 325.9049 358.2609 317.3408 366.8250
## 422       342.0829 325.6747 358.4911 316.9887 367.1771
## 423       342.0829 325.4477 358.7181 316.6415 367.5243
## 424       342.0829 325.2237 358.9421 316.2990 367.8668
## 425       342.0829 325.0026 359.1632 315.9609 368.2049
## 426       342.0829 324.7844 359.3814 315.6272 368.5386
## 427       342.0829 324.5689 359.5969 315.2976 368.8682
## 428       342.0829 324.3560 359.8098 314.9720 369.1938
## 429       342.0829 324.1457 360.0201 314.6503 369.5155
## 430       342.0829 323.9378 360.2280 314.3323 369.8335
plot(prediksi_30,
     main = "Peramalan Penjualan 30 Hari ke Depan",
     xlab = "Hari",
     ylab = "Penjualan")