1.1 Library yang Digunakan

Terdapat beberapa library yang akan digunakan untuk menganalisis data time series, di antaranya :

library(rmarkdown)
library(dLagM)

## Loading required package: nardl

## 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

library(dynlm) 
library(MLmetrics)

## 
## Attaching package: 'MLmetrics'

## The following object is masked from 'package:dLagM':
## 
##     MAPE

## The following object is masked from 'package:base':
## 
##     Recall

library(lmtest)
library(tseries)
library(lmtest)
library(car)

## Loading required package: carData

library(readxl)
library(readr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(corrplot)

## corrplot 0.92 loaded

1.2 Data yang Digunakan

Data yang digunakan merupakan data penjualan properti tahun 2007 — 2019 untuk satu wilayah tertentu per 3 bulan. Data perjualan properti ini berasal dari Kaggle. Data tersebut berisi harga jual properti (MA), yang selanjutnya disebut sebagai peubah Y, dan jumlah kamar tidur (bedrooms), yang selanjutnya disebut sebagai peubah X.

dataforecast <-read_excel("C:/ICA/Semester 5/MPDW/Data P4.xlsx")
dataforecast <- dataforecast[c(1,2,4)]
head(dataforecast)

## # A tibble: 6 x 3
##   saledate                MA bedrooms
##   <dttm>               <dbl>    <dbl>
## 1 2007-09-30 00:00:00 441854        2
## 2 2007-12-31 00:00:00 441854        2
## 3 2008-03-31 00:00:00 441854        2
## 4 2008-06-30 00:00:00 441854        2
## 5 2008-09-30 00:00:00 451583        2
## 6 2008-12-31 00:00:00 440256        2

1.3 Eksplorasi Data

corrplot::corrplot(corr=cor(dataforecast[sapply(dataforecast,is.numeric)]), method = "number", type = "upper")

Dapat dilihat bahwa terdapat hubungan positif antara peubah MA dan Bedrooms. Nilai korelasi antara kedua variabel menunjukkan bahwa korelasi tersebut sangat kuat, yaitu sebesar 0.92.

1.4 Split Data

#SplitData
train <- dataforecast[1:279,]
test <- dataforecast[280:398,]

#DataTimeSeries
train.ts<-ts(train)
test.ts<-ts(test)

#ConvertingData
train$MA <- as.numeric(train$MA)
train$bedrooms <- as.numeric(train$bedrooms)

Data dibagi menjadi 2 bagian, yaitu data training dan data testing yang masing-masing sebesar 80% (279 baris) dan 20% sisanya (118 baris). Setelah itu, data dijadikan time series dengan fungsi ts serta di-convert menjadi numeric.

plot(train$bedrooms,train$MA,pch = 20, col = "red", main = "Scatter Plot Bedrooms VS Price (Test)", ylab = "Price", xlab = "Bedrooms")

plot(test$bedrooms,test$MA,pch = 20, col = "red", main = "Scatter Plot Bedrooms VS Price (Test)", ylab = "Price", xlab = "Bedrooms")

2.1 Regression with Distributed Lag Optimum

#penentuan lag optimum 
finiteDLMauto(formula = MA ~ bedrooms,
              data = data.frame(train), q.min = 1, q.max = 4,
              model.type = "dlm", error.type = "AIC", trace = TRUE) ##q max lag maksimum

##   q - k    MASE      AIC      BIC    GMRAE   MBRAE R.Adj.Sq Ljung-Box
## 4     4 7.13260 6939.388 6964.706 81.38299 1.02662  0.85244         0
## 3     3 7.08890 6964.125 6985.847 78.88007 1.04730  0.85182         0
## 2     2 7.07639 6988.801 7006.921 86.20175 1.07836  0.85123         0
## 1     1 7.06236 7013.419 7027.930 93.55470 0.89357  0.85067         0

Berdasarkan hasil di atas, dapat disimpulkan bahwa semakin besar besar laq maka nilai AIC dan BIC semakin kecil. Akan tetapi, hal ini tidak berarti nilai lag yang besar memberikan model regresi terbaik. Pada penelitian ini, peneliti menggunakan nilai lag maksimum 4 sebagai model regresi lag optimum.

model.dlm = dLagM::dlm(x = train$bedrooms,y = train$MA, q = 4) 
summary(model.dlm)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -148777  -55623    5498   36284  196594 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.803e+05  1.040e+04  17.338  < 2e-16 ***
## x.t          1.026e+05  1.617e+04   6.347  9.3e-10 ***
## x.1         -1.633e-10  2.275e+04   0.000    1.000    
## x.2         -2.895e-11  2.275e+04   0.000    1.000    
## x.3          4.905e-11  2.275e+04   0.000    1.000    
## x.4          2.555e+04  1.617e+04   1.580    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 71950 on 269 degrees of freedom
## Multiple R-squared:  0.8551, Adjusted R-squared:  0.8524 
## F-statistic: 317.6 on 5 and 269 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 6939.388 6964.706

Hasil di atas menunjukkan bahwa peubah hanya Xt (jumlah kamar pada waktu sekarang) yang berpengaruh signifikan terhadap harga properti. Diperoleh pula nilai AIC dan BIC yang masing-masing sebesar 6939.388 dan 6964.706. Dengan begitu, diperoleh model regresi dengan distribusi lag sebagai berikut : Yt = 180300 + 102600Xt − 0.0000002Xt−1 − 0.0000003Xt−2 + 0.0000005Xt−3 + 25550Xt-4

#RamalanLagOptimum
(fore.dlm <- dLagM::forecast(model = model.dlm, x=test$bedrooms, h=24)) #meramalkan 24 periode ke depan

## $forecasts
##  [1] 436633.1 436633.1 436633.1 436633.1 436633.1 436633.1 436633.1 436633.1
##  [9] 436633.1 436633.1 436633.1 436633.1 436633.1 436633.1 436633.1 436633.1
## [17] 436633.1 436633.1 436633.1 539254.5 539254.5 539254.5 539254.5 564808.3
## 
## $call
## forecast.dlm(model = model.dlm, x = test$bedrooms, h = 24)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

Data di atas merupakan data peramalan harga properti untuk 24 periode ke depan menggunakan model regresi dengan distribusi lag sebesar 4.

#AkurasiDataTraining
mape_train <- GoF(model.dlm)["MAPE"]
mape_train

##                MAPE
## model.dlm 0.1007105

#PlotPeramalan
ts.plot(fore.dlm$forecasts, xlab="periode  waktu", ylab="MA", main="Plot Hasil Peramalan MA", col="blue", lty=3)
points(fore.dlm$forecasts)
lines (fore.dlm$forecasts, col="blue",lwd=2)
points(test$MA,col="red")
lines (test$MA, col="red",lwd=2)
legend("topleft",c("aktual", "DLM"), lty=1, col=c("red","blue"), cex=0.8)

Berdasarkan plot peramalan di atas, terlihat bahwa hasil peramalan DLM (garis berwarna biru) cukup mendekati data aktual (garis berwarna merah).

2.2 Model KOYCK

Metode Koyck didasari asumsi bahwa semakin jauh jarak lag pada variabel bebas dari periode sekarang maka semakin kecil pengaruh variabel lag terhadap variabel 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 = dLagM::koyckDlm(x = train$bedrooms, y = train$MA)
summary(model.koyck)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -590819   -2264    1096    6032   31372 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.951e+04  7.735e+03   2.522   0.0122 *  
## Y.1         9.152e-01  3.047e-02  30.041   <2e-16 ***
## X.t         9.407e+03  4.249e+03   2.214   0.0277 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37730 on 275 degrees of freedom
## Multiple R-Squared: 0.9594,  Adjusted R-squared: 0.9591 
## Wald test:  3237 on 2 and 275 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha     beta     phi
## Geometric coefficients:  230026.8 9407.051 0.91519

AIC(model.koyck)

## [1] 6653.179

BIC(model.koyck)

## [1] 6667.689

#RamalanKoyck
fore.koyck <- dLagM::forecast(model = model.koyck, x=test$bedrooms, h=24)
fore.koyck$forecasts

##  [1] 431700.6 433410.8 434975.9 436408.2 437719.1 438918.8 440016.8 441021.6
##  [9] 441941.2 442782.9 443553.1 444258.1 444903.2 445493.6 446034.0 446528.5
## [17] 446981.1 447395.3 447774.4 457528.3 466455.1 474624.7 482101.5 488944.2

#AkurasiDataTraining
mape_train <- dLagM::GoF(model.koyck)["MAPE"]
mape_train

##                   MAPE
## model.koyck 0.01813097

#PlotPeramalan
ts.plot(fore.koyck$forecasts, xlab="periode  waktu", ylab="Harga Properti", main="Plot Hasil Peramalan Harga Properti", col="blue", lty=3)
points(fore.koyck$forecasts)
lines (fore.koyck$forecasts, col="blue",lwd=2)
points(test$MA,col="red")
lines (test$MA, col="red",lwd=2)
legend("topleft",c("aktual", "Koyck"), lty=1, col=c("red","blue"), cex=0.8)

2.3 Model Autoregressive

Autoregressive adalah bentuk regresi yang menghubungkan nilai-nilai sebelumnya pada selang waktu (lag time) yang bermacam-macam. Sehingga model autoregressive, disebut juga dynamic regression, menyatakan suatu ramalan sebagai fungsi nilai-nilai sebelumnya dari data time series periode tertentu (Pawestri et al. 2018).

model.ardl = ardlDlm(x = train$bedrooms, y = train$MA, p = 3 , q = 3)
summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 4, End = 279
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -202730   -1002    2140    6077   24055 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.096e+04  5.204e+03   2.106   0.0362 *  
## X.t          1.177e+05  5.412e+03  21.750   <2e-16 ***
## X.1         -1.227e+05  1.029e+04 -11.933   <2e-16 ***
## X.2          3.191e+03  1.273e+04   0.251   0.8023    
## X.3          1.148e+04  9.103e+03   1.261   0.2085    
## Y.1          1.032e+00  6.094e-02  16.931   <2e-16 ***
## Y.2         -2.354e-02  8.771e-02  -0.268   0.7886    
## Y.3         -7.874e-02  6.091e-02  -1.293   0.1972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23970 on 268 degrees of freedom
## Multiple R-squared:  0.984,  Adjusted R-squared:  0.9836 
## F-statistic:  2356 on 7 and 268 DF,  p-value: < 2.2e-16

Berdasarkan hasil di atas, diperoleh peubah yang signifikan pada taraf nyata 5% adalah Harga Properti periode ini dan 1 periode lalu serta jumlah kamar tidur 1 periode lalu.

AIC(model.ardl)

## [1] 6359.722

BIC(model.ardl)

## [1] 6392.306

Diperoleh nilai AIC dan BIC pada model autoregressive masing-masing bernilai 6359.722 dan 6392.306.

#RamalanAutoregressive
(fore.ardl <- dLagM::forecast(model = model.ardl, x=test$bedrooms, h=24))

## $forecasts
##  [1] 429799.4 429779.4 429687.1 429594.9 429503.5 429418.7 429340.5 429269.1
##  [9] 429203.9 429144.5 429090.4 429041.1 428996.1 428955.2 428917.9 428883.9
## [17] 428852.9 428824.7 428799.0 546483.1 545162.5 544223.2 545495.3 546934.0
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$bedrooms, h = 24)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

