MODEL ARIMA

library(forecast)

## Warning: package 'forecast' was built under R version 4.4.2

## 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(tseries)
library(graphics)
library(readxl)
library(caret)

## Warning: package 'caret' was built under R version 4.4.2

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 4.4.2

## Loading required package: lattice

## Warning: package 'lattice' was built under R version 4.4.2

library(lattice)

Memanggil Data

data<- read_xlsx("C:/Users/ASUS/Documents/mobil_toyota.xlsx")
data

## # A tibble: 149 × 2
##    Bulan   Penjualan
##    <chr>       <dbl>
##  1 2011-01     27619
##  2 2011-02     25532
##  3 2011-03     32275
##  4 2011-04     21128
##  5 2011-05     19554
##  6 2011-06     26175
##  7 2011-07     30149
##  8 2011-08     25207
##  9 2011-09     30655
## 10 2011-10     31109
## # ℹ 139 more rows

Mengubah karakter menjadi numerik

data$Bulan <- as.numeric(gsub("-", "", data$Bulan))
print(data)

## # A tibble: 149 × 2
##     Bulan Penjualan
##     <dbl>     <dbl>
##  1 201101     27619
##  2 201102     25532
##  3 201103     32275
##  4 201104     21128
##  5 201105     19554
##  6 201106     26175
##  7 201107     30149
##  8 201108     25207
##  9 201109     30655
## 10 201110     31109
## # ℹ 139 more rows

Statistika Deskriptif

summary(data)

##      Bulan          Penjualan    
##  Min.   :201101   Min.   :  695  
##  1st Qu.:201402   1st Qu.:25081  
##  Median :201703   Median :29362  
##  Mean   :201678   Mean   :28403  
##  3rd Qu.:202004   3rd Qu.:33344  
##  Max.   :202305   Max.   :40781

Split Data

set.seed(123)
split_index <- sample.int(nrow(data), size = floor(0.75 * nrow(data)), replace = FALSE)
data.train <- data[split_index, ]
data.test <- data[-split_index, ]

print(nrow(data.train))

## [1] 111

print(nrow(data.test))

## [1] 38

print(data.train)

## # A tibble: 111 × 2
##     Bulan Penjualan
##     <dbl>     <dbl>
##  1 201202     33558
##  2 201502     26743
##  3 202010     16345
##  4 201407     28757
##  5 202305     28178
##  6 202303     29471
##  7 201806     18642
##  8 201807     34984
##  9 202210     33740
## 10 201808     31149
## # ℹ 101 more rows

print(data.test)

## # A tibble: 38 × 2
##     Bulan Penjualan
##     <dbl>     <dbl>
##  1 201101     27619
##  2 201102     25532
##  3 201103     32275
##  4 201105     19554
##  5 201111     15195
##  6 201206     37176
##  7 201207     36353
##  8 201304     39668
##  9 201309     40235
## 10 201312     34819
## # ℹ 28 more rows

Mengubah Data Menjadi Time Series

ts.data<-ts(data$Penjualan, frequency=12, start=c(2011,1))
print(ts.data)

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2011 27619 25532 32275 21128 19554 26175 30149 25207 30655 31109 15195 26076
## 2012 29189 33558 33306 34264 34737 37176 36353 25848 33705 35855 36996 34427
## 2013 35923 35318 32726 39668 36282 35125 39210 24899 40235 39246 40781 34819
## 2014 35886 38631 38959 39323 31440 39107 28757 30273 29247 31538 28835 27123
## 2015 27166 26743 31801 30053 23223 23967 13495 31991 32405 29826 29607 21541
## 2016 24859 25447 29995 31059 31662 31519 24718 38803 36994 34116 40510 31888
## 2017 32377 36147 38656 33328 32806 21974 31621 34690 27251 29979 28797 23706
## 2018 25405 27665 31424 29360 28950 18642 34984 31149 29821 36119 31981 26661
## 2019 25081 23443 28732 29408 29103 18541 29362 28929 31831 30944 28970 27453
## 2020 24119 25053 26191  2053   695  3705  7224  8673 13150 16345 15361 18687
## 2021 16033 15144 26258 23301 18167 23271 21763 29898 34046 21730 33450 32707
## 2022 22886 24865 33344 27779 13297 27290 29326 30844 33449 33740 26462 28128
## 2023 28970 27336 29471 21518 28178

Plot Data Time Series

ts.plot(ts.data, xlab="Tahun ", ylab="Penjualan", main= "Time Series Plot Data Penjualan Mobil Toyota")
points(ts.data)

Trend Jangka Panjang

decomp <- decompose(ts.data)
plot(decomp$trend, main="Tren Jangka Panjang Penjualan Mobil Toyota", ylab="Nilai", xlab="Tahun")

Pola Musiman

seasonal.plot <- ggseasonplot(ts.data, main="Pola Musiman Penjualan", 
                            year.labels=TRUE, year.labels.left=TRUE) +
  ylab("Nilai") + xlab("Bulan")
print(seasonal.plot)

Partisi Data

