Data Hortikultura
Rachmat Darmawan
12/30/2021
VAR Model - untuk data Hasil Produksi Cabai Rawit di Sulteng
# load package
library(tsModel)## Warning: package 'tsModel' was built under R version 4.0.5## Time Series Modeling for Air Pollution and Health (0.6)library(vars)## Warning: package 'vars' was built under R version 4.0.5## Loading required package: MASS## Loading required package: strucchange## Warning: package 'strucchange' was built under R version 4.0.5## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.0.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Warning: package 'sandwich' was built under R version 4.0.4## Loading required package: urca## Warning: package 'urca' was built under R version 4.0.3## Loading required package: lmtest## Warning: package 'lmtest' was built under R version 4.0.4library(tseries)## Warning: package 'tseries' was built under R version 4.0.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(dynlm)## Warning: package 'dynlm' was built under R version 4.0.5library(tidyverse)## -- Attaching packages ---------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.2
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## Warning: package 'ggplot2' was built under R version 4.0.5## -- Conflicts ------------------------------------------------------- tidyverse_conflicts() --
## x stringr::boundary() masks strucchange::boundary()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks MASS::select()library(urca)# load data and create time series object ----
dat <- read.csv("D:/produksi_rawit.csv")
head(dat)## Periode Luas.Lahan.Baru Hasil.Produksi
## 1 2014/01 612 9216
## 2 2014/02 369 13547
## 3 2014/03 231 10191
## 4 2014/04 140 12422
## 5 2014/05 334 10156
## 6 2014/06 162 11381tail(dat)## Periode Luas.Lahan.Baru Hasil.Produksi
## 79 2020/07 337 20985
## 80 2020/08 298 22626
## 81 2020/09 362 23329
## 82 2020/10 266 20577
## 83 2020/11 246 21274
## 84 2020/12 200 20654lahanbaru <- ts(dat$Luas.Lahan.Baru, start = c(2014, 2), freq = 12)
produksi <- ts(dat$Hasil.Produksi, start = c(2014, 2), freq = 12)Eksplorasi Data
plot(cbind(produksi, lahanbaru))
Plot data produksi cabai rawit di Sulawesi Tengah menunjukkan hasil produksi cenderung berfluktuatif dari waktu ke waktu. Produksi tertinggi terjadi pada Mei 2019, sedangkan produksi terendah terjadi pada Januari 2016. Plot data luas penanaman baru menunjukkan bahwa luas penanaman baru juga berfluktuatif dari waktu ke waktu dimana luas penanaman baru tertinggi terjadi pada Januari 2014 dan terendah terjadi pada November 2015.
ggplot (data = dat) + geom_point(mapping = aes(x = Luas.Lahan.Baru, y = Hasil.Produksi))
Berdasarkan scatter plot di atas, dapat dilihat bahwa luas lahan baru memiliki hubungan linear poitif dengan hasil produksi. Sebagai dugaan awal, semakin banyak lahan baru yang dibuka untuk penanaman, maka akan menaikkan hasil produksi cabai rawit.
# look at the persistence in the data ----
lahanbaru.acf <- acf(lahanbaru, main = "output")
lahanbaru.pacf <- pacf(lahanbaru, main = "output")
produksi.acf <- acf(produksi, main = "produksi")
produksi.pacf <- pacf(produksi, main = "produksi")
Pengujian Kestasioneran
adf.test(lahanbaru)## Warning in adf.test(lahanbaru): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: lahanbaru
## Dickey-Fuller = -4.8021, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryadf.test(produksi)##
## Augmented Dickey-Fuller Test
##
## data: produksi
## Dickey-Fuller = -2.7759, Lag order = 4, p-value = 0.2576
## alternative hypothesis: stationaryproduksi.dif1 <- diff(produksi, difference=1)
adf.test(produksi.dif1)## Warning in adf.test(produksi.dif1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: produksi.dif1
## Dickey-Fuller = -5.4479, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryHasil pengujian Augmented Dickey Fuller peubah lahan baru menghasilkan p-value sebesar 0.01. Nilai tersebut lebih kecil dari taraf nyata 0.05 sehingga data sudah stasioner dalam rataan. Sementara peubah produksi menghasilkan p-value sebesar 0.257 sehingga perlu dilakukan proses differencing untuk mengatasi ketidakstasioneran tersebut. Pengujian Augmented Dickey Fuller pada deret output produksi cabai rawit setelah differencing satu kali menghasilkan p-value sebesar 0.01. Nilai tersebut lebih kecil dari taraf nyata 0.05 sehingga peubah sudah stasioner dalam rataan. Model yang dapat digunakan karena adanya pembedaan ini adalah model VAR dengan pembedaan untuk data tak terkointegrasi atau model VECM untuk data terkointegrasi.
Pengujian Kointegrasi
cointcy <- dynlm(lahanbaru~produksi.dif1)
ehat <- resid(cointcy)
adf.test(ehat)##
## Augmented Dickey-Fuller Test
##
## data: ehat
## Dickey-Fuller = -3.5405, Lag order = 4, p-value = 0.04369
## alternative hypothesis: stationaryHal tersebut menunjukkan bahwa tidak terdapat hubungan linear jangka panjang antara ketiga peubah yang bersifat stasioner. Hasil ini mengindikasikan bahwa model yang akan digunakan adalah model VAR ordo 2 dengan pembedaan 1 kali.
jotest=ca.jo(data.frame(lahanbaru,produksi), type="trace", K=2, ecdet="none", spec="longrun")
summary(jotest)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.30727884 0.06980445
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 5.93 6.50 8.18 11.65
## r = 0 | 36.04 15.66 17.95 23.52
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## lahanbaru.l2 produksi.l2
## lahanbaru.l2 1.000000000 1.00000000
## produksi.l2 -0.009318023 0.02771273
##
## Weights W:
## (This is the loading matrix)
##
## lahanbaru.l2 produksi.l2
## lahanbaru.d -0.4469099 -0.04465401
## produksi.d 25.1089674 -4.84682552Pemilihan Ordo VAR
Ordo VAR ditentukan melalui perhitungan nilai AIC untuk setiap lag. Dari keseluruhan perhitungan nilai AIC, terlihat bahwa nilai terkecil AIC tercapai pada lag 2 (p=2). Hal ini mengindikasikan bahwa model VAR yang digunakan adalah model VAR ordo ke-2.
# group data and choose lag length
dat.bv <- cbind(lahanbaru, produksi)
colnames(dat.bv) <- c("LahanBaru", "Produksi")info.bv <- VARselect(dat.bv, lag.max = 12, type = "const")
info.bv$selection # All suggest suggests using 2 lags## AIC(n) HQ(n) SC(n) FPE(n)
## 2 1 1 2# group data and choose lag length
bv.est <- VAR(dat.bv, p = 2, type = "const", season = NULL, exog = NULL)
summary(bv.est)##
## VAR Estimation Results:
## =========================
## Endogenous variables: LahanBaru, Produksi
## Deterministic variables: const
## Sample size: 82
## Log Likelihood: -1248.724
## Roots of the characteristic polynomial:
## 0.8559 0.2732 0.2698 0.08473
## Call:
## VAR(y = dat.bv, p = 2, type = "const", exogen = NULL)
##
##
## Estimation results for equation LahanBaru:
## ==========================================
## LahanBaru = LahanBaru.l1 + Produksi.l1 + LahanBaru.l2 + Produksi.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## LahanBaru.l1 4.979e-01 1.143e-01 4.355 4.05e-05 ***
## Produksi.l1 2.103e-03 1.397e-03 1.505 0.13636
## LahanBaru.l2 1.058e-02 9.685e-02 0.109 0.91327
## Produksi.l2 8.241e-04 1.372e-03 0.601 0.54985
## const 6.670e+01 2.357e+01 2.830 0.00593 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 56.92 on 77 degrees of freedom
## Multiple R-Squared: 0.4296, Adjusted R-squared: 0.4
## F-statistic: 14.5 on 4 and 77 DF, p-value: 7.187e-09
##
##
## Estimation results for equation Produksi:
## =========================================
## Produksi = LahanBaru.l1 + Produksi.l1 + LahanBaru.l2 + Produksi.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## LahanBaru.l1 11.3935 9.0437 1.260 0.211535
## Produksi.l1 0.4461 0.1105 4.037 0.000127 ***
## LahanBaru.l2 8.8686 7.6626 1.157 0.250687
## Produksi.l2 0.1856 0.1086 1.710 0.091343 .
## const 1231.9500 1864.4353 0.661 0.510736
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 4503 on 77 degrees of freedom
## Multiple R-Squared: 0.5218, Adjusted R-squared: 0.4969
## F-statistic: 21 on 4 and 77 DF, p-value: 9.762e-12
##
##
##
## Covariance matrix of residuals:
## LahanBaru Produksi
## LahanBaru 3240 9951
## Produksi 9951 20277587
##
## Correlation matrix of residuals:
## LahanBaru Produksi
## LahanBaru 1.00000 0.03882
## Produksi 0.03882 1.00000# testing serial correlation with Portmanteau-Test
bv.serial <- serial.test(bv.est, lags.pt = 12, type = "PT.asymptotic")
bv.serial##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object bv.est
## Chi-squared = 24.108, df = 40, p-value = 0.9778# p-value greater than 5% indicates absence of serial correlationwindows()
plot(bv.serial, names = "LahanBaru")
windows()
plot(bv.serial, names = "Produksi")
Uji Heteroskedastisitas
bv.arch <- arch.test(bv.est, lags.multi = 12, multivariate.only = TRUE)
bv.arch##
## ARCH (multivariate)
##
## data: Residuals of VAR object bv.est
## Chi-squared = 86.583, df = 108, p-value = 0.9358# p-value greater than 5% indicates absence of heteroscedasticityUji Normalitas
bv.norm <- normality.test(bv.est, multivariate.only = TRUE)
bv.norm## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object bv.est
## Chi-squared = 149.19, df = 4, p-value < 2.2e-16
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object bv.est
## Chi-squared = 35.76, df = 2, p-value = 1.717e-08
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object bv.est
## Chi-squared = 113.43, df = 2, p-value < 2.2e-16# p-value indicates pretty normal distributionUji Stabilitas Struktur
# testing for structural stability
bv.cusum <- stability(bv.est, type = "OLS-CUSUM")
windows()
plot(bv.cusum)
# does not break the confidence intervals for structural breakTerlihat bahwa nilai pada grafik berada dalam selang kepercayaan sehingga model yang dimiliki cukup stabil untuk digunakan.
## Granger & instantaneous causality tests ====
bv.cause.lb <- causality(bv.est, cause = "LahanBaru")
bv.cause.lb## $Granger
##
## Granger causality H0: LahanBaru do not Granger-cause Produksi
##
## data: VAR object bv.est
## F-Test = 3.3065, df1 = 2, df2 = 154, p-value = 0.03926
##
##
## $Instant
##
## H0: No instantaneous causality between: LahanBaru and Produksi
##
## data: VAR object bv.est
## Chi-squared = 0.12341, df = 1, p-value = 0.7254bv.cause.pr <- causality(bv.est, cause = "Produksi")
bv.cause.pr## $Granger
##
## Granger causality H0: Produksi do not Granger-cause LahanBaru
##
## data: VAR object bv.est
## F-Test = 2.9072, df1 = 2, df2 = 154, p-value = 0.05763
##
##
## $Instant
##
## H0: No instantaneous causality between: Produksi and LahanBaru
##
## data: VAR object bv.est
## Chi-squared = 0.12341, df = 1, p-value = 0.7254# Null hypothesis of no Granger causality is dismissed in both directions#graphics.off() # close all graphsFungsi Respon Impuls
## Generate irfs ====
irf.lb <- irf(bv.est, impulse = "Produksi", response = "LahanBaru",n.ahead = 40, boot = TRUE)
plot(irf.lb) # reponse of gdp to unemploy shocks
# positive shock to unemployment has a negative affect on ouput (lower spending power)irf.pr <- irf(bv.est, impulse = "LahanBaru", response = "Produksi",n.ahead = 40, boot = TRUE)
plot(irf.pr) # reponse of unemployment to gdp shock
# positive shock to output decreases unemployment by large and persistent amountirf.pr_un <- irf(bv.est, impulse = "Produksi", response = "Produksi",n.ahead = 40, boot = TRUE)
plot(irf.pr, ylab = "Produksi", main = "Shock dari LahanBaru", xlab = "")
plot(irf.pr_un, ylab = "Produksi", main = "Shock dari Produksi")
Ketika terjadi perubahan satu satuan pada hasil produksi pada waktu ke-t akan menyebabkan kenaikan pada produksi untuk peramalan dua bulan ke depan. Kemudian untuk peramalan tiga bulan ke depan, perubahan pada hasil produksi produksi mulai menyebabkan kenaikan pada luas lahan baru. Kondisi ini terus berubah (mengalami penurunan) hingga periode ke-10. Pada periode peramalan bulan ke-11 sampai ke-12 perubahan hasil produksi terhadap luas lahan baru dalam VAR cenderung konstan dan mendekati nol. Hal ini menandakan bahwa perubahan yang terjadi pada peubah hasil produksi sudah tidak lagi memberikan pengaruh pada peubah dependen lainnya.
Dekomposisi Ragam
# generate variance decompositions ----
bv.vardec <- fevd(bv.est, n.ahead = 10)
windows()
plot(bv.vardec)
Terlihat bahwa luas lahan baru dan hasil produksi sebenarnya dipengaruhi peubah itu sendiri, kemudian hasil produksi cabai rawit di Sulawei Tengah dipengaruhi oleh luas penanaman lahan baru sampai tingkat tertentu. Hasil dekomposisi ragam ini cenderung stabil hingga akhir periode. Hasil ini menunjukkan hubungan antar peubah akibat adanya perubahan.
#graphics.off() # close graphics## Forecasting with the VAR model ====
predictions <- predict(bv.est, n.ahead = 8, ci = 0.95)
plot(predictions, names = "LahanBaru")
plot(predictions, names = "Produksi")
# Fanchart for forecasts ----
fanchart(predictions, names = "LahanBaru")
fanchart(predictions, names = "Produksi")
Visualisasi plot di atas menunjukkan bahwa awal tahun 2021 ada perluasan lahan penanaman baru, sementara hasil produksi cabai rawit diperkirakan akan mengalami penurunan disepanjang tahun 2021.
Kesimpulan
Model yang digunakan adalah model VAR dengan pembedaan 1 kali sampai ordo ke-2 atau VARD(2). Hasil analisis VARD(2) secara keseluruhan, termasuk fungsi respon impuls dan dekomposisi ragam menunjukkan bahwa peubah pembukaan luas lahan baru cukup memberikan pengaruh terhadap hasil produksi cabai rawit, demikian pula sebaliknya. Kemudian berdasarkan hasil peramalan diketahui bahwa setahun yang akan datang akan terjadi penurunan produksi cabai rawit di Provinsi Sulawesi Tengah yang berbanding terbalik dengan perluasan daerah penanaman baru.