Tugas Individu STA1341: Regresi Peubah Lag
Afris Setiya Intan Amanda (G1401201018)
2022-09-13
Analisis Regresi dengan Peubah Lag: Penjualan Mempengaruhi Pendapatan Perusahaan (Januari 2015-April 2020)
Latar Belakang
Pada zaman dahulu, analisis runtun waktu menjadi pertentangan antara para ahli ekonometrika dan para ahli runtun waktu. Ahli ekonometrika cenderung menformulasikan model regresi klasik untuk menganalisis perilaku data runtun waktu, menganalisis tentang masalah simultanitas, dan kesalahan autokorelasi. Namun, ahli runtun waktu membuat model perilaku runtun waktu dengan mekanisme sendiri serta tidak begitu memperhatikan peranan peubah bebas (X) dan peubah tak bebas (Y). Perbedaan pendapat ini membuat para ahli ekonometrika mengkaji ulang pendekatannya terutama dalam menganalisis runtun waktu (Jatiningrum 2008).
Tinjauan Pustaka
Analisis Runtun Waktu
Analisis runtun waktu merupakan analisis sekumpulan data dalam suatu periode waktu yang lampau yang digunakan untuk mengetahui atau meramalkan kondisi masa mendatang. Hali ini didasarkan bahwa perilaku manusia banyak dipengaruhi oleh waktu sebelumnya. Oleh katena itu, faktor waktu memiliki peran sangat penting (Handoko 2013).
Data Training dan Testing
Data training merupakan data yang digunakan untuk membangun sebuah model melalui serangkaian algoritma tertentu (tahapan modelling), sedangkan data testing merupakan data yang digunakan untuk mengetahui akurasi performa model data training (tahapan forecasting) (Nasution et al. 2019).
Metode Koyck
Metode Koyck didasari asumsi bahwa semakin jauh jarak lag pada peubah bebas dari periode sekarang maka semakin kecil pengaruh peubah lag terhadap peubah tak bebas. Koyck mengusulkan suatu metode pendugaan model lag terdistribusi didasarkan pada asumsi bahwa koefisien menurun secara eksponensial dari waktu ke waktu (Gujarati 2004). Model Koyck merupakan jenis paling umum dari infinite distributed lag model atau dapat disebut juga dengan geometric lag model.
Model Distribusi Lag
Model regresi yang memuat peubah tak bebas (Y) yang dipengaruhi oleh peubah bebas (X) pada waktu t dan juga peubah bebas pada waktu \(t−1\), \(t−2\), dan seterusnya. Oleh sebab itu, pengaruh dari suatu atau beberapa peubah bebas (X) terhadap peubah tak bebas (Y) menyebar (distributed) ke beberapa periode waktu (Jatiningrum 2008).
Model Autoregressive
Model autoregressive merupakan bentuk regresi yang menghubungkan nilai-nilai sebelumnya pada selang waktu (lag time) yang bermacam-macam. Oleh karena itu, model autoregressive disebut juga dynamic regression, menyatakan suatu ramalan sebagai fungsi nilai-nilai sebelumnya dari data time series periode tertentu (Pawestri et al. 2018).
MAPE
Mean Absolute Percentage Error (MAPE) digunakan untuk merata-ratakan persentase absolut kesalahan dan memberikan gambaran seberapa kesalahan peramalan dibanding dengan nilai sebenarnya. MAPE banyak dipakai dalam perbandingan data-data yang mempunyai skala interval waktu yang berbeda (Robial 2018).
AIC
Akaike information criterion (AIC) merupakan penilaian yang diaplikasikan untuk menentukan model terbaik di antara beberapa model produksi. AIC dapat diaplikasikan sebagai model selektor yang mampu mengaksesi kualitas tiap model secara relatif dengan menggunakan estimasi maximum likelihood sebagai pehitungan yang sesuai (Anderson et al. 1998). Pemilihan model terbaik berdasarkan kriteria AIC dilakukan dengan memilih model yang memiliki nilai AIC terkecil.
BIC
Bayesian Information Criterion (BIC) merupakan kriteria informasi in sampel untuk pemilihan model terbaik dari beberapa kandidat model yang mungkin. BIC sangat efektif untuk memilih model terbaik dan dapat mengatasi komplektifitas model karena BIC memberikan penalti atas penambahan parameter dan cocok untuk data sampel yang besar. Model terbaik adalah model yang memiliki nilai BIC terkecil (Qudratullah dan Farhan 2007).
Analisis Data
Pemanggilan Packages & Library
Packages yang digunakan dalam analisis regresi dengan peubah lag,
sebagai berikut:
1. dLagM
2. dynlm
3. MLmetrics
4. car
5. corrplot
6. dplyr
7. tidyverse
8. lmtest
Deklarasi & Deskripsi Data
Data yang digunakan analisis ini adalah data sekunder bulanan yang
bersumber dari wesite kaggle mengenai penjualan dan
pendapataan produksi sebuah perusahaan periode Januari 2015-April 2020.
Data ini terdiri dari empat peubah. Seleksi peubah hanya memilih peubah
X yaitu peubah Sales_quantity berupa peubah numerik dan
Revenue sebagai peubah Y. Data dapat diakses melalui tautan
sebagai berikut: data
pendapatan-penjualan. Data ini dipilih untuk dianalisis lebih lanjut
mengenai meulakukan pemodelan koyck, distibusi lag, dan model
autoregressive, lalu dilakukan perbandingan untuk mendapatkan
model terbaik yang dapat mengatasi gejala autokorelasi.
dt <- read.csv("C:/Users/HP/Documents/MPDW/DATA/Month_Value_1.csv")
View(dt)
str(dt)## 'data.frame': 96 obs. of 5 variables:
## $ Period : chr "01.01.2015" "01.02.2015" "01.03.2015" "01.04.2015" ...
## $ Revenue : num 16010072 15807587 22047146 18814583 14021480 ...
## $ Sales_quantity : int 12729 11636 15922 15227 8620 13160 17254 8642 16144 18135 ...
## $ Average_cost : num 1258 1359 1385 1236 1627 ...
## $ The_average_annual_payroll_of_the_region: int 30024676 30024676 30024676 30024676 30024676 30024676 30024676 30024676 30024676 30024676 ...dt <- na.omit(dt)#Deteksi missing value
sum(is.na(dt))## [1] 0knitr::kable(dt, align = "c")| Period | Revenue | Sales_quantity | Average_cost | The_average_annual_payroll_of_the_region |
|---|---|---|---|---|
| 01.01.2015 | 16010072 | 12729 | 1257.764 | 30024676 |
| 01.02.2015 | 15807587 | 11636 | 1358.507 | 30024676 |
| 01.03.2015 | 22047146 | 15922 | 1384.697 | 30024676 |
| 01.04.2015 | 18814583 | 15227 | 1235.607 | 30024676 |
| 01.05.2015 | 14021480 | 8620 | 1626.622 | 30024676 |
| 01.06.2015 | 16783929 | 13160 | 1275.375 | 30024676 |
| 01.07.2015 | 19161892 | 17254 | 1110.577 | 30024676 |
| 01.08.2015 | 15204984 | 8642 | 1759.429 | 30024676 |
| 01.09.2015 | 20603940 | 16144 | 1276.260 | 30024676 |
| 01.10.2015 | 20992875 | 18135 | 1157.589 | 30024676 |
| 01.11.2015 | 14993370 | 10841 | 1383.025 | 30024676 |
| 01.12.2015 | 27791808 | 22113 | 1256.809 | 30024676 |
| 01.01.2016 | 28601586 | 15365 | 1861.477 | 27828571 |
| 01.02.2016 | 22367074 | 13153 | 1700.530 | 27828571 |
| 01.03.2016 | 29738609 | 18339 | 1621.605 | 27828571 |
| 01.04.2016 | 28351008 | 13909 | 2038.321 | 27828571 |
| 01.05.2016 | 15264604 | 8553 | 1784.708 | 27828571 |
| 01.06.2016 | 24385658 | 15101 | 1614.837 | 27828571 |
| 01.07.2016 | 29486517 | 15695 | 1878.720 | 27828571 |
| 01.08.2016 | 15270117 | 8314 | 1836.675 | 27828571 |
| 01.09.2016 | 36141028 | 17764 | 2034.510 | 27828571 |
| 01.10.2016 | 27915144 | 18969 | 1471.619 | 27828571 |
| 01.11.2016 | 21272049 | 13433 | 1583.567 | 27828571 |
| 01.12.2016 | 42014160 | 27029 | 1554.410 | 27828571 |
| 01.01.2017 | 36007381 | 16889 | 2132.002 | 27406473 |
| 01.02.2017 | 30396775 | 15864 | 1916.085 | 27406473 |
| 01.03.2017 | 47678131 | 22786 | 2092.431 | 27406473 |
| 01.04.2017 | 27013965 | 17910 | 1508.317 | 27406473 |
| 01.05.2017 | 24948845 | 10777 | 2315.008 | 27406473 |
| 01.06.2017 | 31101346 | 18799 | 1654.415 | 27406473 |
| 01.07.2017 | 33848822 | 17899 | 1891.101 | 27406473 |
| 01.08.2017 | 16454667 | 9649 | 1705.324 | 27406473 |
| 01.09.2017 | 31650093 | 20159 | 1570.023 | 27406473 |
| 01.10.2017 | 31572206 | 19519 | 1617.511 | 27406473 |
| 01.11.2017 | 22446371 | 15360 | 1461.352 | 27406473 |
| 01.12.2017 | 44966126 | 30833 | 1458.377 | 27406473 |
| 01.01.2018 | 44067521 | 19812 | 2224.284 | 28197847 |
| 01.02.2018 | 36020287 | 18424 | 1955.074 | 28197847 |
| 01.03.2018 | 46995990 | 29004 | 1620.328 | 28197847 |
| 01.04.2018 | 35536488 | 22033 | 1612.876 | 28197847 |
| 01.05.2018 | 29699599 | 14959 | 1985.400 | 28197847 |
| 01.06.2018 | 33261065 | 23067 | 1441.933 | 28197847 |
| 01.07.2018 | 35826535 | 18397 | 1947.412 | 28197847 |
| 01.08.2018 | 23268655 | 12045 | 1931.810 | 28197847 |
| 01.09.2018 | 35423490 | 23358 | 1516.546 | 28197847 |
| 01.10.2018 | 39831566 | 22644 | 1759.034 | 28197847 |
| 01.11.2018 | 32999145 | 19765 | 1669.575 | 28197847 |
| 01.12.2018 | 47221828 | 33207 | 1422.044 | 28197847 |
| 01.01.2019 | 36459960 | 24096 | 1513.113 | 29878525 |
| 01.02.2019 | 36546499 | 21624 | 1690.090 | 29878525 |
| 01.03.2019 | 54198707 | 33379 | 1623.737 | 29878525 |
| 01.04.2019 | 32743990 | 22265 | 1470.649 | 29878525 |
| 01.05.2019 | 32531658 | 16967 | 1917.349 | 29878525 |
| 01.06.2019 | 47709702 | 24958 | 1911.600 | 29878525 |
| 01.07.2019 | 45992142 | 21917 | 2098.469 | 29878525 |
| 01.08.2019 | 36933665 | 14431 | 2559.328 | 29878525 |
| 01.09.2019 | 48526260 | 23253 | 2086.882 | 29878525 |
| 01.10.2019 | 44160416 | 26603 | 1659.979 | 29878525 |
| 01.11.2019 | 36374956 | 21987 | 1654.385 | 29878525 |
| 01.12.2019 | 58756474 | 38069 | 1543.420 | 29878525 |
| 01.01.2020 | 56288301 | 27184 | 2070.641 | 29044998 |
| 01.02.2020 | 40225243 | 23509 | 1711.057 | 29044998 |
| 01.03.2020 | 50022165 | 32569 | 1535.883 | 29044998 |
| 01.04.2020 | 52320693 | 26615 | 1965.835 | 29044998 |
Eksplorasi Data
Time Series Plot
yt <- dt$Revenue
xt <- dt$Sales_quantity
#TS: Pendapatan Perusahaan
ts.p1 <-ts(yt)
plot(ts.p1, lwd=2,type='o', main = "Time Series Plot: Pendapatan Perusahaan", ylab="Pendapatan")
points(ts.p1)
Berdasarkan time series plot di atas, pendapatan sebuah perusahaan cenderung meningkat serta menurun yang tidak beraturan dari waktu ke waktu dan terjadi lonjakan yang cukup signifikan. Oleh karena itu, pola data yang terbentuk dari data cenderung tidak stasioner dan mengalami siklis dan tren positif.
#TS: Penjualan Perusahaan
ts.p2 <-ts(xt)
plot(ts.p1, lwd=2,type='o', main = "Time Series Plot: Penjualan Perusahaan", ylab="Penjualan")
points(ts.p2)
Berdasarkan time series plot di atas, sama halnya pada plot sebelumnya, data penjulaan sebuah perusahaan cenderung berpola tidak stasioner dan mengalami siklis dan tren positif.
Scatter plot
#Penjualan vs pendapatan
plot(xt, yt, col="navyblue", lwd=2, type="p",xlab="Penjualan",ylab="Pendapatan",main="Scatterplot: Penjualan vs Pendapatan")
abline(lm(yt ~ xt), col = "red", lwd = 3)
Cross-Correlation Function Plot
ccf(xt,yt)
Cross-correlation function (CCF) menyatakan korelasi antara pendapatan dengan penjualan sebuah perusahaan yang melintasi waktu. Terindikasi bahwa kedua peubah berkorelasi positif pada beberapa lag. Hal ini dapat terlihat dari lag k = 0, 1, 2, 3, 6, 9, 12, dan 15 yang berada di atas batas signifikan. Lag k = 0 menyatakan korelasi antara pendapatan dengan penjualan sebuah perusahaan pada bulan yang sama. Lag-lag tersebut menyatakan korelasi antara pendapatan k periode ke depan dengan penjualan k periode sebelumnya.
Nilai Korelasi
cor(xt, yt)## [1] 0.8875704Peubah Revenue dengan peubah Sales_quantity
memiliki hubungan linier positif kuat karena nilai
korelasi kedua peubah sebesar 0.8875704 dimana mendekati 1. Artinya,
dari waktu ke waktu (per bulan), semakin tinggi penjualan sebuah
perusahaan, maka semkain tinggi juga pendapatannya. Selain itu,
interpretasi ini didukung oleh analisis secara eksploratif, time
series plot, yang berpola tren positif dimana tebaran data bergerak
dari kiri bawah ke kanan atas dan juga sama halnya dengan
Cross-correlation function plot.
Splitting dan Setting Data
Terdapat observasi yang missing value dengan jumlahnya masih
sekitar 30% dari data keseluruhan, sehingga observasi tersebut dapat
dihapus (omitted). Data terlebih dahulu di-setting
menjadi data time series, sebelum dilakukan pemisahan menjadi
data training dan data testing.
Data training digunakan sebesar 48/64 (75%)
dari keseluruhan amatan data dan data testing yang digunakan
sebesar 16/64 (25%) dari keseluruhan amatan data.
dt <- na.omit(dt)
train <- dt[1:48,]
test <- dt[49:64,]
#data time series
train.ts <- ts(train)
test.ts <- ts(test)
data.ts <- ts(dt)Model Koyck
model.koyck <- dLagM::koyckDlm(x = train$Sales_quantity, y = train$Revenue)
summary(model.koyck)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -8526324 -3486529 -312743 3129736 13059304
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.848e+06 1.923e+07 0.148 0.8829
## Y.1 3.246e-01 1.621e-01 2.002 0.0515 .
## X.t 9.613e+02 1.330e+03 0.723 0.4735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4809000 on 44 degrees of freedom
## Multiple R-Squared: 0.7579, Adjusted R-squared: 0.7469
## Wald test: 16.09 on 2 and 44 DF, p-value: 5.7e-06
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: 4217003 961.3282 0.3245635
Alt text
Berdasarkan persamaan model Koyck di atas, koefisien β dari lag 0 sampai lag 2 berpengaruh positif. Artinya, semakin meningkat penjualan sebuah perusahaan dari bulan ke bulan yang semakin lama, maka dugaan rata-rata pendapatan juga akan meningkat secara signifikan.
#Nilai AIC dan BIC Model Koyck
AIC(model.koyck)## [1] 1584.555BIC(model.koyck)## [1] 1591.956#Forecasting Koyck
(fore.koyck <- forecast(model = model.koyck, x=test$Sales_quantity, h=16))## $forecasts
## [1] 41338966 37053199 46962609 39494640 31977693 37219940 35997983 28404878
## [9] 34421270 39594421 36835946 51400725 45663865 40269008 47227667 43762446
##
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Sales_quantity,
## h = 16)
##
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"Output di atas merupakan nilai forecasting model Koyck untuk meramalkan 16 bulan ke depan.
#MAPE data testing Koyck
mape.koyck <- MAPE(fore.koyck$forecasts, test$Revenue)
#MAPEdata training Koyck
mape_train <- dLagM::GoF(model.koyck)["MAPE"]
c("MAPE_testing" = mape.koyck, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.1321049
##
## $MAPE_training.MAPE
## [1] 0.1517361Model Distribusi Lag
DLM: OLS (lag=2)
#Estimasi parameter menggunakan ordinary least square (ols)
model.dlm = dLagM::dlm(x = train$Sales_quantity, y = train$Revenue, q = 2)
summary(model.dlm)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7048662 -2810627 -116322 2280240 11951227
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5379707.4 3360124.2 -1.601 0.11686
## x.t 1383.9 115.5 11.984 3.87e-15 ***
## x.1 426.8 125.1 3.412 0.00144 **
## x.2 165.9 124.9 1.328 0.19130
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4348000 on 42 degrees of freedom
## Multiple R-squared: 0.8032, Adjusted R-squared: 0.7892
## F-statistic: 57.15 on 3 and 42 DF, p-value: 7.066e-15
##
## AIC and BIC values for the model:
## AIC BIC
## 1 1542.589 1551.732
Alt text
Berdasarkan persamaan model distribusi lag (lag=2) di atas, koefisien β dari lag 0 sampai lag 2 berpengaruh positif. Artinya, semakin meningkat penjualan sebuah perusahaan dari bulan ke bulan yang semakin lama, maka dugaan rata-rata pendapatan juga akan meningkat secara signifikan.
#Nilai AIC dan BIC model distribusi lag
AIC(model.dlm)## [1] 1542.589BIC(model.dlm)## [1] 1551.732#Forecasting DLM OLs
(fore.dlm <- forecast(model = model.dlm, x=test$Sales_quantity, h=16))## $forecasts
## [1] 45419614 40340433 54040679 43267669 33142795 40095565 38418908 28087370
## [9] 36595879 43755097 40261044 61101954 52136978 45073897 54236769 49254536
##
## $call
## forecast.dlm(model = model.dlm, x = test$Sales_quantity, h = 16)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"Output di atas merupakan nilai forecasting model distribusi lag dengan lag=2 untuk meramalkan 16 bulan ke depan.
#MAPE data testing DLM OLS
mape.dlm <- MAPE(fore.dlm$forecasts, test$Revenue)
#MAPE data training DLM OLs
mape_train <- GoF(model.dlm)["MAPE"]
c("MAPE_testing" = mape.dlm, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.1247095
##
## $MAPE_training.MAPE
## [1] 0.1192191DLM: Lag Optimum (Lag=4)
#Penentuan lag optimum
finiteDLMauto(formula = Revenue ~ Sales_quantity,
data = data.frame(train), q.min = 1, q.max = 4,
model.type = "dlm", error.type = "AIC", trace= TRUE)## q - k MASE AIC BIC GMRAE MBRAE R.Adj.Sq Ljung-Box
## 4 4 0.34315 1471.639 1484.128 0.33811 -0.18777 0.81664 0.3510478
## 3 3 0.34468 1503.788 1514.628 0.31416 0.75893 0.81854 0.2916490
## 2 2 0.39416 1542.589 1551.732 0.38225 -0.18747 0.78917 0.3731823
## 1 1 0.40860 1575.311 1582.711 0.41880 0.00501 0.79208 0.2074719Lag optimum ditentukan dengan melihat niali AIC dan BIC terkecil, tetapi nilai tersebut terus menurun seiring pertambahan range jumlah lag. Oleh karena itu, penentuan range lag (q) di antara 1-4. Dengan demikian, laq optimumnya ialah q=4 yaitu pada 4 bulan sebelumnya.
#Model DLM: lag optimum
model.dlm2 = dLagM::dlm(x = train$Sales_quantity, y = train$Revenue, q = 4)
summary(model.dlm2)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7708223 -2269385 -395873 1993129 11261592
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6386571.8 3537621.3 -1.805 0.07895 .
## x.t 1021.6 171.0 5.976 6.15e-07 ***
## x.1 516.7 181.8 2.842 0.00718 **
## x.2 151.8 119.0 1.275 0.21011
## x.3 501.9 180.4 2.783 0.00835 **
## x.4 -133.2 189.2 -0.704 0.48579
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4067000 on 38 degrees of freedom
## Multiple R-squared: 0.838, Adjusted R-squared: 0.8166
## F-statistic: 39.3 on 5 and 38 DF, p-value: 5.127e-14
##
## AIC and BIC values for the model:
## AIC BIC
## 1 1471.639 1484.128
Alt text
Berdasarkan persamaan model distribusi lag optimum (lag=4) di atas, koefisien β dari lag 0 sampai lag 3 berpengaruh positif, kecuali lag 4. Artinya, untuk koefisen lag 0, 1, 2, dan 3, semakin meningkat penjualan sebuah perusahaan dari bulan ke bulan yang semakin lama, maka dugaan rata-rata pendapatan juga akan meningkat secara signifikan. Namun, koefisien di lag 4 bernilai negatif, maka semakin meningkat penjualan sebuah perusahaan pada 4 bulan sebelumnya, akan menurunkan dugaan rata-rata pendapatan pada periode tersebut.
#Nilai AIC dan BIC DLM lag optimum
AIC(model.dlm2)## [1] 1471.639BIC(model.dlm2)## [1] 1484.128#Forecasting DLM lag optimum
(fore.dlm2 <- forecast(model = model.dlm2, x=test$Sales_quantity, h=16))## $forecasts
## [1] 46643032 40099978 56578764 44560854 35162147 45130519 38205390 29019943
## [9] 38418989 42673304 37675132 57652269 54648258 44947279 59338376 49774824
##
## $call
## forecast.dlm(model = model.dlm2, x = test$Sales_quantity, h = 16)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"Output di atas merupakan nilai forecasting model distribusi lag dengan lag=4 (optimum) untuk meramalkan 16 bulan ke depan.
#MAPE data testing DLM lag optimum
mape.dlm2 <- MAPE(fore.dlm2$forecasts, test$Revenue)
#MAPE data training DLM lag optimum
mape_train <- GoF(model.dlm2)["MAPE"]
c("MAPE_testing" = mape.dlm2, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.123609
##
## $MAPE_training.MAPE
## [1] 0.1053313Model Autoregressive
ARDL: p=1 & q=1
#Menggunakan library "dLagM"
model.ardl <- ardlDlm(x = train$Sales_quantity, y = train$Revenue, p = 1 , q = 1)
#p:lag x, q:lag y
#Model untuk p=1, q=1: yt=b0+b1yt-1+b2xt+b3xt-1
summary(model.ardl)##
## Time series regression with "ts" data:
## Start = 2, End = 48
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7583581 -2812209 -219750 2254501 11108666
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.393e+06 2.739e+06 -1.239 0.2221
## X.t 1.375e+03 1.126e+02 12.209 1.45e-15 ***
## X.1 8.142e+01 2.319e+02 0.351 0.7272
## Y.1 2.413e-01 1.306e-01 1.847 0.0716 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4243000 on 43 degrees of freedom
## Multiple R-squared: 0.8157, Adjusted R-squared: 0.8029
## F-statistic: 63.46 on 3 and 43 DF, p-value: 7.771e-16
Alt text
Dari persamaan di atas, koefisien β dari lag 0 sampai lag 1 berpengaruh positif. Artinya, semakin banyak penjualan sebuah perusahaan dari bulan ke bulan yang semakin lama, maka pendapatkan juga akan meningkat secara signifikan. Begitupula, koefisien gamma (γ) bernilai positif, maka semakin banyak pendapatan prusahaan pada 1 bulan sebelumnya, akan meningkatkan pendapatan pada periode tersebut.
#Nilai AIC dan BIC ARDL1
AIC(model.ardl)## [1] 1573.721BIC(model.ardl)## [1] 1582.972#Forecasting ARDL1
(fore.ardl <- forecast(model = model.ardl, x=test$Sales_quantity, h=16))## $forecasts
## [1] 43843137 38886413 53654886 42891427 32103223 40058052 38446181 27514639
## [9] 36399541 43868865 39595961 60305493 51642896 43612287 53834814 48851145
##
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Sales_quantity,
## h = 16)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"Output di atas merupakan nilai forecasting model autoregressive dengan p=1 dan q=1 untuk meramalkan 16 bulan ke depan.
#MAPE testing ARDL1
mape.ardl <- MAPE(fore.ardl$forecasts, test$Revenue) #data testing
#MAPE data training ARDL1
mape_train <- GoF(model.ardl)["MAPE"]
c("MAPE_testing" = mape.ardl, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.1162367
##
## $MAPE_training.MAPE
## [1] 0.1169336#Penentuan lag optimum
(pq.optim <- ardlBoundOrders(data = data.frame(train) , formula = Revenue ~ Sales_quantity))## $p
## Sales_quantity
## 1 12
##
## $q
## [1] 15
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7 q = 8
## p = 1 1539.423 1506.149 1473.930 1443.553 1409.260 1378.191 1346.718 1314.258
## p = 2 1505.139 1506.888 1474.923 1444.514 1410.578 1378.998 1348.158 1316.240
## p = 3 1473.304 1473.304 1475.300 1444.886 1411.383 1379.818 1348.122 1315.563
## p = 4 1442.888 1443.731 1443.731 1445.355 1411.949 1380.386 1349.319 1316.362
## p = 5 1411.321 1412.048 1413.862 1413.862 1406.953 1374.866 1344.321 1313.916
## p = 6 1370.548 1370.081 1370.995 1372.979 1372.979 1374.976 1344.088 1313.437
## p = 7 1344.896 1344.470 1346.047 1347.820 1342.378 1342.378 1344.370 1313.533
## p = 8 1310.856 1307.618 1309.077 1309.617 1309.832 1311.832 1311.832 1313.127
## p = 9 1273.371 1271.992 1272.510 1273.418 1271.506 1271.794 1273.451 1273.451
## p = 10 1241.099 1239.116 1234.701 1236.158 1235.965 1235.151 1237.146 1237.834
## p = 11 1207.689 1207.061 1204.530 1206.487 1207.413 1205.688 1207.684 1207.411
## p = 12 1174.498 1175.159 1172.173 1173.540 1174.975 1172.125 1173.451 1174.751
## p = 13 1137.161 1137.758 1132.308 1134.307 1135.005 1131.655 1128.810 1130.804
## p = 14 1107.164 1107.738 1096.215 1095.716 1097.585 1088.750 1090.051 1091.408
## p = 15 1068.651 1070.228 1065.288 1067.083 1068.236 1057.926 1059.474 1061.452
## q = 9 q = 10 q = 11 q = 12 q = 13 q = 14 q = 15
## p = 1 1270.439 1236.525 1204.387 1174.856 1144.025 1108.114 1078.7286
## p = 2 1271.897 1237.895 1206.192 1176.652 1145.740 1103.387 1074.0663
## p = 3 1273.300 1239.705 1207.912 1178.314 1147.737 1105.284 1075.9351
## p = 4 1273.127 1240.581 1209.143 1179.569 1148.897 1106.953 1077.4509
## p = 5 1272.870 1239.108 1207.749 1178.320 1146.938 1107.513 1077.4872
## p = 6 1272.022 1237.069 1206.992 1177.598 1145.679 1108.403 1079.0119
## p = 7 1272.076 1236.989 1207.228 1177.778 1145.423 1108.731 1079.4488
## p = 8 1273.812 1238.932 1209.165 1179.699 1147.348 1107.760 1078.2867
## p = 9 1273.413 1236.767 1207.313 1177.424 1145.658 1095.315 1062.2816
## p = 10 1237.834 1238.552 1209.240 1179.394 1147.612 1097.252 1064.2802
## p = 11 1209.265 1209.265 1211.236 1181.070 1148.059 1087.712 978.6481
## p = 12 1176.745 1175.288 1175.288 1177.288 1147.545 1088.853 956.7385
## p = 13 1132.673 1134.100 1135.330 1135.330 1136.128 1088.787 NA
## p = 14 1090.808 1091.837 1092.904 1094.723 1094.723 1090.340 NA
## p = 15 1063.366 1058.724 1060.643 NA NA NA NA
##
## $min.Stat
## [1] 956.7385c(p = pq.optim$p$Sales_quantity, q = pq.optim$q)## p q
## 12 15Berdasarkan output di atas, lag optimum bagi peubah
Sales_quantity ialah 12 bulan sebelumnya dan lag optimum
bagi peubah Revenue ialah 15 bulan sebelumnya.
ARDL Lag Optimum: p=12 & q=15
#Menggunakan library "dLagM"
model.ardl2 <- ardlDlm(x = train$Sales_quantity, y = train$Revenue, p = 12 , q = 15)
summary(model.ardl2)##
## Time series regression with "ts" data:
## Start = 16, End = 48
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## 16 17 18 19 20 21 22 23
## 9414 -1216069 -518874 1222241 643550 2297266 -2234424 -1332398
## 24 25 26 27 28 29 30 31
## -1364093 2151237 1005638 1435934 -1113344 -690050 1163257 1453116
## 32 33 34 35 36 37 38 39
## -1277779 -157101 -959945 937891 -891035 -1384140 47811 119313
## 40 41 42 43 44 45 46 47
## 912147 741433 -1726068 -2676872 1793511 -1109758 3010255 -1377705
## 48
## 1085641
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.098e+07 1.660e+07 1.866 0.1355
## X.t 3.417e+03 1.482e+03 2.306 0.0824 .
## X.1 2.026e+03 9.482e+02 2.137 0.0994 .
## X.2 3.689e+02 8.418e+02 0.438 0.6838
## X.3 8.194e+02 7.780e+02 1.053 0.3517
## X.4 -1.043e+03 7.890e+02 -1.322 0.2567
## X.5 5.347e+01 8.410e+02 0.064 0.9524
## X.6 7.830e+02 8.851e+02 0.885 0.4264
## X.7 1.775e+03 8.840e+02 2.008 0.1151
## X.8 2.813e+02 8.678e+02 0.324 0.7620
## X.9 -7.709e+02 9.071e+02 -0.850 0.4433
## X.10 -1.630e+03 9.555e+02 -1.706 0.1631
## X.11 -8.836e+02 1.032e+03 -0.856 0.4402
## X.12 -2.127e+03 1.644e+03 -1.294 0.2653
## Y.1 -5.739e-01 3.693e-01 -1.554 0.1951
## Y.2 -2.598e-01 3.246e-01 -0.801 0.4682
## Y.3 -3.984e-01 4.270e-01 -0.933 0.4036
## Y.4 -7.289e-01 4.037e-01 -1.806 0.1453
## Y.5 -4.618e-01 3.773e-01 -1.224 0.2881
## Y.6 -1.117e-01 3.616e-01 -0.309 0.7728
## Y.7 5.905e-01 3.344e-01 1.766 0.1522
## Y.8 4.223e-02 3.076e-01 0.137 0.8974
## Y.9 -1.369e-01 3.141e-01 -0.436 0.6854
## Y.10 -2.637e-01 3.700e-01 -0.713 0.5154
## Y.11 1.824e-01 3.495e-01 0.522 0.6293
## Y.12 -8.852e-02 3.549e-01 -0.249 0.8153
## Y.13 2.273e-01 3.514e-01 0.647 0.5530
## Y.14 2.526e-01 3.778e-01 0.669 0.5404
## Y.15 -2.912e-01 3.697e-01 -0.788 0.4750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4023000 on 4 degrees of freedom
## Multiple R-squared: 0.9744, Adjusted R-squared: 0.7956
## F-statistic: 5.448 on 28 and 4 DF, p-value: 0.05499
Alt text
#Nilai AIC dan BIC ARDL2
AIC(model.ardl2)## [1] 1087.712BIC(model.ardl2)## [1] 1132.608#Forecasting ARDL2
(fore.ardl2 <- forecast(model = model.ardl2, x=test$Sales_quantity, h=16))## $forecasts
## [1] 40483034 28777109 46218740 32747217 39517176 50694727 41702983 18899763
## [9] 16834881 38032246 25387412 75003288 65338291 29716280 20624133 12285394
##
## $call
## forecast.ardlDlm(model = model.ardl2, x = test$Sales_quantity,
## h = 16)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"Output di atas merupakan nilai forecasting model autoregressive dengan p=12 dan q=15 untuk meramalkan 16 bulan ke depan.
#MAPE testing ARDL2
mape.ardl2 <- MAPE(fore.ardl2$forecasts, test$Revenue) #data testing
#MAPE data training ARDL2
mape_train2 <- GoF(model.ardl2)["MAPE"]
c("MAPE_testing" = mape.ardl2, "MAPE_training" = mape_train2)## $MAPE_testing
## [1] 0.2796527
##
## $MAPE_training.MAPE
## [1] 0.04035136Perbandingan Akurasi (MAPE) Model Forecasting
#Perbandingan MAPE
akurasi <- matrix(c(mape.koyck, mape.dlm, mape.dlm2, mape.ardl, mape.ardl2))
row.names(akurasi)<- c("Koyck","DLM lag=2","DLM lag=4 (optimum)","ARDL lag x=1 & lag y=1","ARDL lag x=12 & lag y=15 (optimum)")
colnames(akurasi) <- c("MAPE")
knitr::kable(akurasi, align = "c")| MAPE | |
|---|---|
| Koyck | 0.1321049 |
| DLM lag=2 | 0.1247095 |
| DLM lag=4 (optimum) | 0.1236090 |
| ARDL lag x=1 & lag y=1 | 0.1162367 |
| ARDL lag x=12 & lag y=15 (optimum) | 0.2796527 |
Berdasarkan tabel perbandingan Mean Absolute Percentage Error (MAPE) di atas, model autoregressive dengan lag x =1 dan lag y = 1 merupakan model terbaik dalam forecasting model data penjualan-pendapatan karena model ini memiliki nilai MAPE terkecil dibandingkan oleh model lainnya.
Plot Perbandingan Model Forecasting
par(mfrow=c(1,1))
plot(test$Sales_quantity, test$Revenue, type="b", col="black",pch=20, ylim=c(12000000,75000000), main="Perbandingan Model Forecasting vs Aktual")
points(test$Sales_quantity, fore.koyck$forecasts,col="red",pch=20)
lines(test$Sales_quantity, fore.koyck$forecasts,col="red")
points(test$Sales_quantity, fore.dlm$forecasts,col="blue",pch=20)
lines(test$Sales_quantity, fore.dlm$forecasts,col="blue")
points(test$Sales_quantity, fore.dlm2$forecasts,col="orange",pch=20)
lines(test$Sales_quantity, fore.dlm2$forecasts,col="orange")
points(test$Sales_quantity, fore.ardl$forecasts,col="green",pch=20)
lines(test$Sales_quantity, fore.ardl$forecasts,col="green")
points(test$Sales_quantity, fore.ardl2$forecasts,col="purple",pch=20)
lines(test$Sales_quantity, fore.ardl2$forecasts,col="purple")
legend("bottomright",c("aktual", "koyck", "DLM 1", "DLM 2", "ARDL 1", "ARDL 2"), lty=1, col=c("black","red","blue","orange", "green", "purple"), cex=0.8)
Analisis eksploratif dengan plot di atas, terlihat kurang efektif karena keterbatasan penglihatan mata manusia dan bersifat subjektif. Menurut penulis, terlihat bahwa plot model distribusi lag dengan lag = 4 (optimum) (yang ditandai dengan warna orange) paling mendekati data aktual (yang ditandai dengan warna hitam) dibandingkan oleh metode lainnya. Hal ini menandakan bahwa model distribusi lag dengan lag = 4 (optimum) lebih baik dalam peramalan data yang digunakan.Namun, perbandingan akurasi model dengan MAPE cenderung lebih objektif.
Uji Formal Asumsi Klasik Model
Uji Non Autokorelasi
Hipotesis yang diuji:
H0: E[𝜀𝑖,𝜀𝑗]=0, ∀𝑖≠𝑗 (sisaan saling bebas atau tidak ada
autokorelasi atau 𝜌 = 0)
H1: E[𝜀𝑖, 𝜀𝑗]≠0, ∀𝑖≠𝑗 (sisaan tidak saling bebas atau ada
autokorelasi atau 𝜌 ≠ 0)
dengan menggunakan α=5%.
bgtest(model.dlm2$model)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model.dlm2$model
## LM test = 0.84226, df = 1, p-value = 0.3588bgtest(model.ardl$model)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model.ardl$model
## LM test = 2.9791, df = 1, p-value = 0.08434Berdasarkan uji formal: Uji Breusch-Godfrey pada model autoregressive dengan lag x = 1 & lag y = 1 dan model distribusi lag dengan lag = 4 (optimum) menunjukkan p-value > α, maka tak tolak H0. Artinya, cukup bukti untuk menyatakan bahwa sisaan saling bebas atau tidak ada autokorelasi pada taraf nyata 5%.
Uji Homoskedastisitas
Hipotesis yang diuji:
H0: Ragam sisaan homogen
H1: Ragam sisaan tidak homogen
dengan menggunakan α=5%.
#Breush-Pagan Test
bptest(model.dlm2$model)##
## studentized Breusch-Pagan test
##
## data: model.dlm2$model
## BP = 4.6253, df = 5, p-value = 0.4633bptest(model.ardl$model)##
## studentized Breusch-Pagan test
##
## data: model.ardl$model
## BP = 1.063, df = 3, p-value = 0.786Berdasarkan uji formal: Uji Breush-Pagan pada model autoregressive dengan lag x = 1 & lag y = 1 dan model distribusi lag dengan lag = 4 (optimum) menunjukkan p-value > α, maka tak tolak H0. Artinya, cukup bukti untuk menyatakan bahwa ragam sisaan homogen pada taraf nyata 5%.
Uji Normalitas
Hipotesis yang diuji:
H0: Sisaan menyebar normal
H1: Sisaan tidak menyebar normal
dengan menggunakan α=5%.
#Shapiro Wilk Test
shapiro.test(residuals(model.dlm2))## 1 2 3 4 5 6 7
## -2827866.4 -3480751.9 -5708815.6 -448650.5 -2043472.9 -7708223.2 -3555979.6
## 8 9 10 11 12 13 14
## -3718220.4 -732574.9 994817.0 -1393306.1 4288579.3 -1612492.7 1361591.1
## 15 16 17 18 19 20 21
## 6199045.6 320675.6 11261591.8 -1384326.0 -645024.5 3159320.1 1979874.0
## 22 23 24 25 26 27 28
## 3531517.8 8248825.8 -3955451.8 1899683.7 672015.6 4645402.4 -2142179.0
## 29 30 31 32 33 34 35
## 1739698.4 -343095.1 -2463564.7 121652.8 4838715.8 3558124.7 -2204658.5
## 36 37 38 39 40 41 42
## -4206982.2 -1591646.7 -7095073.9 2032893.9 -230992.8 -653917.1 3025537.1
## 43 44
## 352592.2 -4084888.4##
## Shapiro-Wilk normality test
##
## data: residuals(model.dlm2)
## W = 0.97801, p-value = 0.5558shapiro.test(residuals(model.ardl))## Time Series:
## Start = 2
## End = 48
## Frequency = 1
## 2 3 4 5 6 7
## -1701196.04 -1218019.58 -5349421.51 -219749.56 -2006276.33 -6294767.92
## 8 9 10 11 12 13
## 684645.69 -2577331.42 -6840123.72 -3064662.70 -3726204.58 2357483.61
## 14 15 16 17 18 19
## -480972.22 1443148.40 3946761.86 -1078361.08 2631612.70 4181457.29
## 20 21 22 23 24 25
## -1163595.68 10742872.37 -4945891.49 -2088871.69 2009380.55 3835210.92
## 26 27 28 29 30 31
## 1909327.89 11108665.83 -7583581.20 5544224.53 1743747.13 3591113.83
## 32 33 34 35 36 37
## -3047086.87 2563719.80 -1156529.37 -4491872.84 -710306.50 6853487.83
## 38 39 40 41 42 43
## 1829288.55 309909.13 -5072826.09 2151518.91 -3452917.96 4015290.67
## 44 45 46 47 48
## -45934.19 98407.40 1634216.39 -2244471.25 -4624519.50##
## Shapiro-Wilk normality test
##
## data: residuals(model.ardl)
## W = 0.97031, p-value = 0.2723Berdasarkan uji formal: Uji Shapiro Wilk pada model autoregressive dengan lag x =1 dan lag y = 1 dan model distribusi lag dengan lag = 4 (optimum) menunjukkan p-value > α, maka tak tolak H0. Artinya, cukup bukti untuk menyatakan bahwa sisaan menyebar normal pada taraf nyata 5%.
Simpulan
Berdasarkan analisis regresi deret waktu yang sudah dilakukan,
diperoleh beberapa simpulan sebagai berikut:
1. Analisis regresi untuk menangani autokorelasi jika data dipengaruhi
oleh data pada periode sebelumnya (peubah lag) dengan model
Koyck, model distribusi lag, dan model
autoregressive.
2. Ada perbedaan hasil forecasting data yang
digunakan berdasarkan perbandingan secara eksploratif (plot)
dan nilai akurasi model secara objektif (MAPE). Berdasarkan
plot, model distribusi lag dengan lag
= 4 (optimum) menjadi model terbaik, sedangkan model
autoregressive dengan lag x =1 dan lag y = 1 merupakan
model terbaik berdasarkan perbandingan nilai
MAPE yang terkecil.
3. Kedua model memenuhi uji asumsi non autokorelasi, homoskedastisitas,
dan normalitas yang diindikasikan dengan p-value > α.
Daftar Pustaka
Anderson DR, Burnham KP, White GC. 1998. Comparison of akaike information criterion and consistent akaike information criterion for model selection and statistical inference from capture-recapture studies. Journal of Applied Statistics. 25(2): 263-282.
Gujarati DN. 2004. Basic Econometrics. Fourth Edition. New York: The McGrawHill.
Handoko A. 2013. Model distribusi lag dan distribusi autoregressive [skripsi]: Universitas Lambung Mangkurat Banjarmasin.
Jatiningrum N. Model dinamis: autoregressive dan distribusi lag [skripsi]: Universitas Negeri Yogyakarta.
Nasution DA, Khotimah HH, Chamidah N. 2019. Perbandingan normalisasi data untuk klasifikasi wine menggunakan algoritma k-nn. CESS (Journal of Computer Engineering, System and Science). 4(1):78-82.
Pawestri V, Setiawan A, Linawati L. 2019. Pemodelan data penjualan mobil menggunakan model autoregressive moving average berdasarkan metode bayesian. Jurnal Sains dan dan Eduka Sains. 2(1): 26-35.
Qudratullah, Farhan M. 2007. Bayesian information criterion (BIC) dalam pemilihan model terbaik feed forward neural network (FFNN) : Studi kasus data posisi dana simpanan tabungan bank Umum dan BPR di provinsi D.I. Yogyakarta [skripsi]: Universitas Gajah Mada.
Robial SM. 2018. Perbandingan model statistik pada analisis peramalan time series (studi kasus: PT. Telekomunikasi Indonesia, Tbk Kandatel Sukabumi). Jurnal Ilmiah SANTIKA. 8(2): 1-17.