Tugas Individu 2 MPDW: Regresi Peubah Lag
Regresi Peubah Lag
Analisis regresi merupakan salah satu teknik analisis data dalam statistika untuk menentukan hubungan sebab-akibat antara variabel takbebas (y) dengan satu atau lebih variabel bebas (x). Analisis regresi dengan serangkaian pengamatan terhadap suatu peristiwa diambil dari waktu ke waktu merupakan analisis regresi deret waktu. Analisis regresi deret waktu dilakukan untuk melihat atau meramalkan kondisi pada masa depan (Gujarati 2006). Model regresi menggunakan data deret waktu tidak hanya menggunakan pengaruh perubahan variabel bebas (x) terhadap variabel takbebas (y) dalam kurun waktu yang sama dan selama periode pengamatan yang sama, namun juga menggunakan periode waktu sebelumnya. Waktu yang diperlukan bagi variabel bebas (x) dalam memengaruhi variabel takbebas (y) disebut beda kala atau lag (Sarwoko 2005).
Pada penelitian ini, akan dilakukan analisis regresi pada data penjualan mingguan Walmart dengan data Indeks Harga Konsumen. Analisis regresi dilakukan dengan menggunakan metode model koyck, model distributed lag, dan model autoregressive distributed lag.
Data Preparation
Memanggil Package
library(dLagM) ## Warning: package 'dLagM' was built under R version 4.1.3## Loading required package: nardl## Warning: package 'nardl' was built under R version 4.1.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Loading required package: dynlm## Warning: package 'dynlm' was built under R version 4.1.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.1.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(dynlm)
library(MLmetrics) ## Warning: package 'MLmetrics' was built under R version 4.1.3##
## Attaching package: 'MLmetrics'## The following object is masked from 'package:dLagM':
##
## MAPE## The following object is masked from 'package:base':
##
## Recalllibrary(lmtest)## Warning: package 'lmtest' was built under R version 4.1.3library(car)## Warning: package 'car' was built under R version 4.1.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.1.2library(caTools)## Warning: package 'caTools' was built under R version 4.1.3library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.1.3## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --## v ggplot2 3.3.6 v purrr 0.3.4
## v tibble 3.1.8 v dplyr 1.0.10
## v tidyr 1.2.1 v stringr 1.4.1
## v readr 2.1.2 v forcats 0.5.2## Warning: package 'ggplot2' was built under R version 4.1.3## Warning: package 'tibble' was built under R version 4.1.3## Warning: package 'tidyr' was built under R version 4.1.3## Warning: package 'readr' was built under R version 4.1.2## Warning: package 'purrr' was built under R version 4.1.2## Warning: package 'dplyr' was built under R version 4.1.3## Warning: package 'stringr' was built under R version 4.1.3## Warning: package 'forcats' was built under R version 4.1.3## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::recode() masks car::recode()
## x purrr::some() masks car::some()library(hrbrthemes)## Warning: package 'hrbrthemes' was built under R version 4.1.2## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
## Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
## if Arial Narrow is not on your system, please see https://bit.ly/arialnarrowlibrary(knitr)## Warning: package 'knitr' was built under R version 4.1.3Menginput Data
Data yang digunakan adalah data Walmart Sales Forecast dari Kaggle dengan 2 peubah uji yaitu Weekly Sales (Penjualan mingguan Walmart) sebagai peubah respon dan CPI (Indeks Harga Konsumen) sebagai peubah penjelas dengan rentang waktu mingguan dari 31 Desember 2010 hingga 9 Desember 2011 dengan jumlah amatan sebanyak 50 amatan.
#Import data
data1<- read.csv("C:\\Users\\tiasa\\Downloads\\Data1.csv")
data2<- read.csv("C:\\Users\\tiasa\\Downloads\\Data2.csv")
t<-data1[192:241,]$Date
Xt<-data2[231:280,]$CPI
Yt<-data1[192:241,]$Weekly_Sales
df <- as.data.frame(cbind(t, Xt, Yt))
df$t <- as.Date(df$t)
df$Xt <- as.numeric(df$Xt)
df$Yt <- as.numeric(df$Yt)
str(df)## 'data.frame': 50 obs. of 3 variables:
## $ t : Date, format: "2011-01-07" "2011-01-14" ...
## $ Xt: num 211 211 211 212 212 ...
## $ Yt: num 43202 39857 41597 41651 46829 ...head(df)## t Xt Yt
## 1 2011-01-07 211.0649 43202.29
## 2 2011-01-14 211.1177 39857.40
## 3 2011-01-21 211.4865 41596.99
## 4 2011-01-28 211.8553 41651.01
## 5 2011-02-04 212.2241 46829.12
## 6 2011-02-11 212.5929 47595.35Splitting Data
Data dibagi menjadi data train dan data test dengan proporsi masing-masing 80% dan 20% dari jumlah seluruh amatan.
#Split data
train<-df[1:40,]
test<-df[41:50,]
#Data time series
train.ts<-ts(train)
test.ts<-ts(test)
df.ts<-ts(df)Analisis Eksplorasi Data
Plot Time Series
data.ts1<-ts(Yt)
plot(data.ts1, main = "Time Series Plot of Weekly Sales", xlab = "Period", ylab="Weekly Sales")
points(data.ts1)
Berdasarkan plot, terlihat bahwa peubah Weekly Sales atau penjualan mingguan Walmart memiliki pola data yang stasioner.
data.ts<-ts(Xt)
plot(data.ts, main = "Time Series Plot of Customer Price Index", xlab = "Period", ylab="Customer Price Index")
points(data.ts)
Berdasarkan plot, terlihat bahwa peubah CPI atau Indeks Harga Konsumen memiliki pola data yang tren positif.
Korelasi dan Scatter Plot Peubah X dan Y
Korelasi
cor(Xt, Yt)## [1] 0.1675282Scatter Plot
plot(Xt, Yt, pch = 20, col = "blue", main = "Scatter Plot CPI & Weekly Sales")
Berdasarkan korelasi dan scatter plot antara peubah CPI dan Weekly sales, terlihat bahwa terdapat hubungan positif yang lemah dengan korelasi sebesar 0.1675282.
Model Regresi Awal
model1 <- lm(df$Yt~df$Xt, data = df)
summary(model1)##
## Call:
## lm(formula = df$Yt ~ df$Xt, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4737.8 -2305.2 -477.7 1985.2 7848.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8318.6 45804.4 -0.182 0.857
## df$Xt 250.6 212.9 1.177 0.245
##
## Residual standard error: 2783 on 48 degrees of freedom
## Multiple R-squared: 0.02807, Adjusted R-squared: 0.007817
## F-statistic: 1.386 on 1 and 48 DF, p-value: 0.2449Diperoleh model regresi linier data deret waktu yaitu: \(y = -8318.6 + 250.6x\)
Berdasarkan persamaan regresi, nilai dugaan rataan Y bertambah sebesar 250.6 jika X bertambah satu satuan.
Berdasarkan uji t yang dilakukan untuk menguji signifikansi peubah penjelas terhdap variabel dependen, diperoleh bahwa tidak ada peubah yang berpengaruh signifikan pada taraf nyata 5%.
Nilai R-Squared pada model regresi linier deret waktu yaitu 0.007817, artinya keragaman peubah respon yang mampu dijelaskan oleh peubah penjelas adalah sebesar 0.007817.
Uji Autokorelasi
Uji autokorelasi merupakan uji yang dilakukan untuk dapat melihat apakah terjadi korelasi di antara suatu periode dengan periode-periode sebelumnya. Uji statistik yang digunakan untuk pendeteksian autokorelasi uji Breusedch grodfrey.
Hipotesis:
H0: tidak ada autokorelasi
H1: ada autokorelasi
bgtest(model1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model1
## LM test = 4.4659, df = 1, p-value = 0.03458Diperoleh p-value < 5% : Tolak H0, sehingga cukup bukti untuk menyatakan bahwa terdapat autokorelasi (positif) pada taraf 5%. Karena terdeteksi adanya autokorelasi, akan dilakukan penanganan dengan menambahkan lag tertentu. Metode yang digunakan adalah Model Koyck, Distributed Lag, dan Autoregressive Distributed Lag.
Model KOYCK
Model Koyck merupakan suatu model yang mengasumsikan bahwa koefisien dari variabel yang mengalami keterlambatan dapat menurun secara geometris (Gujarati 2006).
model.koyck = dLagM::koyckDlm(x = train$Xt, y = train$Yt)
summary(model.koyck)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -4338.8 -2224.1 -224.7 1833.0 7422.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.061e+04 7.260e+04 -0.284 0.7782
## Y.1 2.793e-01 1.637e-01 1.706 0.0966 .
## X.t 2.498e+02 3.429e+02 0.728 0.4711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2700 on 36 degrees of freedom
## Multiple R-Squared: 0.1056, Adjusted R-squared: 0.05588
## Wald test: 2.012 on 2 and 36 DF, p-value: 0.1485
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: -28593.59 249.7842 0.2793467AIC(model.koyck)## [1] 731.8415BIC(model.koyck)## [1] 738.4957Forecasting Model Koyck
(fore.koyck <- forecast(model = model.koyck, x=test$Xt, h=10))## $forecasts
## [1] 47532.33 46916.33 46784.26 46787.38 46828.26 46895.01 46975.10 47058.92
## [9] 47143.78 47221.63
##
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 10)
##
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"Akurasi Model Koyck
#mape data testing
mape.koyck <- MAPE(fore.koyck$forecasts, test$Yt)
#akurasi data training
mape_train <- dLagM::GoF(model.koyck)["MAPE"]
c("MAPE_testing" = mape.koyck, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.0650684
##
## $MAPE_training.MAPE
## [1] 0.04656572Distributed Lag Model (DLM)
Model distributed lag disebut sebagai model dinamis karena efek perubahan satu unit dalam nilai variabel bebas terdistribusi pada sejumlah periode waktu (Gujarati 2006).
lagop<-finiteDLMauto(formula=Yt ~ Xt,
data = data.frame(train),model.type = "dlm",error.type = "AIC");lagop## q - k MASE AIC BIC GMRAE MBRAE R.Adj.Sq Ljung-Box
## 10 10 0.65378 566.3786 584.5942 0.70362 0.70537 -0.27964 0.6270995Diperolah lag optimum untuk peubah Weekly Sales adalah 10 hari sebelumnya.
model.dlm = dLagM::dlm(x = train$Xt,y = train$Yt , q = 10)
summary(model.dlm)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2901.7 -1122.1 -443.3 1057.4 4235.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -198439.0 206993.6 -0.959 0.350
## x.t 61.2 5771.3 0.011 0.992
## x.1 1952.8 12994.3 0.150 0.882
## x.2 -3805.8 14664.3 -0.260 0.798
## x.3 5919.9 14905.5 0.397 0.696
## x.4 2245.9 15095.4 0.149 0.883
## x.5 -5499.2 15249.6 -0.361 0.723
## x.6 -5225.8 15136.3 -0.345 0.734
## x.7 9801.5 15052.2 0.651 0.523
## x.8 -871.0 14671.3 -0.059 0.953
## x.9 -7023.4 11996.7 -0.585 0.566
## x.10 3576.4 5068.7 0.706 0.489
##
## Residual standard error: 2547 on 18 degrees of freedom
## Multiple R-squared: 0.2057, Adjusted R-squared: -0.2796
## F-statistic: 0.4239 on 11 and 18 DF, p-value: 0.9256
##
## AIC and BIC values for the model:
## AIC BIC
## 1 566.3786 584.5942Diperoleh tidak terdapat peubah yang berpengaruh terhadap model dalam tingkat kepercayaan 95%. Berdasarkan R-Squared, model mengindikasikan performa yang kurang baik dengan nilai sebesar 0.2796.
AIC(model.dlm)## [1] 566.3786BIC(model.dlm)## [1] 584.5942Forecasting DLM
(fore.dlm <- forecast(model = model.dlm, x=test$Xt, h=10)) #meramalkan 12 periode ke depan## $forecasts
## [1] 48160.00 49250.94 48698.56 48894.70 49564.72 49094.15 49023.44 49853.33
## [9] 49515.78 49763.11
##
## $call
## forecast.dlm(model = model.dlm, x = test$Xt, h = 10)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"Akurasi DLM
#mape data testing
mape.dlm <- MAPE(fore.dlm$forecasts, test$Yt)
#akurasi data training
mape_train <- GoF(model.dlm)["MAPE"]
c("MAPE_testing" = mape.dlm, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.09150085
##
## $MAPE_training.MAPE
## [1] 0.035039Autoregressve Distributed Lag Model (ARDLM)
Model regresi dengan memasukkan nilai variabel yang menjelaskan nilai masa kini atau nilai masa lalu dari variabel bebas (X) sebagai tambahan pada model yang memasukkan lag dari variabel takbebas (Y) disebut autoregressive distributed lag (ARDL).
pqop<-ardlBoundOrders(data = data.frame(df) , formula = Yt ~ Xt )
c(p=pqop$p$Xt, q=pqop$q)## p q
## 14 15Diperoleh lag optimum untuk peubah Xt adalah 14 hari sebelumnya, dan lag optimum untuk peubah Yt adalah 15 hari sebelumnya.
model.ardl = ardlDlm(x = train$Xt, y = train$Yt, p = 14 , q = 9) #seharusnya nilai q optimum 15, namun karena menghasilkan nilai NaN pada output, q yang optimum pada model bernilai 9
summary(model.ardl)##
## Time series regression with "ts" data:
## Start = 15, End = 40
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## 15 16 17 18 19 20 21 22
## -309.050 1014.842 -807.012 -278.759 552.627 -366.918 92.173 240.044
## 23 24 25 26 27 28 29 30
## -151.681 -94.538 -20.120 176.476 8.622 -203.076 -79.241 233.208
## 31 32 33 34 35 36 37 38
## 89.059 -155.741 -246.013 184.908 185.278 67.239 210.621 -706.694
## 39 40
## 531.407 -167.661
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.745e+05 5.576e+05 -0.492 0.709
## X.t -2.153e+03 7.476e+03 -0.288 0.821
## X.1 8.040e+02 1.550e+04 0.052 0.967
## X.2 6.702e+03 1.769e+04 0.379 0.769
## X.3 -8.296e+03 1.800e+04 -0.461 0.725
## X.4 7.618e+03 1.790e+04 0.426 0.744
## X.5 3.343e+03 1.745e+04 0.192 0.880
## X.6 -1.221e+04 1.536e+04 -0.795 0.572
## X.7 1.696e+04 1.589e+04 1.067 0.479
## X.8 -1.869e+04 1.569e+04 -1.191 0.445
## X.9 1.391e+04 1.477e+04 0.942 0.519
## X.10 2.370e+03 1.394e+04 0.170 0.893
## X.11 -1.709e+04 1.496e+04 -1.143 0.458
## X.12 6.508e+03 2.134e+04 0.305 0.812
## X.13 -7.316e+03 1.569e+04 -0.466 0.722
## X.14 9.303e+03 6.497e+03 1.432 0.388
## Y.1 -5.030e-01 5.776e-01 -0.871 0.544
## Y.2 -4.523e-01 4.825e-01 -0.938 0.521
## Y.3 -4.634e-01 4.185e-01 -1.107 0.468
## Y.4 -2.748e-01 4.539e-01 -0.605 0.653
## Y.5 -4.995e-01 5.155e-01 -0.969 0.510
## Y.6 -1.763e-01 3.848e-01 -0.458 0.727
## Y.7 6.261e-03 4.225e-01 0.015 0.991
## Y.8 1.978e-01 4.191e-01 0.472 0.719
## Y.9 8.569e-01 3.799e-01 2.256 0.266
##
## Residual standard error: 1883 on 1 degrees of freedom
## Multiple R-squared: 0.9737, Adjusted R-squared: 0.3415
## F-statistic: 1.54 on 24 and 1 DF, p-value: 0.5717Tidak ada peubah yang berpengaruh terhadap model dalam tingkat kepercayaan 95%. Model tersebut merupakan model yang kurang baik karena memiliki nilai Adjusted R-Squared sebesar 34,15%.
AIC(model.ardl)## [1] 433.1892BIC(model.ardl)## [1] 465.8997Forecasting ARDLM
#ramalan
(fore.ardl <- forecast(model = model.ardl, x=test$Xt, h=10))## $forecasts
## [1] 46697.03 49962.22 49452.46 47756.07 42233.90 44908.22 48147.13 52560.02
## [9] 50015.07 48202.69
##
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 10)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"Akurasi ARDLM
#akurasi testing
mape.ardl <- MAPE(fore.ardl$forecasts, test$Yt) #data testing
#akurasi data training
mape_train <- GoF(model.ardl)["MAPE"]
c("MAPE_testing" = mape.ardl, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.08337594
##
## $MAPE_training.MAPE
## [1] 0.006089235Perbandingan Metode Penanganan Autokorelasi
Perbandingan Secara Eksploratif
#plot
par(mfrow=c(1,1))
plot(test$Xt, test$Yt, type="b", col="black")
points(test$Xt, fore.koyck$forecasts,col="red")
points(test$Xt, fore.dlm$forecasts,col="blue")
points(test$Xt, fore.ardl$forecasts,col="green")
legend("topleft",c("aktual", "koyck","DLM", "autoregressive"), lty=1, col=c("black","red","blue","green"), cex=0.8)
Secara eksploratif, belum terlihat metode yang paling sesuai dan mendekati pola aktual.
Perbandingan Keakuratan Forecasting
akurasi <- matrix(c(mape.koyck,mape.dlm, mape.ardl))
row.names(akurasi)<- c("Koyck","DLM","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi## MAPE
## Koyck 0.06506840
## DLM 0.09150085
## Autoregressive 0.08337594Berdasarkan akurasi ramalan, diperoleh nilai error model Koyck yang terendah daripada model lainnya. Maka, dapat disimpulkan bahwa model Koyck merupakan model terbaik untuk modelling dan forecasting.
Uji Diagnostik Model
Uji Autokorelasi
Hipotesis:
H0: tidak ada autokorelasi
H1: ada autokorelasi
bgtest(model.koyck$model)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model.koyck$model
## LM test = 0.25269, df = 1, p-value = 0.6152Diperoleh p-value > 0,05, artinya dalam tingkat signifikansi 95%, terdapat cukup bukti untuk menyatakan bahwa tidak terdapat autokorelasi, artinya penanganan autokorelasi dengan menambahkan lag dari peubah respon berhasil.
Daftar Pustaka
Kaggle. 2013. Walmart Sales Forecast. https://www.kaggle.com/datasets/aslanahmedov/walmart-sales-forecast
Gujarati DN. 2006. Essentials of Econometrics. Third Edition. New York: The McGraw-Hill.
Sarwoko. 2005. Dasar-Dasar Ekonometrika. Yogyakarta: Penerbit Andi.