Data di atas merupakan data peramalan harga properti untuk 24 periode ke depan menggunakan model regresi dengan Model Autoregressive.

#AkurasiDataTraining
mape_train <- dLagM::GoF(model.koyck)["MAPE"]
mape_train

##                   MAPE
## model.koyck 0.01813097

#PlotPeramalan
ts.plot(fore.ardl$forecasts, xlab="periode  waktu", ylab="Harga Properti", main="Plot Hasil Peramalan Harga Properti", col="blue", lty=3, ylim=c(-2000,40000))
points(fore.ardl$forecasts)
lines (fore.ardl$forecasts, col="blue",lwd=2)
points(test$MA,col="red")
lines (test$MA, col="red",lwd=2)
legend("topleft",c("aktual", "autoregressive"), lty=1, col=c("red", "blue"), cex=0.8)

2.4 Uji Diagnostik Model

#sama dengan model dlm q=1
cons_lm1 <- dynlm(MA ~ bedrooms+L(bedrooms),data = train.ts)

#sama dengan model dlm q=4
cons_lm2 <- dynlm(MA ~ bedrooms+L(bedrooms)+L(bedrooms,2)+L(bedrooms,3)+L(bedrooms,4),data = train.ts)

#sama dengan ardl p=3 q=3
cons_lm3 <- dynlm(MA ~ bedrooms+L(bedrooms)+L(bedrooms,2)+L(bedrooms,3)+L(MA)+L(MA,2)+L(MA,3),data = train.ts)

