Praktikum 03 - Regresi dengan Peubah Lag
Angga Fathan Rofiqy (G1401211006)
I. Setting
Abaikan bagian ini. Lansung saja loncat ke
1. Tentang Data . Function biar gak perlu
ganti backslash (\) jadi slash
(/).
path <- function() gsub ( "\\\\", "/", readClipboard () )
#Copy path, Panggil function di console
#Copy r path, paste ke var yang diinginkanFunction biar gak perlu repot buat
install("") atau load() package.
# -=( Install & Load Package Function )=-
install_load <- function (package1, ...) {
# convert arguments to vector
packages <- c(package1, ...)
# start loop to determine if each package is installed
for(package in packages){
# if package is installed locally, load
if(package %in% rownames(installed.packages()))
do.call('library', list(package))
# if package is not installed locally, download, then load
else {
install.packages(package)
do.call("library", list(package))
}
}
}A. Tentang Data
Dataset ini dibangun dari dataset awal yang terdiri dari \(600.000\) data transaksi yang dikumpulkan dalam waktu 6 tahun (periode 2014-2019), yang mencatat tanggal dan waktu penjualan, nama merek obat farmasi, dan jumlah yang terjual. Data ini diekspor dari sistem Point-of-Sale di apotek individu. Kelompok obat yang dipilih dari dataset (57 obat) diklasifikasikan ke dalam kategori berikut dalam Sistem Klasifikasi Anatomi Terapeutik Kimia (ATC):
M01AB- Produk antiinflamasi dan antirheumatik, non-steroid, Turunan asam asetat dan substansi terkaitM01AE- Produk antiinflamasi dan antirheumatik, non-steroid, Turunan asam propionatN02BA- Analgesik dan antipiretik lainnya, Asam salisilat dan turunannyaN02BE/B- Analgesik dan antipiretik lainnya, Pirazolona dan AnilidaN05B- Obat psikoleptik, Obat ansiolitikN05C- Obat psikoleptik, Obat hipnotik dan sedatifR03- Obat untuk penyakit saluran napas obstruktifR06- Antihistamin untuk penggunaan sistemik
Data penjualan diambil ulang dalam periode per jam, per hari, per minggu, dan per bulan. Data sudah diproses sebelumnya, termasuk deteksi dan penanganan outlier serta pengisian data yang hilang. Data yang akan saya gunakan adalah data dalam periode per bulan.
Beberapa waktu yang lalu saya merasa sulit untuk tidur. Kebetulan ada
obat N05C yang merupakan obat tidur dalam dataset saya.
N05C merupakan obat penenang, sifat tenang ini yang mungkin
membantu seseorang untuk tidur. Ternyata ada obat tipe penenang lainnya,
yakni N05B.
Kedua obat ini memiliki kesamaan yakni termasuk dalam klasifikasi ATC N (Nervous System), yang berarti keduanya berhubungan dengan sistem saraf. Kedua obat ini termasuk dalam kelompok “Psycholeptics drugs,” yang berarti keduanya dapat memengaruhi sistem saraf dan suasana hati pasien.
Namun ada beberapa perbedaan, diantaranya: N05B lebih
fokus pada obat penenang dan pengobatan
gangguan kecemasan. Sementara N05C lebih
fokus pada obat tidur dan pengobatan gangguan
tidur. Contoh obat-obatan yang termasuk dalam kategori ini
berbeda, dengan N05B memiliki contoh obat-obatan penenang
seperti benzodiazepin, sementara N05C memiliki
obat tidur seperti zolpidem.
Tujuan Praktikum kali ini adalah membuat Model Regresi dengan Peubah Lag terbaik. Lag secara bahasa sendiri memiliki arti keterlambatan. Bisa dikatakan Lag mengacu pada nilai suatu variabel pada waktu sebelumnya dalam rangkaian data waktu.
Kali ini saya akan melihat “Apakah penjualan
obat penenang N05B pada
periode waktu saat
inidipengaruhi oleh penjualan obat tidur
N05C pada periode waktu
sebelumnya?” Atau “Apakah penjualan
obat tidur N05Cperiodesebelumnya memiliki
pengaruh signifikan terhadap penjualan
obat penenang N05B pada
periode waktu saat ini
?.
Import Data
install_load('rio')
raw.data <- import("https://raw.githubusercontent.com/Zen-Rofiqy/STA1341-MPDW/main/Pertemuan%203/salesmonthly.csv")Pembagian Data
#SPLIT DATA
data <- raw.data
train <- data[1:round(nrow(data) *80/100),] #80% data
test <- data[(nrow(train)+1) : nrow(data),] #20% data#data time series
train.ts <- ts(train)
test.ts <- ts(test)
data.ts <- ts(data)1. Model Koyck
Model Koyck didasarkan pada asumsi bahwa semakin jauh jarak lag peubah independen dari periode sekarang maka semakin kecil pengaruh peubah lag terhadap peubah dependen.
Koyck mengusulkan suatu metode untuk menduga model dinamis distributed lag dengan mengasumsikan bahwa semua koefisien \(\beta\) mempunyai tanda sama.
Model kyock merupakan jenis paling umum dari model infinite distributed lag dan juga dikenal sebagai geometric lag
\[ y_t=a(1-\lambda)+\beta_0X_t+\beta_1Z_t+\lambda Y_{t-1}+V_t \]
dengan \[V_t=u_t-\lambda u_{t-1}\]
Pemodelan
Pemodelan model Koyck dengan R dapat menggunakan
dLagM::koyckDlm() . Fungsi umum dari koyckDlm
adalah sebagai berikut.
koyckDlm(x , y , intercept)Fungsi koyckDlm() akan menerapkan model lag
terdistribusi dengan transformasi Koyck satu prediktor. Nilai
x dan y tidak perlu sebagai objek time
series (ts). intercept dapat dibuat
TRUE untuk memasukkan intersep ke dalam model.
install_load('dLagM')## Warning: package 'dLagM' was built under R version 4.2.3## Loading required package: nardl## Warning: package 'nardl' was built under R version 4.2.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Loading required package: dynlm## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric#MODEL KOYCK
model.koyck <- koyckDlm(x = train$N05C, y = train$N05B)
summary(model.koyck)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -418.505 -56.352 -3.085 56.335 365.217
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 146.5938 150.3615 0.975 0.3341
## Y.1 1.0619 0.6181 1.718 0.0917 .
## X.t -9.8906 17.9684 -0.550 0.5844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 118.9 on 52 degrees of freedom
## Multiple R-Squared: -0.6422, Adjusted R-squared: -0.7054
## Wald test: 8.872 on 2 and 52 DF, p-value: 0.0004841
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: -2367.792 -9.890557 1.061912AIC(model.koyck)## [1] 686.6325BIC(model.koyck)## [1] 694.6618Dari hasil tersebut, didapat bahwa peubah \(x_t\) dan \(y_{t-1}\) memiliki nilai \(P.Value(0.0004841)<0.05\). Hal ini menunjukkan bahwa peubah \(x_t\) dan \(y_{t-1}\) berpengaruh signifikan terhadap \(y\). Adapun model keseluruhannya adalah sebagai berikut
\[ \hat{Y_t}=146.5938+ 1.0619X_t-9.8906 Y_{t-1} \]
Peramalan dan Akurasi
Berikut adalah hasil peramalan y untuk \(14\) periode kedepan menggunakan model
koyck
install_load('MLmetrics')##
## Attaching package: 'MLmetrics'## The following object is masked from 'package:dLagM':
##
## MAPE## The following object is masked from 'package:base':
##
## Recallfore.koyck <- forecast(model = model.koyck, x=test$N05C, h=14)
fore.koyck## $forecasts
## [1] 296.57072 214.26173 156.52854 45.76820 -32.28726 -6.37913
## [7] -48.10089 -122.07737 -240.19609 -355.73717 -428.97880 -566.09828
## [13] -612.80147 -573.38112
##
## $call
## forecast.koyckDlm(model = model.koyck, x = test$N05C, h = 14)
##
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"mape.koyck <- MAPE(fore.koyck$forecasts, test$N05C)
#akurasi data training
GoF(model.koyck)## n MAE MPE MAPE sMAPE MASE MSE MRAE
## model.koyck 55 80.96454 -7.627032 7.892919 0.3323228 1.823982 13371.69 3.485585
## GMRAE
## model.koyck 1.824132Pada perhitungan keakuratan model menggunakan metode
Koyck didapatkan nilai MAPE \(7.89\%\). Nilai akurasi model ini
kurang dari \(10\%\) sehingga
dapat dikategorikan sangat baik.
2. Regression with Distributed Lag
Pemodelan model Regression with Distributed Lag dengan R
dapat menggunakan dLagM::dlm() . Fungsi umum dari
dlm adalah sebagai berikut.
dlm(formula , data , x , y , q , remove )Fungsi dlm() akan menerapkan model lag terdistribusi
dengan satu atau lebih prediktor. Nilai x dan
y tidak perlu sebagai objek time series
(ts). \(q\) adalah integer
yang mewakili panjang lag yang terbatas.
Pemodelan (Lag=2)
model.dlm <- dlm(x = train$N05C, y = train$N05B , q = 2)
summary(model.dlm)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.311 -45.216 0.153 35.738 172.852
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 106.58449 29.86292 3.569 0.000802 ***
## x.t 6.86841 1.30848 5.249 3.13e-06 ***
## x.1 2.65332 1.32657 2.000 0.050931 .
## x.2 -0.03826 1.12790 -0.034 0.973071
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 70.1 on 50 degrees of freedom
## Multiple R-squared: 0.4427, Adjusted R-squared: 0.4093
## F-statistic: 13.24 on 3 and 50 DF, p-value: 1.749e-06
##
## AIC and BIC values for the model:
## AIC BIC
## 1 618.0751 628.02AIC(model.dlm)## [1] 618.0751BIC(model.dlm)## [1] 628.02Dari hasil diatas, didapat bahwa \(P.value\) dari intercept dan \(x_{t-1}(1.749\times 10^{-6})<0.05\). Hal ini menunjukkan bahwa intercept dan \(x_{t-1}\) berpengaruh signifikan terhadap \(y\). Adapun model keseluruhan yang terbentuk adalah sebagai berikut
\[ \hat{Y_t}=106.58449 + 6.86841X_t + 2.65332X_{t-1} - 0.03826X_{t-2} \]
Peramalan dan Akurasi
Berikut merupakan hasil peramalan \(y\) untuk \(14\) periode kedepan
fore.dlm <- forecast(model = model.dlm, x=test$N05C, h=14)
fore.dlm## $forecasts
## [1] 241.3449 309.3695 323.5635 349.4482 335.3560 248.9988 268.0442 307.6436
## [9] 342.8094 346.4395 309.2910 337.2731 284.7003 196.1217
##
## $call
## forecast.dlm(model = model.dlm, x = test$N05C, h = 14)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"mape.dlm <- MAPE(fore.dlm$forecasts, test$N05C)
#akurasi data training
GoF(model.dlm)## n MAE MPE MAPE sMAPE MASE MSE MRAE
## model.dlm 54 51.6059 -2.854632 3.00373 0.2119 1.198559 4549.502 3.57092
## GMRAE
## model.dlm 1.143447Pada perhitungan keakuratan model menggunakan metode
Regression with Distributed Lag didapatkan nilai
MAPE \(3.0037\%\). Nilai
akurasi model ini kurang dari \(10\%\) sehingga dapat dikategorikan
sangat baik.
Lag Optimum
#penentuan lag optimum
finiteDLMauto(formula = N05B ~ N05C,
data = data.frame(train), q.min = 1, q.max = 6,
model.type = "dlm", error.type = "AIC", trace = FALSE)## q - k MASE AIC BIC GMRAE MBRAE R.Adj.Sq Ljung-Box
## 6 6 1.10173 563.8085 581.0168 1.39547 0.31176 0.56603 9.088995e-05Berdasarkan output tersebut, lag optimum didapatkan ketika
lag=\(6\). Selanjutnya
dilakukan pemodelan untuk lag=\(6\).
#model dlm dengan lag optimum
model.dlm2 <- dlm(x = train$N05C,y = train$N05B , q = 6)
summary(model.dlm2)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -137.341 -45.021 -6.279 41.848 123.681
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.0039 36.4634 -0.165 0.8700
## x.t 6.7690 1.1996 5.643 1.29e-06 ***
## x.1 2.4839 1.2163 2.042 0.0474 *
## x.2 0.3845 1.1924 0.322 0.7487
## x.3 2.8511 1.1958 2.384 0.0217 *
## x.4 2.1521 1.2057 1.785 0.0815 .
## x.5 0.8164 1.2037 0.678 0.5014
## x.6 1.1695 1.1205 1.044 0.3026
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 61.95 on 42 degrees of freedom
## Multiple R-squared: 0.628, Adjusted R-squared: 0.566
## F-statistic: 10.13 on 7 and 42 DF, p-value: 2.374e-07
##
## AIC and BIC values for the model:
## AIC BIC
## 1 563.8085 581.0167AIC(model.dlm2)## [1] 563.8085BIC(model.dlm2)## [1] 581.0167Dari hasil tersebut terdapat beberapa peubah yang berpengaruh signifikan terhadap taraf nyata \(5\%\) yaitu \(x_t\) , \(x_{t-1}\) , \(x_{t-3}\). Adapun keseluruhan model yang terbentuk adalah
\[ \hat{Y_t}=-6.0039 + 6.769 X_t + 2.4839 X_{t-1} + 0.3845 X_{t-2} + 2.8511 X_{t-3} + 2.1521 X_{t-4} + 0.8164 X_{t-5} + 1.1695 X_{t-6} \]
Adapun hasil peramalan \(14\) periode kedepan menggunakan model tersebut adalah sebagai berikut
#peramalan dan akurasi
fore.dlm2 <- forecast(model = model.dlm2, x=test$N05C, h=14)
mape.dlm2<- MAPE(fore.dlm2$forecasts, test$N05B)
#akurasi data training
GoF(model.dlm2)## n MAE MPE MAPE sMAPE MASE MSE MRAE
## model.dlm2 50 47.66669 -2.765934 2.921334 0.2131324 1.10173 3223.728 3.232601
## GMRAE
## model.dlm2 1.395469Didapatkan nilai MAPE sebesar \(2.92\%\). Model tersebut merupakan
model yang sangat baik dengan nilai MAPE
yang kurang dari \(10\%\).
3. Model Autoregressive
Peubah dependen dipengaruhi oleh peubah independen pada waktu sekarang, serta dipengaruhi juga oleh peubah dependen itu sendiri pada satu waktu yang lalu maka model tersebut disebut autoregressive (Gujarati 2004).
Pemodelan
Pemodelan Autoregressive dilakukan menggunakan fungsi
dLagM::ardlDlm() . Fungsi tersebut akan menerapkan
autoregressive berordo \((p,q)\) dengan satu prediktor. Fungsi umum
dari ardlDlm() adalah sebagai berikut.
ardlDlm(formula = NULL , data = NULL , x = NULL , y = NULL , p = 1 , q = 1 ,
remove = NULL )Dengan \(p\) adalah integer yang mewakili panjang lag yang terbatas dan \(q\) adalah integer yang merepresentasikan ordo dari proses autoregressive.
model.ardl <- ardlDlm(x = train$N05C, y = train$N05B, p = 1 , q = 1)
summary(model.ardl)##
## Time series regression with "ts" data:
## Start = 2, End = 56
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -172.28 -29.32 11.79 32.68 96.38
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.42282 22.22861 1.414 0.1635
## X.t 4.99458 0.95928 5.207 3.47e-06 ***
## X.1 -1.88602 0.96878 -1.947 0.0571 .
## Y.1 0.68407 0.09714 7.042 4.67e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 50.21 on 51 degrees of freedom
## Multiple R-squared: 0.713, Adjusted R-squared: 0.6961
## F-statistic: 42.22 on 3 and 51 DF, p-value: 7.407e-14AIC(model.ardl)## [1] 592.7032BIC(model.ardl)## [1] 602.7399Dari hasil tersebut, didapat bahwa peubah \(x_t\) , dan \(y_{t−1}\) memiliki nilai \(P.Value < 0.05\) Hal ini menunjukkan bahwa kedua peubah tersebut berpengaruh signifikan terhadap \(y_t\) pada taraf nyata \(5\%\). Model keseluruhannya adalah sebagai berikut:
\[ \hat{Y}=31.42282 + 4.99458 X_t - 1.88602 X_{t-1} + 0.68407 Y_{t-1} \]
Peramalan dan Akurasi
fore.ardl <- forecast(model = model.ardl, x=test$N05C, h=14)
fore.ardl## $forecasts
## [1] 226.7060 288.7369 291.6681 324.3042 317.2211 264.9794 284.9509 300.3944
## [9] 325.2790 329.7630 309.7435 335.4464 291.7670 235.7963
##
## $call
## forecast.ardlDlm(model = model.ardl, x = test$N05C, h = 14)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"Data di atas merupakan hasil peramalan untuk 5 periode ke depan menggunakan Model Autoregressive dengan \(p=1\) dan \(q=1\).
mape.ardl <- MAPE(fore.ardl$forecasts, test$N05B)
mape.ardl## [1] 0.3127644#akurasi data training
GoF(model.ardl)## n MAE MPE MAPE sMAPE MASE MSE MRAE
## model.ardl 55 36.87358 -3.136187 3.268783 0.1712242 0.830694 2337.242 1.717624
## GMRAE
## model.ardl 0.9326488Berdasarkan akurasi di atas, terlihat bahwa nilai MAPE
keduanya tidak jauh berbeda. Artinya, model regresi
dengan distribusi lag ini
tidak overfitted atau underfitted
Lag Optimum
#penentuan lag optimum
model.ardl.opt <- ardlBoundOrders(data = data.frame(data), ic = "AIC",
formula = N05B ~ N05C )
min_p=c()
for(i in 1:6){
min_p[i]=min(model.ardl.opt$Stat.table[[i]])
}
q_opt=which(min_p==min(min_p, na.rm = TRUE))
p_opt=which(model.ardl.opt$Stat.table[[q_opt]] ==
min(model.ardl.opt$Stat.table[[q_opt]], na.rm = TRUE))
data.frame("q_optimum" = q_opt, "p_optimum" = p_opt,
"AIC"=model.ardl.opt$min.Stat)## q_optimum p_optimum AIC
## 1 3 15 568.2968Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=15\) dan \(q=3\), yaitu sebesar \(568.2968\). Artinya, model autoregressive optimum didapat ketika \(p=15\) dan \(q=3\).
Selanjutnya dapat dilakukan pemodelan dengan nilai \(p\) dan \(q\) optimum seperti inisialisasi di langkah sebelumnya.
4. Pemodelan DLM & ARDL dengan Library dynlm
Pemodelan regresi dengan peubah lag tidak hanya dapat
dilakukan dengan fungsi pada packages dLagM ,
tetapi terdapat packages dynlm yang dapat
digunakan. Fungsi dynlm secara umum adalah sebagai
berikut.
dynlm(formula, data, subset, weights, na.action, method = "qr",
model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE,
contrasts = NULL, offset, start = NULL, end = NULL, ...)Untuk menentukan formula model yang akan digunakan,
tersedia fungsi tambahan yang memungkinkan spesifikasi dinamika (melalui
d() dan L()) atau pola linier/siklus dengan
mudah (melalui trend(), season(), dan
harmon()). Semua fungsi formula baru mengharuskan
argumennya berupa objek deret waktu (yaitu, "ts" atau
"zoo").
install_load('dynlm')
#sama dengan model dlm q=1
cons_lm1 <- dynlm(N05B ~ N05C+L(N05C),data = train.ts)
#sama dengan model ardl p=1 q=0
cons_lm2 <- dynlm(N05B ~ N05C+L(N05B),data = train.ts)
#sama dengan ardl p=1 q=1
cons_lm3 <- dynlm(N05B ~ N05C+L(N05C)+L(N05B),data = train.ts)
#sama dengan dlm p=2
cons_lm4 <- dynlm(N05B ~ N05C+L(N05C)+L(N05C,2),data = train.ts)Ringkasan Model
summary(cons_lm1)##
## Time series regression with "ts" data:
## Start = 2, End = 56
##
## Call:
## dynlm(formula = N05B ~ N05C + L(N05C), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -153.213 -47.059 0.404 31.168 184.815
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 119.862 25.509 4.699 1.96e-05 ***
## N05C 6.688 1.292 5.178 3.68e-06 ***
## L(N05C) 1.908 1.120 1.704 0.0943 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 69.83 on 52 degrees of freedom
## Multiple R-squared: 0.4338, Adjusted R-squared: 0.4121
## F-statistic: 19.92 on 2 and 52 DF, p-value: 3.77e-07summary(cons_lm2)##
## Time series regression with "ts" data:
## Start = 2, End = 56
##
## Call:
## dynlm(formula = N05B ~ N05C + L(N05B), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -178.614 -28.110 8.402 30.603 94.482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.84116 22.81509 1.352 0.182
## N05C 4.76028 0.97690 4.873 1.07e-05 ***
## L(N05B) 0.57888 0.08287 6.986 5.22e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 51.53 on 52 degrees of freedom
## Multiple R-squared: 0.6916, Adjusted R-squared: 0.6798
## F-statistic: 58.31 on 2 and 52 DF, p-value: 5.202e-14summary(cons_lm3)##
## Time series regression with "ts" data:
## Start = 2, End = 56
##
## Call:
## dynlm(formula = N05B ~ N05C + L(N05C) + L(N05B), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -172.28 -29.32 11.79 32.68 96.38
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.42282 22.22861 1.414 0.1635
## N05C 4.99458 0.95928 5.207 3.47e-06 ***
## L(N05C) -1.88602 0.96878 -1.947 0.0571 .
## L(N05B) 0.68407 0.09714 7.042 4.67e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 50.21 on 51 degrees of freedom
## Multiple R-squared: 0.713, Adjusted R-squared: 0.6961
## F-statistic: 42.22 on 3 and 51 DF, p-value: 7.407e-14summary(cons_lm4)##
## Time series regression with "ts" data:
## Start = 3, End = 56
##
## Call:
## dynlm(formula = N05B ~ N05C + L(N05C) + L(N05C, 2), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.311 -45.216 0.153 35.738 172.852
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 106.58449 29.86292 3.569 0.000802 ***
## N05C 6.86841 1.30848 5.249 3.13e-06 ***
## L(N05C) 2.65332 1.32657 2.000 0.050931 .
## L(N05C, 2) -0.03826 1.12790 -0.034 0.973071
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 70.1 on 50 degrees of freedom
## Multiple R-squared: 0.4427, Adjusted R-squared: 0.4093
## F-statistic: 13.24 on 3 and 50 DF, p-value: 1.749e-06SSE
cat(" Model DLM q=1 : ", deviance(cons_lm1), "\n",
"Model ARDL p=1 q=0 : ", deviance(cons_lm1), "\n",
"Model ARDL p=1 q=1 : ", deviance(cons_lm1), "\n",
"Model DLM p=2 : ", deviance(cons_lm1), "\n")## Model DLM q=1 : 253542.4
## Model ARDL p=1 q=0 : 253542.4
## Model ARDL p=1 q=1 : 253542.4
## Model DLM p=2 : 253542.4Uji Diagnostik
#uji model
install_load('lmtest')
encomptest(cons_lm1, cons_lm2)## Encompassing test
##
## Model 1: N05B ~ N05C + L(N05C)
## Model 2: N05B ~ N05C + L(N05B)
## Model E: N05B ~ N05C + L(N05C) + L(N05B)
## Res.Df Df F Pr(>F)
## M1 vs. ME 51 -1 49.59 4.671e-09 ***
## M2 vs. ME 51 -1 3.79 0.05708 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Autokorelasi
#durbin watson
dwtest(cons_lm1)##
## Durbin-Watson test
##
## data: cons_lm1
## DW = 0.62301, p-value = 9.3e-10
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm2)##
## Durbin-Watson test
##
## data: cons_lm2
## DW = 1.7446, p-value = 0.142
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm3)##
## Durbin-Watson test
##
## data: cons_lm3
## DW = 2.1159, p-value = 0.6111
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm4)##
## Durbin-Watson test
##
## data: cons_lm4
## DW = 0.65465, p-value = 5.485e-09
## alternative hypothesis: true autocorrelation is greater than 0Heterogenitas
bptest(cons_lm1)##
## studentized Breusch-Pagan test
##
## data: cons_lm1
## BP = 4.2819, df = 2, p-value = 0.1175bptest(cons_lm2)##
## studentized Breusch-Pagan test
##
## data: cons_lm2
## BP = 6.1905, df = 2, p-value = 0.04526bptest(cons_lm3)##
## studentized Breusch-Pagan test
##
## data: cons_lm3
## BP = 5.5969, df = 3, p-value = 0.133bptest(cons_lm4)##
## studentized Breusch-Pagan test
##
## data: cons_lm4
## BP = 7.2047, df = 3, p-value = 0.06565Kenormalan
shapiro.test(residuals(cons_lm1))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm1)
## W = 0.96605, p-value = 0.1219shapiro.test(residuals(cons_lm2))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm2)
## W = 0.95002, p-value = 0.02304shapiro.test(residuals(cons_lm3))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm3)
## W = 0.93754, p-value = 0.006708shapiro.test(residuals(cons_lm4))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm4)
## W = 0.97604, p-value = 0.3498B. Perbandingan Model
akurasi <- matrix(c(mape.koyck, mape.dlm, mape.dlm2, mape.ardl))
row.names(akurasi)<- c("Koyck","DLM 1","DLM 2","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi## MAPE
## Koyck 17.5008687
## DLM 1 15.0606495
## DLM 2 0.4461025
## Autoregressive 0.3127644Berdasarkan nilai MAPE, model paling
optimum didapat pada Model Autoregressive
karena memiliki nilai MAPE yang terkecil.
Plot
par(mfrow=c(1,1))
plot(test$N05C, test$N05B, type="b", col="black",
ylim=c(-600,380), xlim=c(5,30))
points(test$N05C, fore.koyck$forecasts, col="red")
lines(test$N05C, fore.koyck$forecasts, col="red")
points(test$N05C, fore.dlm$forecasts, col="blue")
lines(test$N05C, fore.dlm$forecasts, col="blue")
points(test$N05C, fore.dlm2$forecasts, col="orange")
lines(test$N05C, fore.dlm2$forecasts, col="orange")
points(test$N05C, fore.ardl$forecasts, col="green")
lines(test$N05C, fore.ardl$forecasts, col="green")
legend("topleft",c("aktual", "koyck","DLM 1","DLM 2", "autoregressive"),
lty=1, col=c("black","red","blue","orange","green"), cex=0.8)
Agak chaos yak hiks :“)
Berdasarkan plot tersebut, terlihat bahwa plot yang paling mendekati data aktualnya adalah Model Autoregressive, sehingga dapat disimpulkan model terbaik dalam hal ini adalah Model Reresi Autoregressive.
C. Kesimpulan
Penjualan obat tidur N05C periode sebelumnya
memiliki pengaruh signifikan terhadap penjualan obat
penenang N05B pada periode waktu saat ini obat. Model
Regresi terbaik degan Peubah Lag nya adalah Model Regresi
Autoregressive. Dengan urutan kebaikan model dari yang terbaik
adalah Autoregressive (MAPE =
\(0.3\) ; “Sangat Baik”),
DLM 2 (MAPE = \(0.45\) ; “Sangat Baik”),
DLM 1 (MAPE = \(15.06\) ;”Baik”), dan
Koyck (MAPE = \(17.5\) ; “Baik”).