ts.datatrain <- ts(data.train$Penjualan, frequency=12, start=c(2011,1))
ts.datatest <- ts(data.test$Penjualan, frequency=12, start=c(2011,1))
print(ts.datatrain)

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2011 33558 26743 16345 28757 28178 29471 18642 34984 33740 31149 13297 28732
## 2012 31888 35318 30149 21974 27251 27336 29362 13150 33328 33306 24899 30944
## 2013 24119 33344 30655 31440 36147 36996 32726 21541 23223 23301 21730 15361
## 2014 26661 38631 28950 39246 29821 36994 29326 34046 29995 29189 29979 25081
## 2015 33449 35923  2053 33705 31621 15144 28835 26462 28128 34264 26175 28970
## 2016 27665 26191 38959 39210 22886 21518 25053 21128 28929 27453 29408 30053
## 2017 35855 21763 18541 40781 39323 35125 26076 29360 36119 31059 28970 18687
## 2018 35886 25207 31801 31109 23706 39107 30273 28797 25447 32806 40510 25848
## 2019 31538 34737  7224 16033 32377 31831 34690 29898 34427 31519 32707  8673
## 2020 36282 27123 29607

Plot Data Partisi

plot.ts(ts.datatrain, main="Data Training Penjualan Mobil Toyota", xlab="Periode", ylab="Nilai Penjualan")

plot.ts(ts.datatest, main="Data Testing Penjualan Mobil Toyota", xlab="Periode", ylab="Nilai Penjualan")

Uji Stasioneritas

adf.test(ts.datatrain)

## Warning in adf.test(ts.datatrain): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts.datatrain
## Dickey-Fuller = -5.0467, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

ACF dan PACF Plot

acf(ts.datatrain, main="ACF Plot", lag.max=10)

pacf(ts.datatrain, main="PACF Plot", lag.max=10)

EACF

eacf(ts.datatrain)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o 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 o x o o o o o o o o o  o  o  o 
## 4 x x o x o o o o o o o  o  o  o 
## 5 o x o x x o o o o o o  o  o  o 
## 6 x x o x o x o o o o o  o  o  o 
## 7 o x o x o x o o o o o  o  o  o

Penentuan Model Terbaik

arima(ts.datatrain, order=c(1,0,0), method="ML")

## 
## Call:
## arima(x = ts.datatrain, order = c(1, 0, 0), method = "ML")
## 
## Coefficients:
##           ar1   intercept
##       -0.0245  28692.5220
## s.e.   0.0947    657.7485
## 
## sigma^2 estimated as 50376532:  log likelihood = -1141.8,  aic = 2287.59

arima(ts.datatrain, order=c(1,0,1), method="ML")

## 
## Call:
## arima(x = ts.datatrain, order = c(1, 0, 1), method = "ML")
## 
## Coefficients:
##         ar1      ma1   intercept
##       0.859  -1.0000  28794.6017
## s.e.  0.054   0.0351    122.3568
## 
## sigma^2 estimated as 46895129:  log likelihood = -1138.94,  aic = 2283.89

arima(ts.datatrain, order=c(2,0,1), method="ML")

## 
## Call:
## arima(x = ts.datatrain, order = c(2, 0, 1), method = "ML")
## 
## Coefficients:
##           ar1      ar2     ma1   intercept
##       -0.2095  -0.1316  0.1852  28690.4243
## s.e.   1.5781   0.0935  1.6016    591.3853
## 
## sigma^2 estimated as 49511446:  log likelihood = -1140.85,  aic = 2289.71

ARIMA.data <- Arima(ts.datatrain, order=c(1,0,1), method="ML")

Eksploratif

arima<-arima(ts.datatrain,   order=c(1,0,1), include.mean = TRUE,method="ML")
sisaan <- arima$residuals

par(mfrow=c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(c(1:length(sisaan)),sisaan)
acf(sisaan)
pacf(sisaan)

Ljung-Box Test

Box.test(sisaan, lag = 12 ,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 8.904, df = 12, p-value = 0.7111

Periode Prediksi

tahun.awal <- 2023
jumlah.periode <- 6
periode <- seq(6, jumlah.periode, by = 1)
tahun <- tahun.awal + floor((periode - 1) / 12)
angka_periode <- ((periode - 1) %% 12) + 1
periode_baru <- tahun * 100 + angka_periode
prediksi.periode <- data.frame(Bulan = periode_baru)

Forecast

forecast <- forecast(ARIMA.data, h=1)

Mengambil Nilai Forecast

nilai.forecast <- forecast$mean

Hasil Forecast

hasil.forecast <- data.frame(Periode = prediksi.periode, Forecast = nilai.forecast)
print(hasil.forecast)

##    Bulan Forecast
## 1 202306 27837.23

Plot Forecast

plot(forecast)

Accuracy Matrix

MPE.arima <- mean(abs((forecast$mean - ts.datatest[length(ts.datatest)]) / ts.datatest[length(ts.datatest)]) * 100)
print(MPE.arima)

## [1] 9.748311