#Ringkasan Model
summary(cons_lm1)

## 
## Time series regression with "ts" data:
## Start = 2, End = 279
## 
## Call:
## dynlm(formula = MA ~ bedrooms + L(bedrooms), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -148691  -53584    4906   35902  199739 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   184784      10236  18.052  < 2e-16 ***
## bedrooms      101857      16206   6.285 1.28e-09 ***
## L(bedrooms)    24789      16206   1.530    0.127    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 72130 on 275 degrees of freedom
## Multiple R-squared:  0.8517, Adjusted R-squared:  0.8507 
## F-statistic:   790 on 2 and 275 DF,  p-value: < 2.2e-16

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 5, End = 279
## 
## Call:
## dynlm(formula = MA ~ bedrooms + L(bedrooms) + L(bedrooms, 2) + 
##     L(bedrooms, 3) + L(bedrooms, 4), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -148777  -55623    5498   36284  196594 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.803e+05  1.040e+04  17.338  < 2e-16 ***
## bedrooms        1.026e+05  1.617e+04   6.347  9.3e-10 ***
## L(bedrooms)    -1.633e-10  2.275e+04   0.000    1.000    
## L(bedrooms, 2) -2.895e-11  2.275e+04   0.000    1.000    
## L(bedrooms, 3)  4.905e-11  2.275e+04   0.000    1.000    
## L(bedrooms, 4)  2.555e+04  1.617e+04   1.580    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 71950 on 269 degrees of freedom
## Multiple R-squared:  0.8551, Adjusted R-squared:  0.8524 
## F-statistic: 317.6 on 5 and 269 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 4, End = 279
## 
## Call:
## dynlm(formula = MA ~ bedrooms + L(bedrooms) + L(bedrooms, 2) + 
##     L(bedrooms, 3) + L(MA) + L(MA, 2) + L(MA, 3), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -202730   -1002    2140    6077   24055 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.096e+04  5.204e+03   2.106   0.0362 *  
## bedrooms        1.177e+05  5.412e+03  21.750   <2e-16 ***
## L(bedrooms)    -1.227e+05  1.029e+04 -11.933   <2e-16 ***
## L(bedrooms, 2)  3.191e+03  1.273e+04   0.251   0.8023    
## L(bedrooms, 3)  1.148e+04  9.103e+03   1.261   0.2085    
## L(MA)           1.032e+00  6.094e-02  16.931   <2e-16 ***
## L(MA, 2)       -2.354e-02  8.771e-02  -0.268   0.7886    
## L(MA, 3)       -7.874e-02  6.091e-02  -1.293   0.1972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23970 on 268 degrees of freedom
## Multiple R-squared:  0.984,  Adjusted R-squared:  0.9836 
## F-statistic:  2356 on 7 and 268 DF,  p-value: < 2.2e-16

#SSE
deviance(cons_lm1)

## [1] 1.43057e+12

deviance(cons_lm2)

## [1] 1.392503e+12

deviance(cons_lm3)

## [1] 153932187214

2.5 Uji Model

dwtest(cons_lm1)

## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 0.12602, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm2)

## 
##  Durbin-Watson test
## 
## data:  cons_lm2
## DW = 0.12883, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm3)

## 
##  Durbin-Watson test
## 
## data:  cons_lm3
## DW = 2.0172, p-value = 0.5073
## alternative hypothesis: true autocorrelation is greater than 0

bptest(cons_lm1)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 59.858, df = 2, p-value = 1.005e-13

bptest(cons_lm2)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 56.371, df = 5, p-value = 6.814e-11

bptest(cons_lm3)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm3
## BP = 9.1479, df = 7, p-value = 0.2422

shapiro.test(residuals(cons_lm1))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm1)
## W = 0.97121, p-value = 2.212e-05

shapiro.test(residuals(cons_lm2))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.97324, p-value = 5.035e-05

shapiro.test(residuals(cons_lm3))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm3)
## W = 0.33548, p-value < 2.2e-16

2.5 Perbandingan Keakuratan Ramalan

Koyck_Tabel <- matrix(c(AIC(model.koyck),BIC(model.koyck)))

## [1] 6653.179
## [1] 6667.689

row.names(Koyck_Tabel)<- c("AIC", "BIC")
colnames(Koyck_Tabel) <- c("Nilai")
Koyck_Tabel

##        Nilai
## AIC 6653.179
## BIC 6667.689

DLM_Tabel <- matrix(c(AIC(model.dlm),BIC(model.dlm)))

## [1] 6939.388
## [1] 6964.706

row.names(DLM_Tabel)<- c("AIC", "BIC")
colnames(DLM_Tabel) <- c("Nilai")
DLM_Tabel

##        Nilai
## AIC 6939.388
## BIC 6964.706

ARDL_Tabel <- matrix(c(AIC(model.ardl),BIC(model.ardl)))

## [1] 6359.722
## [1] 6392.306

row.names(ARDL_Tabel)<- c("AIC", "BIC")
colnames(ARDL_Tabel) <- c("Nilai")
ARDL_Tabel

##        Nilai
## AIC 6359.722
## BIC 6392.306

Berdasarkan ketiga tabel di atas, dapat dilihat bahwa nilai AIC dan BIC terkecil dimiliki oleh Model Autoregressive. Artinya, Model Autoregressive merupakan model paling ideal dibandingkan model lainnya.

Daftar Pustaka

Gujarati DN. 2004. Basic Econometrics. Fourth Edition. New York: The McGrawHill.

Pawestri V, Setiawan A, Linawati L. 2019. Pemodelan data penjualan mobil menggunakan model autoregressive moving average berdasarkan metode bayesian. Jurnal Sains dan Eduka Sains. 2(1): 26-35.