PEMODELAN LUAS PANEN PADI NASIONAL DENGAN METODE GARCH (Generalized Autoregressive Conditional Heteroscedastic) TEUKU ACHMAD IQBAL

PEMODELAN LUAS PANEN PADI NASIONAL DENGAN METODE GARCH (Generalized Autoregressive Conditional Heteroscedastic)

TEUKU ACHMAD IQBAL

SEKOLAH PASCASARJANA INSTITUT PERTANIAN BOGOR BOGOR 2014

PERNYATAAN MENGENAI TESIS DAN SUMBER INFORMASI SERTA PELIMPAHAN HAK CIPTA Dengan ini saya menyatakan bahwa tesis berjudul Pemodelan Luas Panen Padi Nasional dengan Metode GARCH (Generalized Autoregressive Conditional Heteroscedastic) adalah benar karya saya dengan arahan dari komisi pembimbing dan belum diajukan dalam bentuk apa pun kepada perguruan tinggi mana pun. Sumber informasi yang berasal atau dikutip dari karya yang diterbitkan maupun tidak diterbitkan dari penulis lain telah disebutkan dalam teks dan dicantumkan dalam Daftar Pustaka di bagian akhir disertasi ini. Dengan ini saya melimpahkan hak cipta dari karya tulis saya kepada Institut Pertanian Bogor. Bogor, Januari 2014 Teuku Achmad Iqbal NIM G152100061

RINGKASAN TEUKU ACHMAD IQBAL. Pemodelan Luas Panen Padi Nasional dengan Metode GARCH (Generalized Autoregressive Conditional Heteroscedastic). Dibimbing oleh KUSMAN SADIK dan I MADE SUMERTAJAYA. Data deret waktu luas panen padi nasional memiliki volatilitas yang tinggi dan ragam yang tidak homogen menurut waktu. Data deret waktu dengan ragam yang tidak homogen di setiap waktunya dinamakan data deret waktu dengan heteroskedastisitas bersyarat. Metode analisis deret waktu yang dapat digunakan untuk mengatasi heteroskedastisitas diantaranya adalah model GARCH. Akan tetapi pada data luas panen padi terdapat kemungkinan asimetris dalam volatilitasnya. Untuk mengatasi pengaruh asimetri, beberapa model GARCH sisaan asimetri dapat digunakan, antara lain: model GARCH sisaan eksponensial asimetris (EGARCH), model GARCH sisaan kuadratik asimetri model (QGARCH), model T-GARCH dan model GARCH sisaan non-linier asimetri (NAGARCH). Penelitian ini bertujuan untuk memodelkan luas panen padi nasional dengan cara memasukkan unsur keheterogenan ragam dan pengaruh asimetri pada data dengan menggunakan lima jenis model GARCH simetri, asimetri, dan non-linier, kemudian mendapatkan model terbaik dari lima jenis model GARCH tersebut. Model yang sesuai untuk luas panen padi nasional adalah model ARIMA(2,0,0)(1,1,0)12 ‒ GARCH(1,2) dan model ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2). Berdasarkan nilai mean absolute percentage error (MAPE) hingga dua puluh dua periode ke depan, model ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2) lebih baik dibandingkan dengan model ARIMA(2,0,0)(1,1,0)12 ‒ GARCH(1,2). Namun nilai MAPE untuk kedua model tersebut cukup tinggi karena terdapat beberapa nilai prediksi yang menyimpang cukup jauh dari nilai aktual. Akan tetapi, nilai MAPE hingga dua belas periode ke depan bernilai cukup baik untuk model ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2), yaitu 16.88%. Selain itu, berdasarkan nilai MAD dan MSE terlihat juga bahwa model ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2) lebih baik daripada model ARIMA(2,0,0)(1,1,0)12 ‒ GARCH(1,2). Dengan demikian, dapat disimpulkan bahwa model ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2) merupakan model prediksi luas panen padi nasional yang sesuai dengan hasil prediksi yang cukup baik. Kata kunci: data deret waktu, luas panen, GARCH, asimetri

SUMMARY TEUKU ACHMAD IQBAL. Modelling National Area Harvested of Paddy Using GARCH Methods (Generalized Autoregressive Conditional Heteroscedastic) Model. Supervised by KUSMAN SADIK and I MADE SUMERTAJAYA. Time series data of national area harvested of paddy has high volatility and non homogeneous variance. A time series data with non homogeneous variance at time is called time series data with conditional heteroscedasticity. Time series analysis methods that can be used to overcome heteroskedasticity are GARCH models. However, the data of area harvested of paddy contained in the possibility of asymmetric volatility. To overcome the influence of asymmetry, some asymmetry GARCH models can be used, such as: exponential asymmetric GARCH model (EGARCH), quadratic asymmetric GARCH model (QGARCH), TGARCH model, and non-linear asymmetry GARCH model (NAGARCH). This study aims to model the national area harvested of paddy by incorporating elements of varians heterogeneity and the influence of asymmetry on its data using five types of symmetry, asymmetry, and non-linear GARCH models, and find the best models of those five types of GARCH models. Model for national area harvested of paddy are ARIMA(2,0,0)(1,1,0)12 ‒ GARCH(1,2) and ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2). Based on the mean absolute percentage error (MAPE) value to twenty two periods ahead, ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2) better than ARIMA(2,0,0)(1,1,0)12 ‒ GARCH(1,2) but MAPE values for both models is quite high because there is some predicted value deviates quite far from the actual value. However, the value of MAPE to twelve periods ahead is low, 16.88% for ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2). Furthermore, based on the value of mean absolute deviation (MAD) and mean square error (MSE), ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2) also seems to be the better model than ARIMA(2,0,0)(1,1,0)12 ‒ GARCH(1,2). Thus, it can be concluded that the quadratic GARCH model is a fit model of national area harvested of paddy with a fairly good prediction results. Keywords: time series data, area harvested, GARCH, asymmetric

© Hak Cipta Milik IPB, Tahun 2014 Hak Cipta Dilindungi Undang-Undang Dilarang mengutip sebagian atau seluruh karya tulis ini tanpa mencantumkan atau menyebutkan sumbernya. Pengutipan hanya untuk kepentingan pendidikan, penelitian, penulisan karya ilmiah, penyusunan laporan, penulisan kritik, atau tinjauan suatu masalah; dan pengutipan tersebut tidak merugikan kepentingan IPB Dilarang mengumumkan dan memperbanyak sebagian atau seluruh karya tulis ini dalam bentuk apa pun tanpa izin IPB

PEMODELAN LUAS PANEN PADI NASIONAL DENGAN METODE GARCH (Generalized Autoregressive Conditional Heteroscedastic)

TEUKU ACHMAD IQBAL

Tesis sebagai salah satu syarat untuk memperoleh gelar Magister Sains pada Program Studi Statistika Terapan

SEKOLAH PASCASARJANA INSTITUT PERTANIAN BOGOR BOGOR 2014

Judul Tesis Nama NIM

: Pemodelan Luas Panen Padi Nasional dengan Metode GARCH (Generalized Autoregressive Conditional Heteroscedastic) : Teuku Achmad Iqbal : G152100061

Disetujui oleh Komisi Pembimbing

Dr Ir Kusman Sadik, Msi Ketua

Dr Ir I Made Sumertajaya, MSi Anggota

Diketahui oleh

Ketua Program Studi Statistika Terapan

Dekan Sekolah Pascasarjana

Dr Ir Anik Djuraidah, MS

Dr Ir Dahrul Syah, MScAgr

Tanggal Ujian: 9 Januari 2014

Tanggal Lulus:

PRAKATA Puji dan syukur penulis panjatkan kepada Allah subhanahu wa ta’ala atas segala karunia-Nya sehingga karya ilmiah ini berhasil diselesaikan. Tema yang dipilih dalam penelitian yang dilaksanakan sejak bulan Maret 2012 ini ialah luas panen padi, dengan judul Pemodelan Luas Panen Padi Nasional dengan Metode GARCH (Generalized Autoregressive Conditional Heteroscedastic). Terima kasih penulis ucapkan kepada Bapak Dr Ir Kusman Sadik, MSi dan Bapak Dr Ir I Made Sumertajaya, MSi selaku pembimbing. Terima kasih diucapkan juga kepada Bapak Dr Farit Mochamad Afendi, MSi selaku penguji dan Dr Ir Anik Djuraidah, MS selaku penguji dan ketua program studi Statistika Terapan. Di samping itu, penghargaan penulis sampaikan kepada Sub Bagian Data dan Informasi, Sekretariat Direktorat Jenderal Tanaman Pangan, Kementerian Pertanian yang telah membantu selama pengumpulan data. Ungkapan terima kasih juga disampaikan kepada istri, mama, ayah, serta seluruh keluarga, atas segala doa dan kasih sayangnya. Karya ilmiah ini akan diterbitkan di Jurnal Penelitian Pertanian Tanaman Pangan pada tahun 2014. Semoga karya ilmiah ini bermanfaat.

Jakarta, Januari 2014 Teuku Achmad Iqbal

DAFTAR ISI

DAFTAR TABEL

v

DAFTAR GAMBAR

v

DAFTAR LAMPIRAN

v

PENDAHULUAN Latar Belakang Tujuan Penelitian Manfaat Penelitian

1 1 2 2

METODOLOGI Data Metode Analisis

2 2 3

HASIL DAN PEMBAHASAN Deskripsi Data Pembangunan Model Rataan Pembangunan Model Ragam Pemeriksaan Model Ragam Prediksi dan Validasi Penerapan Model

12 12 12 15 20 21 22

SIMPULAN DAN SARAN Simpulan Saran

23 23 23

DAFTAR PUSTAKA

24

LAMPIRAN

25

RIWAYAT HIDUP

28

DAFTAR TABEL 1 2 3 4 5 6 7 8 9 10

Uji ADF data luas panen padi Ringkasan hasil pendugaan parameter model ARIMA Hasil uji Ljung-Box model tentatif Hasil uji langrange multiplier (LM) hingga lag 12 Pendugaan Parameter Model GARCH (1,2) Pendugaan Parameter Model EGARCH (1,1) Pendugaan Parameter Model QGARCH (1,2) Pendugaan Parameter Model TGARCH (1,1) Pendugaan Parameter Model NAGARCH (1,1) Uji kehomogenan ragam galat baku pada model GARCH, QGARCH, dan TGARCH 11 Ringkasan hasil validasi dua puluh dua periode kedepan

13 14 15 16 17 18 18 19 20 20 22

DAFTAR GAMBAR 1 Skema dari metode analisis 2 Plot data luas panen bulanan padi nasional periode Januari 2000 sampai dengan Februari 2012 3 Plot ACF data luas panen padi nasional setelah dilakukan pembedaan terhadap musiman 4 Plot PACF data luas panen padi nasional setelah dilakukan pembedaan terhadap musiman 5 Plot sisaan data series luas panen padi nasional 6 Prediksi dan validasi model GARCH dan QGARCH 7 Penerapan model ARIMA(2,0,0)(1,1,0)12 – QGARCH(1,2) untuk prediksi luas panen padi nasional periode Januari 2014 hingga Desember 2014

11 12 13 14 16 21 22

DAFTAR LAMPIRAN 1 Plot data bulanan luas panen padi nasional setelah dilakukan pembedaan terhadap musiman 2 Plot residual untuk model tentatif ARIMA (2,0,0)(1,1,0)12 3 Hasil pendugaan parameter model-model ARIMA overfitting 4 Pemilihan model GARCH 5 Pemilihan model EGARCH 6 Pemilihan model QGARCH 7 Pemilihan model TGARCH 8 Pemilihan model NAGARCH

25 25 26 26 26 27 27 27

1 PENDAHULUAN Latar Belakang Padi merupakan salah satu komoditas hasil pertanian tanaman pangan yang sangat strategis di Indonesia. Semua kebijakan pemerintah terkait komoditas ini berdampak luas, tidak hanya secara sosial dan ekonomi, tetapi juga politik (BPSRI 2012). Karena itu, pengambilan kebijakan pada komoditas padi perlu didukung dengan data yang lengkap, akurat, dan terkini agar kebijakan tersebut lebih fokus dan tepat sasaran. Salah satu informasi penting sebagai dasar pengambilan kebijakan terkait komoditas padi adalah data luas panen. Dalam sepuluh tahun terakhir, Indonesia mengalami peningkatan luas panen padi yang mencapai 1,96 juta ha. Peningkatan luas panen tersebut berfluktuasi di beberapa periode dengan volatilitas yang tinggi. Hal ini ditunjukkan oleh suatu fase di mana fluktuasinya relatif tinggi dan kemudian diikuti fluktuasi yang relatif rendah dan kembali tinggi, seperti yang terjadi pada periode tahun 2000 – 2011 (Ditjen Tanaman Pangan 2012). Peramalan yang dapat mengakomodir pengaruh volatilitas pada data luas panen padi tersebut dapat digunakan sebagai masukan dalam pengambilan kebijakan. Model-model deret waktu telah banyak digunakan untuk tujuan peramalan. Permasalahan volatilitas data deret waktu menyebabkan asumsi ragam satu periode proyeksi konstan pada model deret waktu tradisional tidak terpenuhi atau ragam sisaan menjadi tidak konstan (heteroskedastisitas). Oleh karena itu, Engle (1982) memperkenalkan proses stokastik yang disebut autoregressive conditional heteroscedastic (ARCH) model. Namun, seringkali pada saat sedang menentukan model ARCH dibutuhkan orde yang besar agar didapatkan model yang tepat. Karenanya, Bollerslev (1986) mengembangkan model ARCH ke dalam model GARCH untuk menghindari orde ARCH yang besar. Kedua model tersebut telah terbukti bermanfaat untuk pemodelan berbagai fenomena deret waktu karena banyak peubah-peubah deret waktu menunjukkan adanya autokorelasi dan heterokedastik yang dinamik. Akan tetapi, pada data luas panen padi terdapat kemungkinan asimetris dalam volatilitasnya (Ditjen Tanaman Pangan 2012). Ramirez dan Shonkwiler (2001) mengemukakan bahwa model GARCH sisaan simetri cenderung mengecilkan nilai sisaan dari model rataan ketika distribusi sisaan yang sebenarnya adalah asimetri. Untuk itu beberapa model GARCH sisaan asimetri telah dikembangkan untuk mengatasi permasalahan tersebut, antara lain: model GARCH sisaan eksponensial asimetris (EGARCH) yang dikembangkan oleh Nelson (1991), model GARCH sisaan kuadratik asimetri model (QGARCH) oleh Sentana (1995), model threshold asimetri GARCH (TGARCH) oleh Zakoian (1994) dan model GARCH sisaan non-linier asimetri (NAGARCH) oleh Engle dan Ng (1993). Model EGARCH menjamin ragam bersyarat yang dihasilkan selalu positif dan bebas dari tanda parameter estimasi dalam model. Model QGARCH memungkinkan gejolak positif dan negatif memiliki dampak yang berbeda dengan periode sebelumnya. Model T-GARCH memiliki kendala positif pada parameter-parameternya yang menjamin ragam

2 bersyarat selalu bernilai positif. Model NAGARCH mampu mengukur efek pengungkit dan efek ukuran sampel. Beberapa penelitian bidang terapan juga telah menggunakan model-model GARCH dengan hipotesis asimetri tersebut, antara lain: Zheng et al. (2008) dalam penelitian pasar makanan di AS dan Rezitis dan Stavropoulos (2007a, b) pada penelitian industri ayam pedaging dan daging domba di Yunani GARCH dengan hipotesis volatilitas yang asimetris. Karena adanya kemungkinan asimetri pada data luas panen padi nasional, perlu dilakukan penelitian dengan menggunakan model-model GARCH sisaan simetris, sisaan asimetris, dan sisaan non-linear untuk mendapatkan prediksi luas panen padi nasional yang tepat dan akurat. Model-model tersebut dibandingkan dengan cara diuji dan dievaluasi, kemudian yang paling tepat akan dipilih untuk menggambarkan volatilitas luas panen padi nasional dan untuk mendapatkan persamaan yang dapat memprediksi luas panen padi nasional.

Tujuan Penelitian 1. Memodelkan luas panen padi nasional dengan cara memasukkan unsur keheterogenan ragam dan pengaruh keasimetrikan pada data luas panen padi nasional, dengan menggunakan lima jenis model GARCH sisaan simetri, asimetri, dan non linier. 2. Menentukan model terbaik dari data luas panen padi nasional.

Manfaat Penelitian Penelitian ini dapat digunakan untuk menghasilkan prediksi luas panen padi nasional berdasarkan pemodelan data deret waktu.

2 METODOLOGI Data Penelitian ini menggunakan data sekunder luas panen padi nasional dengan satuan hektar. Data yang digunakan merupakan data bulanan yang diambil dari bulan Januari tahun 2000 sampai dengan bulan Desember tahun 2013. Data diperoleh dari Direktorat Jenderal Tanaman Pangan, Kementerian Pertanian dan BPS-RI. Data primer luas panen padi dikumpulkan setiap bulan oleh mantri tani atau petugas dari Dinas Pertanian Kabupaten/Kota yang ada di setiap kecamatan. Kemudian data tersebut direkapitulasi secara berjejang, mulai dari level kabupaten, level provinsi hingga akhirnya direkapitulasi secara Nasional di Direktorat Jenderal Tanaman Pangan, Kementerian Pertanian. Pengumpulan data Statistik Pertanian (SP) tanaman pangan, termasuk padi, dilakukan secara lengkap melalui pendekatan area di seluruh kecamatan. Data luas panen padi diperoleh dengan cara penaksiran menggunakan metode pandangan mata (Eye estimate).

3 Metode ini dilakukan dengan cara perkiraan berdasarkan pencatatan yang dilakukan oleh pegawai/petugas desa, dengan syarat bahwa luas baku lahan telah diketahui terlebih dahulu dan yang melakukan taksiran sudah berpengalaman.

Metode Analisis Tahapan analisis dalam penelitian ini adalah: 1 Melakukan analisis data secara deskriptif dengan cara membuat plot data untuk mempelajari karakteristiknya. 2 Membangun model rataan yang berupa model Box-Jenkins. Model Box-Jenkins merupakan salah satu teknik prediksi model deret waktu yang hanya berdasarkan perilaku data peubah yang diamati. Model Box-Jenkins ini secara teknis dikenal sebagai model autoregressive integrated moving average (ARIMA). Model ini berbeda dengan model struktural baik model kausal maupun simultan di mana persamaan model tersebut menunjukkan hubungan antara beberapa peubah. Model Box-Jenkins ini terdiri dari beberapa model, yaitu: autoregressive (AR), moving average (MA), autoregressive-moving average (ARMA), dan autoregressive integrated moving average (ARIMA). a Identifikasi model rataan Sebelum menentukan model ARIMA tentatif, perlu dilakukan pengujian kestasioneran terhadap rataan. Pemerikasaan kestasioneran terhadap rataan secara deskriptif dilakukan dengan menggunakan plot autocorrelation function (ACF) dan partial autocorrelation function (PACF). Kemudian pemeriksaan dilanjutkan menggunakan uji Augmented Dickey Fuller (ADF) yang merupakan uji formal yang digunakan untuk melihat kestasioneran dari set data. Uji tersebut merupakan pengembangan dari uji Dickey Fuller (Enders 2004). Uji ADF menggunakan proses higher order autoregressive untuk peubah terikat. Proses ini memungkinkan pengujian pada ordo tinggi. Misal persamaan autoregressive ordo ke – p : Y𝑡 = ∅1Y𝑡 – 1 + ∅2Y𝑡 – 2 + ⋯ + ∅𝑝Y𝑡−𝑝 + u𝑡 Pendekatan ADF mengontrol korelasi ordo lebih tinggi dengan menambahkan lag periode pembedaan dari peubah terikat Y terhadap terhadap sisi kanan persamaan sehingga diperoleh: ΔY𝑡 = γY𝑡 – 1 + ∑ + u𝑡 dengan ∑ ∅ ) dan ∑ ∅ ( Hipotesis yang digunakan untuk uji ADF adalah : H0 : γ = 0 (Data belum stasioner dalam rataan) H1 : γ 0 (Data sudah stasioner dalam rataan) dengan statistik uji : di mana n adalah banyaknya amatan yang digunakan. Hipotesis nol ditolak jika statistik uji ADF ( ) lebih kecil dari nilai kritis Dickey-Fuller pada taraf nyata tertentu. Dengan demikian data dapat dikatakan sudah stasioner dalam

4 rataan (Hamilton 1994). Selanjutnya, berdasarkan ACF dan PACF ditentukan model ARIMA tentatif. b Pendugaan parameter model rataan Setelah berhasil identifikasi model ARIMA tentatif selanjutnya dilakukan pendugaan parameter model. Model rataan yang memiliki penduga parameter yang nyata dipilih sebagai model tentatif. c Pemeriksaan model rataan i Mempelajari secara deskriptif nilai sisaan Nilai sisaan dipelajari secara deskriptif untuk melihat apakah masih terdapat beberapa pola yang belum diperhitungkan. Selanjutnya, dilakukan pemeriksaan kebebasan pada sisaan (tidak autokorelasi) menggunakan Uji Ljung-Box. Statistik uji Ljung-Box dinyatakan sebagai berikut (Enders 2004): 𝑟 ∑ dengan 𝑟 adalah autokorelasi sisaan ke–j, n adalah banyaknya pengamatan, dan k adalah lag maksimum yang diinginkan. Hipotesis yang akan diuji adalah: H0 : Tidak terdapat autokorelasi antar sisaan di semua lag k H1 : Terdapat autokorelasi antar sisaan di semua lag k Statistik uji Ljung-Box menyebar Khi-kuadrat dengan derajat bebas k-p-q, di mana p dan q merupakan orde pada model. Jika nilai QLB > maka hipotesis nol (H0) ditolak dan artinya model yang dibangun tidak layak (Cryer 2008). ii Mendeteksi adanya ketidakhomogenan ragam sisaan pada model rataan Langkah sederhana untuk pemeriksaan ini adalah melalui plot deret waktu data sisaan. Selanjutnya, dilakukan pengujian keheterogenan ragam bersyarat untuk mendeteksi keberadaan proses ARCH/GARCH dengan menggunakan uji langrange multiplier (LM). Sisaan yang diperoleh dari model ARIMA dikuadratkan. Kemudian dilanjutkan dengan meregresikan kuadrat sisaan dengan menggunakan konstanta sampai lag ke q, sehingga membentuk persamaan regresi sebagai berikut: ⋯ Jika nilai dugaan sampai dengan bernilai nol, maka dapat disimpulkan bahwa tidak memiliki autokorelasi yang nyata atau dengan kata lain tidak terdapat pengaruh ARCH. Sehingga hipotesis yang digunakan dalam pengujian ini adalah: H0 : ⋯ (Tidak ada pengaruh ARCH/GARCH) H1 : minimal ada satu , untuk i = 1,...,q (Ada pengaruh ARCH/GARCH) dengan statistik uji LM sebagai berikut : LM = nR2 di mana n merupakan jumlah amatan dan R2 merupakan koefisien determinasi dari model regresi kuadrat sisaan diatas. Statistik uji LM ini mengikuti sebaran khi-kuadrat dengan derajat bebas q yang merupakan ordo dari ARCH. Hipotesis nol (H0) akan ditolak jika statistik uji LM lebih besar dari nilai tabel dengan taraf nyata tertentu.

5 iii Pemeriksaan kemungkinan adanya asimetri dalam model ragam Pemeriksaan kemungkinan adanya asimetri dalam model ragam dilakukan dengan melakukan pendugaan parameter empat jenis model GARCH asimetri dan Non-linier. Pendugaan parameter dilakukan dengan menggunakan metode kemungkinan maksimum. 3 Pembangunan model ragam Model ragam dapat dibangun apabila terdapat ketidakhomogenan ragam sisaan atau heteroskedastisitas pada model rataan. Model analisis deret waktu yang memperbolehkan adanya heteroskedastisitas adalah model ARCH yang diperkenalkan pertama kali oleh Engle (1982). Model ARCH dipakai untuk memodelkan ragam sisaan yang tergantung pada kuadrat sisaan pada periode sebelumnya secara autoregresi (regresi diri sendiri), atau dengan kata lain model ini digunakan untuk memodelkan ragam bersyarat. Pada pemodelan ARCH, ada dua model yang disusun, yaitu model rataan dan model ragam. Model rataan disusun berdasarkan identifikasi awal. Bentuk model rataan dapat saja berupa model regresi, model ARIMA, konstanta, dan sebagainya. Model ragam menyatakan hubungan antara ragam sisaan pada waktu t dengan besarnya kuadrat sisaan pada waktu sebelumnya. Misalkan model rataan adalah ARIMA (p,d,q) sebagai berikut: ∅p(𝐵)(1-B)dYt = 𝜃𝑞(𝐵)ut di mana pada analisis deret waktu ut diasumsikan sebagai white noise, ut ~ N(0,σ2). Karena data deret waktu seringkali bersifat heteroskedastis maka ragam bersyarat akan mengikuti model berikut: ht = k + α1 + ... + αq + vt (1) proses white noise (ut) yang mengikuti persamaan (1) didefinisikan sebagai model ARCH dengan orde-q (ARCH(q)), dengan vt merupakan peubah acak yang independen dan identik dengan rataan nol dan ragam 1, k > 0 dan αi ≥ 0 untuk i = 1,...,q atau vt ~ N(0,1). Bentuk lain dari ARCH(q) adalah: ut = vt√ di mana: ht = k + α1 + ... + αq dengan q>0, k>0 dan αi≥0 untuk i = 1,...,q. Syarat k>0 dan αi≥0 dibutuhkan agar ragam bersyarat ht > 0. Seringkali pada saat sedang menentukan model ARCH, dibutuhkan orde yang besar agar didapatkan model yang tepat untuk data deret waktu. Oleh karena itu, Bollerslev (1986) mengembangkan model ARCH ke dalam model GARCH untuk menghindari orde ARCH yang besar dan memberikan hasil yang lebih praktis (parsimonious) daripada model ARCH. Dalam model GARCH, perubahan ragam bersyaratnya selain dipengaruhi oleh kuadrat sisaan, juga dipengaruhi oleh ragam bersyarat periode sebelumnya. Secara umum ragam sisaan dalam model GARCH(p,q) mengikuti model berikut: ht = k + α1 + ... + αq + β1ht-1 + ... + βpht-p + vt (2) di mana vt ~ N(0,1).

6 Bentuk lain dari GARCH(p,q) adalah: ut = vt√ di mana ht = k + α1 + ... + αq + β1ht-1 + ... + βpht-p dengan q>0, k>0, αi≥0, βj≥0 untuk i = 1,...,q dan j = 1,...,p. Dan seperti pada ARCH, syarat k>0, αi≥0, dan βj≥0 dibutuhkan agar ragam bersyarat ht > 0. Model GARCH dapat dinyatakan dalam bentuk fungsi antara nilai ragam bersyarat dengan nilai kuadrat sisaan waktu-waktu sebelumnya. Dalam hal ini, model GARCH menghasilkan model yang mendefinisikan bahwa ragam bersyarat adalah fungsi dari kuadrat sisaan dari lag time yang sangat panjang sedangkan model ARCH hanya melibatkan fungsi dari kuadrat sisaan pada laglag awal saja. Hal ini dapat dijelaskan melalui ilustrasi pada model GARCH (1,1) berikut: ht = k + α1 + β1ht-1 namun + β1ht-2 ht-1 = k + α1 sehingga ht = k + α1 + β1(k + α1 + β1ht-2) dan ht = k + α1 + β1(k + α1 + β1 (k + α1 + β1ht-3)) dengan demikian ht = k* + α1* + α2* + α3* + ... + αq* Pada pemodelan GARCH klasik diatas, vt positif dan negatif masa lalu memiliki efek yang sama pada volatilitas saat ini. Disamping itu, penduga parameter model ragam (GARCH) mendekati independen terhadap model rataan (ARIMA) pasangannya jika vt memiliki distribusi simetris (misalnya, normal atau distribusi-t) namun jika vt memiliki distribusi miring maka penduga GARCH dan penduga ARIMA berkorelasi. Distribusi miring tersebut terjadi karena adanya kemungkinan asimetri, yaitu berbeda volatilitas dicatat dalam hal penurunan dari kenaikan dengan jumlah yang sama. Dengan demikian, apabila terdapat kemungkinan efek asimetri maka model GARCH klasik tidak dapat menjelaskannya dengan baik. Oleh karena itu, beberapa model GARCH sisaan asimetri lebih tepat untuk digunakan, antara lain: model GARCH sisaan eksponensial asimetris (EGARCH), model GARCH sisaan kuadratik asimetri (QGARCH), model T-GARCH, dan model GARCH sisaan non-linier asimetri (NAGARCH). Model GARCH klasik hanya dapat menjelaskan volatilitas tetapi model GARCH klasik tidak dapat menjelaskan efek pengungkit karena ragam bersyarat hanya merupakan fungsi magnitude dari nilai-nilai masa lalu dan bukan tanda mereka (Nelson 1991). Model EGARCH mengakomodasi adanya gejolak asimetri tersebut. Misalkan model rataan adalah ARIMA (p,d,q) sebagai berikut: ∅p(𝐵)(1-B)dyt = 𝜃𝑞(𝐵)ut Maka spesifikasi untuk ragam bersyarat model EGARCH adalah: ut = vt√ di mana ∑ ∑ (3)

7 dengan 𝜃

| |

| | dan

| |

( ⁄ )

jika vt ~ N(0,1)

Berbeda dengan model GARCH, model EGARCH tidak memiliki pembatasan parameter dalam model. Penggunaan ln menjamin model EGARCH selalu menghasilkan ragam positif bersyarat yang bebas dari tanda parameter estimasi dalam model dan tidak ada pembatasan diperlukan. Hal ini lebih baik karena pembatasan-pembatasan dalam GARCH model kadangkadang membuat masalah ketika parameter estimasi melanggar ketidaksamaan kendala. Selanjutnya, untuk menduga dampak dari efek asimetris pada gejolak volatilitas, Sentana (1995) memperkenalkan model kuadratik asimetris GARCH (QGARCH). Adanya persamaan tambahan γut - 1 memungkinkan gejolak positif dan negatif memiliki dampak yang berbeda dengan periode sebelumnya. Misalkan model rataan adalah ARIMA (p,d,q) sebagai berikut: ∅p(𝐵)(1-B)dyt = 𝜃𝑞(𝐵)ut Maka secara umum ragam bersyarat dalam model QGARCH mengikuti model berikut: ut = vt√ dengan (4) Perbedaan proses QGARCH dengan GARCH yaitu pada persamaan γut-1 yang memperkenalkan asimetri. Model T-GARCH simetris yang diajukan oleh Zakoian (1994) juga dapat menduga dampak dari efek asimetris pada gejolak volatilitas. Misalkan model rataan adalah ARIMA (p,d,q) sebagai berikut: ∅p(𝐵)(1-B)dyt = 𝜃𝑞(𝐵)ut maka ragam bersyarat pada model T-GARCH simetris yang diajukan oleh (Zakoian 1994), dituliskan mengikuti persamaan berikut: ut = vt√ dengan (5) Parameter-parameter ( dalam persamaan ragam bersyarat tergantung pada peubah ambang batas yt sebagai berikut: jika jika di mana St ditentukan oleh peubah ambang batas yt-1 yang dapat dianggap sebagai peubah exogen maupun endogen, dan nilai ambang batas y0 menentukan peluang 𝑝 𝑝 . Dengan asumsi peubah ambang batas adalah bebas terhadap . Dan model non-linear asimetris GARCH (NAGARCH) yang diusulkan oleh Engle dan Ng (1993) mampu mengukur efek pengungkit dan efek ukuran sampel. Misalkan model rataan adalah ARIMA (p,d,q) sebagai berikut: ∅p(𝐵)(1-B)dyt = 𝜃𝑞(𝐵)ut maka secara umum ragam bersyarat dalam model NAGARCH mengikuti model berikut:

8 ut = vt√ dengan ht = k + βht-1 + α(ut-1 + γ√ )2 (6) di mana ht adalah ragam bersyarat pada saat t dan α, β, γ, k adalah parameter yang akan diduga. Parameter model GARCH bisa diduga dengan metode kemungkinan maksimum. Parameter GARCH dapat diduga dengan metode quasi maximum likelihood yang memaksimalkan logaritma fungsi kemungkinan apabila asumsi vt merupakan peubah acak yang independen dan identik dengan rataan nol dan ragam 1 terpenuhi. Fungsi kemungkinannya adalah sebagai berikut: 𝜃

∏

𝜃

𝑝(

√

)

dengan 𝜃 𝜃 𝜃 untuk t ≥ 1, maka ht didefinisikan sebagai berikut: 𝜃

∑

∑

dan quasi maximum likelihood didefinisikan sebagai berikut: yang setara dengan

𝜃̂

𝜃

𝜃̂

𝜃

di mana ∑

𝜃

𝜃

dan

sehingga: 𝜃 𝜃

}

∑{

𝜃

dan 𝜃 ) 𝜃

𝜃 ( dengan wt =

.

Sedangkan, pendugaan parameter pada EGARCH dilakukan dengan metode kemungkinan maksimum. Jika asumsi kenormalan vt terpenuhi, logaritma fungsi kemungkinan dari model EGARCH adalah sebagai berikut : ∑ ∑ Lt = dengan ∑{

(|

|

| |)}

∑

Misalkan β = (α0, α1, ..., αq, ψ1, ..., ψq, β1, ..., βp). Turunan pertama terhadap parameter-parameter EGARCH adalah:

9

( )∑

∑ dengan ∑{

|

|}

∑

di mana |

|

| |

|

| |

|

Berdasarkan teori Bayes, 𝜃| dengan n observasi dituliskan sebagai berikut: 𝜃| |𝜃 𝜃 di mana |𝜃 adalah fungsi kemungkinan. 𝜃 adalah fungsi kepekatan untuk θ. Dengan asumsi bahwa 𝜃 adalah konstan, fungsi kemungkinannya dituliskan sebagai berikut : L(y|θ)

√

𝑝

di mana θ = (k,α,β,γ)’ merupakan parameter-parameter QGARCH. Dengan menggunakan 𝜃 , parameter QGARCH diinferensiakan sebagai nilai harapan yang dituliskan sebagai berikut: 〈𝜃〉

∫𝜃

𝜃|

𝜃

di mana Z = ∫ 𝜃| 𝜃 adalah normalisasi konstan yang irrelevant untuk estimasi Markov Chain Monte Carlo (MCMC). Teknik MCMC memberikan metode untuk estimasi persamaan di atas secara numerik. Prosedur dasar dari metode MCMC adalah sebagai berikut: ambil sambel θ dari distribusi peluang 𝜃| dengan menggunakan teknik Markov Chain. Setelah mengambil sampel beberapa data, nilai ekpektasi dievaluasi nilai rataan dari data sampel θ(i), di mana 〈𝜃〉

∑𝜃

dengan k adalah jumlah sampel. Sisaan dari k independen data adalah proporsional terhadap . Akan tetapi, secara umum data yang dihasilkan oleh √ metode MCMC saling berkorelasi. Dengan demikian, sisaan akan proporsional terhadap √ , di mana

T

adalah autokorelasi waktu diantara data sampel.

Autokorelasi waktu tergantung pada metode MCMC yang digunakan. Sehingga diharapkan untuk memilih metode MCMC yang dapat menghasilkan data dengan T yang kecil. Selanjutnya, parameter T-GARCH diduga dengan metode kemungkinan maksimum. Logaritma fungsi kemungkinan dari model TGARCH adalah:

10 𝜃

∑

√

∑

𝜃̂ 𝜃 di mana θ = (k0, k1, α0, α1, β0, β1)’. Untuk menduga θ, diperlukan nilai ambang batas y0 sehingga fungsi kemungkinan diatas dapat diformulasi. Pendugaan parameter pada NAGARCH dilakukan dengan metode kemungkinan maksimum dengan logaritma fungsi kemungkinan p(Y; ) didefinisikan sebagai berikut: ∑

∑{

(

√

)

}

di mana  = (λ, ht)’ merupakan vektor dari parameter. Logaritma fungsi kemungkinan tersebut dapat diperoleh apabila model GARCH yang digunakan adalah model GARCH tipe Gaussian seperti dalam persamaan berikut: ht - 0.5ht + vt Yt = ut + apabila vt merupakan peubah acak yang independen dan identik dengan rataan nol dan ragam 1 maka ragam bersyarat ht dapat dihitung sebagai berikut:





2

ht = k + a vt -1  c ht -1 + bht-1 4 Pemeriksaan model ragam Setelah didapat model GARCH simetri, asimetri, atau non-linier dengan penduga parameter yang nyata, selanjutnya pemeriksaan model dilakukan dengan melakukan pemeriksaan pada galat baku. Pemeriksaan model yang dilakukan adalah pemeriksaan kehomogenan galat baku dengan menggunakan uji LM. 5 Dari model yang didapat dilakukan simulasi prediksi dengan langkah-langkah sebagai berikut: a Dilakukan prediksi sebanyak tiga belas periode (tiga belas bulan). b Hitung nilai mean percentage absolute error (MAPE), mean absolute deviation (MAD), dan mean square error (MSE) periode prediksi. Perangkat lunak yang digunakan untuk menentukan model ARCH/GARCH dan menguji kemungkinan adanya asimetri dalam perilaku volatilitas luas panen padi nasional dalam penelitian ini adalah SAS 9.1.

11 Gambar 1 menunjukan proses dari metode analisis untuk mendapatkan prediksi luas panen padi nasional. Analisis data secara deskriptif Pembangunan model rataan Pemeriksaan model rataan

Model rataan sudah sesuai ?

Ragam sudah homogen ?

Pemeriksaan kesimetrikan Pembangunan model ragam Pemeriksaan model ragam

Model ragam sudah sesuai ?

Melakukan prediksi Gambar 1. Skema dari metode analisis

12

3 HASIL DAN PEMBAHASAN Deskripsi Data Data bulanan luas panen padi nasional sebanyak 146 pengamatan. Gambar 2 merupakan plot antara luas panen padi nasional dengan waktu. 3000000

Luas Panen (Ha)

2500000 2000000 1500000 1000000 500000 Jan-00 Mei-00 Sep-00 Jan-01 Mei-01 Sep-01 Jan-02 Mei-02 Sep-02 Jan-03 Mei-03 Sep-03 Jan-04 Mei-04 Sep-04 Jan-05 Mei-05 Sep-05 Jan-06 Mei-06 Sep-06 Jan-07 Mei-07 Sep-07 Jan-08 Mei-08 Sep-08 Jan-09 Mei-09 Sep-09 Jan-10 Mei-10 Sep-10 Jan-11 Mei-11 Sep-11 Jan-12

0

Gambar 2. Plot data bulanan luas panen padi nasional periode Januari 2000 hingga Februari 2012 Perkembangan luas panen bulanan padi periode Januari 2000 sampai dengan Februari 2012 menunjukkan pola musiman yang cenderung meningkat. Periode puncak panen sebagian besar terjadi pada bulan Maret, sedangkan periode panen terendah sebagian besar terjadi pada bulan Desember. Luas panen padi tertinggi sebesar 2.41 juta ha terjadi pada periode ke-111, yaitu pada Maret 2009. Sedangkan luas panen padi terendah sebesar 0.33 juta ha terjadi pada periode ke12, yaitu pada Desember 2000. Pembangunan Model Rataan Hal pertama yang perlu diperhatikan adalah bahwa kebanyakan deret waktu bersifat non-stasioner dan bahwa aspek-aspek AR dan MA dari model ARIMA hanya berkenaan dengan deret waktu yang stasioner. Oleh karena itu, sebelum menentukan model tentatif, perlu dilakukan pengujian kestasioneran terhadap ragam dan nilai tengah. Dari pemeriksaan secara deskriptif (Gambar 2) terlihat bahwa fluktuasi luas panen padi nasional periode 1 sampai 146 tidak konstan pada suatu nilai tertentu dan cenderung menunjukkan pola musiman. Selain itu, simpangan lokal data menunjukkan adanya keheterogenan. Hal ini menunjukkan bahwa data belum stasioner terhadap nilai tengah sehingga data series harus dilakukan pembedaan. Hasil uji Augmented Dickey Fuller (ADF) pada Tabel 1 menunjukkan nilai statistik uji ADF ( ) untuk lag 1 nyata pada α = 5% dengan Pr < Rho sebesar 0.0002 untuk zero mean dan 0.0001 untuk single mean dan trend. Maka Hipotesis nol ditolak, yang artinya data sudah stasioner rataan untuk lag 1. Sedangkan nilai

13 statistik uji ADF ( ) untuk lag 12 tidak nyata pada α = 5% dengan Pr < Rho sebesar 0.6936 untuk zero mean, 0.0991 untuk single mean dan 0.9999 untuk trend. Maka Hipotesis nol diterima, yang artinya data belum stasioner rataan untuk lag 12. Dengan demikian pembedaan yang perlu dilakukan adalah pembedaan musiman untuk lag 12. Tabel 1. Uji ADF data luas panen padi Tipe Lag Rho Pr < Rho Zero Mean 1 -25.47 0.0002 12 0.05 0.6936 Single Mean 1 -261.72 0.0001 12 -11.01 0.0991 Trend 1 -264.52 0.0001 12 285.13 0.9999 Setelah dilakukan pembedaan musiman lag 12, data sudah tidak menunjukkan pola musiman (Lampiran 1). Selanjutnya, Pemerikasaan kestasioneran terhadap nilai tengah dilakukan dengan menggunakan plot ACF dan PACF. Gambar 3 dan Gambar 4 menunjukkan bahwa data produksi padi nasional telah stasioner terhadap nilai tengah.

Gambar 3. Plot ACF data luas panen padi nasional setelah dilakukan pembedaan terhadap musiman

14

Gambar 4. Plot PACF data luas panen padi nasional setelah dilakukan pembedaan terhadap musiman Selanjutnya, dapat ditentukan model tentatif sebagai berikut: apabila ACF dianggap cut off maka didapat model ARIMA (2,0,0)(1,1,0)12, apabila PACF dianggap cut off maka didapat model ARIMA (0,0,2)(0,1,1)12, serta model ARIMA (2,0,2)(1,1,1)12. Setelah berhasil menetapkan identifikasi model ARIMA tentatif selanjutnya dilakukan pengukuran kebaikan model dan pendugaan parameter model. Tabel 2 menunjukkan model ARIMA (2,0,2)(1,1,1)12 memiliki nilai AIC terkecil yaitu sebesar 26.63, artinya memiliki ukuran kebaikan model terbaik. Akan tetapi, koefisien AR(1) dan SAR(12) model ARIMA (2,0,2)(1,1,1)12 tidak nyata sehingga model ini tidak dapat digunakan. Disamping itu, koefisien MA(2) model ARIMA (0,0,2)(0,1,1)12 tidak nyata sehingga model ini juga tidak dapat digunakan. Sedangkan model ARIMA (2,0,0)(1,1,0)12 semua koefisiennya nyata, maka selanjutnya dapat dilakukan pemeriksaan model untuk model tentatif ini. Tabel 2. Ringkasan hasil pendugaan parameter model-model ARIMA tentatif AIC No Model ARIMA Tipe Koefisien Nilai-p 1

(0,0,2)(0,1,1)12

26.68

2

(2,0,0)(1,1,0)12

26.89

3

(2,0,2)(1,1,1)12

26.63

MA(1) MA(2) SMA(12) AR(1) AR(2) SAR(12) AR(1) AR(2) SAR(12) MA(1) MA(2) SMA(12)

0.5426 -0.1111 -0.8972 0.3942 -0.4912 -0.5568 -0.1460 -0.5016 0.0011 0.7610 0.4009 -0.9066

0.0000 0.2029 0.0000 0.0000 0.0000 0.0000 0.3929 0.0000 0.9900 0.0001 0.0093 0.0000

15 Langkah pertama pemeriksaan adalah mempelajari nilai sisaan untuk melihat apakah masih terdapat beberapa pola yang belum diperhitungkan. Secara deskriptif, Lampiran 2 menampilkan pola sisaan dari model tentatif ARIMA (2,0,0)(1,1,0)12. Dari Lampiran 2 terlihat bahwa tidak terdapat pola pada sisaan model tentatif ARIMA (2,0,0)(1,1,0)12. Selanjutnya, dilakukan uji modifikasi BoxPierce (Ljung-Box) untuk membuktikan bahwa model tentatif tersebut sudah sesuai. Hasil uji Ljung-Box pada model tentatif ARIMA (2,0,0)(1,1,0)12 diringkas pada Tabel 3 berikut: Tabel 3. Hasil uji Ljung-Box model tentatif Nilai-p Model tentatif Lag 12 Lag 24 Lag 36 ARIMA (2,0,0)(1,1,0)12

0.172

0.079

0.534

Lag 48 0.325

Berdasarkan hasil uji Ljung-Box pada Tabel 2, model tentatif ARIMA (2,0,0)(1,1,0)12 mempunyai p-value > 0.05 pada lag 12 sampai 48, yang artinya memiliki residual yang saling bebas sehingga model tentatif ini merupakan model yang memadai. Setelah didapatkan model tentatif yang memadai, selanjutnya dilakukan overfitting. Model tentatif yang memadai yang telah didapatkan adalah ARIMA (2,0,0)(1,1,0)12 maka overfittingnya adalah ARIMA (3,0,0)(1,1,0)12, ARIMA (2,0,1)(1,1,0)12, ARIMA (2,0,0)(2,1,0)12, ARIMA (2,0,0)(1,1,1)12. Dari ringkasan hasil pendugaan parameter untuk model-model ARIMA tersebut diketahui bahwa pada model-model ARIMA tersebut ada koefisien yang tidak nyata sehingga model-model tersebut tidak dapat digunakan (Lampiran 3). Dengan demikian model ARIMA (2,0,0)(1,1,0)12 dapat ditetapkan sebagai model rataan yang memadai. Model rataan dapat dituliskan dalam persamaan berikut: (1-∅1𝐵-∅2𝐵2)(1-Φ12𝐵12)(1-B12)1Yt = ut Pembangunan Model Ragam Model ragam dapat dibangun apabila terdapat ketidakhomogenan ragam sisaan pada model rataan. Langkah sederhana untuk pemeriksaan ini adalah melalui plot deret waktu data sisaan. Plot sisaan pada Gambar 5 menunjukkan bahwa ragam sisaan tidak homogen, di mana terdapat periode dengan fluktuasi sisaan yang tinggi dan periode dengan fluktuasi sisaan yang rendah.

16

Gambar 5. Plot sisaan data series luas panen padi nasional Selanjutnya, pemeriksaan apakah terdapat proses ARCH pada sisaan dapat dilakukan melalui uji lagrange multiplier (LM). Hasil uji keberadaan pengaruh ARCH menggunakan uji LM pada Tabel 4 menunjukkan bahwa nilai p signifikan pada α = 0.05 untuk ordo 1-12. Maka hipotesis nol (H0) ditolak, artinya ada pengaruh ARCH/GARCH pada galat model rataan. Banyaknya ordo yang signifikan menunjukkan banyaknya ordo ARCH yang diperlukan untuk memodelkan fungsi ragam. Tabel 4. Hasil uji LM hingga lag 12 ARIMA (2,0,0)(1,1,0)12 Lag LM Nilai p 1 72.95 <.0001 2 77.47 <.0001 3 142.41 <.0001 4 146.33 <.0001 5 148.85 <.0001 6 165.24 <.0001 7 165.31 <.0001 8 165.69 <.0001 9 165.86 <.0001 10 166.05 <.0001 11 167.13 <.0001 12 167.93 <.0001 Model ARCH adalah proses short memory yang hanya memasukkan q kuadrat galat yang digunakan untuk menduga perubahan ragam. Sedangkan model GARCH adalah proses long memory yang menggunakan semua kuadrat galat pada waktu sebelumnya untuk menduga ragam saat ini. Berdasarkan uji LM pada Tabel 4, ordo yang panjang hingga ordo 12 ini mengindikasikan adanya proses GARCH.

17 Pendugaan Parameter Model GARCH Model rataan pada model GARCH adalah sebagai berikut: (1-0.50𝐵+0.35𝐵2)(1+0.32𝐵12)(1-B12)1Yt = ut Model ragam yang sesuai adalah Model GARCH (1,2) dengan parameter k, α1, α2, β masing-masing 0.00, 0.63, 0.64 dan -0.96 yang dapat diformulasikan sebagai berikut: ht = 0.636 + 0.64 – 0.96ht-1 Model GARCH (1,2) dipilih karena memiliki parameter yang signifikan pada α = 0.05. Pemilihan model GARCH dapat dilihat pada Lampiran 4.

Peubah ar1 ar2 ar12 arch0 arch1 arch2 garch1

Tabel 5. Pendugaan Parameter Model GARCH (1,2) Dugaan Standard Error Nilai t 0.50 -0.35 -0.32 0.00 0.63 0.64 -0.96

0.08 0.07 0.07 0.00 0.31 0.29 0.02

6.12 -4.81 -4.63 2.14 2.03 2.17 -54.42

Nilai p <0.0001 <0.0001 <0.0001 0.0342 0.0443 0.0314 <0.0001

Pada model rataan, nilai dugaan ∅1 bernilai positif, artinya ∅1 memiliki pengaruh positif terhadap Yt. Sedangkan nilai dugaan ∅2 dan Φ12 bernilai negatif, artinya ∅2 dan Φ12 memiliki pengaruh negatif terhadap Yt. Pada model ragam, nilai dugaan α1 dan α2 bernilai positif, artinya dan memiliki pengaruh positif terhadap ht. Sedangkan nilai dugaan β bernilai negatif, artinya ht-1 memiliki pengaruh negatif terhadap ht. Disamping itu, diketahui juga bahwa parameter ∅1, ∅2, Φ12, dan β signifikan pada α = 0.01, namun parameter k, α1, dan α2 tidak signifikan pada α = 0.01 (Tabel 5). Model GARCH standar mengasumsikan bahwa gejolak terhadap volatilitas adalah simetris. Untuk melihat apakah perilaku volatilitas luas panen padi menunjukkan adanya efek asimetris, pada penelitian ini akan dicoba empat jenis model GARCH sisaan asimetris dan sisaan non-linear. Model-model tersebut antara lain, sebagai berikut: model EGARCH, model QGARCH, model TGARCH dan model NAGARCH. Pendugaan Parameter Model EGARCH Model rataan pada model EGARCH adalah sebagai berikut: (1-0.39𝐵+0.45𝐵2)(1+0.49𝐵12)(1-B12)1Yt = ut Model EGARCH (1,1) merupakan model ragam yang digunakan, dengan parameter k, α, β, θ masing-masing 10.51, 1.62, 0.57 dan -0.02 yang dapat dirumuskan sebagai berikut: log(ht) = 10.51 + 0.57log(ht-1) + 1.62 g(vt-1) | |] g(vt-1) = -0.02 vt-1 + [| | vt-1 = ut-1/√

18 Peubah ar1 ar2 ar12 earch0 earch1 egarch1 theta

Tabel 6. Pendugaan Parameter Model EGARCH (1,1) Dugaan Standard Error Nilai t 0.39 -0.45 -0.49 10.51 1.62 0.57 -0.02

0.11 0.05 0.04 1.83 0.33 0.08 0.09

3.63 -9.25 -13.72 5.75 4.85 7.40 -0.20

Nilai p 0.0004 <0.0001 <0.0001 <0.0001 <0.0001 <0.0001 0.8441

Tabel 6 menunjukkan nilai dugaan parameter k, α, β pada model EGARCH (1,1) sebesar <0.0001 nyata pada α = 0.05. Namun nilai dugaan parameter θ pada model EGARCH (1,1) sebesar 0.8441 tidak nyata pada α = 0.05. Overfitting model EGARCH (1,1), yaitu EGARCH (1,2), EGARCH (2,1), dan EGARCH (2,2) juga menghasilkan nilai dugaan yang tidak nyata (Lampiran 5). Dengan demikian model EGARCH merupakan model ragam yang tidak sesuai. Pendugaan Parameter Model QGARCH Model rataan pada model QGARCH adalah sebagai berikut: (1-0.29𝐵+0.37𝐵2)(1+0.33𝐵12)(1-B12)1Yt = ut Model QGARCH (1,2) merupakan model ragam yang sesuai, dengan parameter k, α1, α2, β, masing-masing 0.00, 1.02, 1.04, -0.97, dan -1.39 yang dapat dirumuskan sebagai berikut: Model QGARCH (1,2) dipilih karena memiliki parameter yang signifikan pada α = 0.05. Pemilihan model QGARCH dapat dilihat pada Lampiran 6.

Peubah ar1 ar2 ar12 arch0 arch1 arch2 garch1 phi

Tabel 7. Pendugaan Parameter Model QGARCH (1,2) Dugaan Standard Error Nilai t 0.29 -0.37 -0.33 0.00 1.02 1.04 -0.97 -1.39

0.05 0.03 0.10 0.00 0.23 0.22 0.01 0.50

5.87 -14.29 -3.27 3.06 4.49 4.75 -90.70 -2.77

Nilai p <0.0001 <0.0001 0.0014 0.0026 <0.0001 <0.0001 <0.0001 0.0064

Pada model rataan, nilai dugaan ∅1 bernilai positif, artinya ∅1 memiliki pengaruh positif terhadap Yt. Sedangkan nilai dugaan ∅2 dan Φ12 bernilai negatif, artinya ∅2 dan Φ12 memiliki pengaruh negatif terhadap Yt. Pada model ragam, nilai dugaan α1 dan α2 bernilai positif, artinya dan memiliki pengaruh positif terhadap ht. Sedangkan nilai dugaan β dan bernilai negatif, artinya ht-1 dan memiliki pengaruh negatif terhadap ht. Disamping itu, diketahui juga bahwa parameter model QGARCH (1,2) signifikan pada α = 0.01 (Tabel 7).

19 Pendugaan Parameter Model TGARCH Model rataan pada model TGARCH adalah sebagai berikut: (1-0.45𝐵+0.43𝐵2)(1+0.66𝐵12)(1-B12)1Yt = ut Model TGARCH (1,1) merupakan model ragam yang sesuai, dengan parameter masing-masing 167.61, 0.16, dan 0.80 yang dapat dirumuskan sebagai berikut: Model TGARCH (1,1) dipilih karena memiliki parameter yang signifikan pada α = 0.05. Pemilihan model TGARCH dapat dilihat pada Lampiran 7.

Peubah

Tabel 8. Pendugaan Parameter Model TGARCH (1,1) Dugaan Standard Error Nilai t

ar1 ar2 ar12 arch0 arch1_plus arch1_minus garch1

0.45 -0.43 -0.66 167.61 0.16 0.48 0.80

0.0005 0.0004 0.0002 1.116 0.001 0.002 0.000

953.86 -1080.40 -2906.10 150.19 200.26 218.11 3145.83

Nilai p <0.0001 <0.0001 <0.0001 <0.0001 <0.0001 <0.0001 <0.0001

Pada model rataan, nilai dugaan ∅1 bernilai positif, artinya ∅1 memiliki pengaruh positif terhadap Yt. Sedangkan nilai dugaan ∅2 dan Φ12 bernilai negatif, artinya ∅2 dan Φ12 memiliki pengaruh negatif terhadap Yt. Pada model ragam, nilai dugaan parameter-parameter model TGARCH (1,2) bernilai positif, artinya dan memiliki pengaruh positif terhadap ht. Disamping itu, diketahui juga bahwa parameter-parameter model TGARCH (1,2) signifikan pada α = 0.01 (Tabel 8). Pendugaan Parameter Model NAGARCH Model rataan pada model NAGARCH adalah sebagai berikut: (1-0.57𝐵+0.52𝐵2)(1+0.45𝐵12)(1-B12)1Yt = ut Model NAGARCH (1,1) merupakan model ragam yang digunakan, dengan parameter k, α, β, γ masing-masing 10.6497, 1.7028, 0.0435 dan 0.4967 yang dapat dirumuskan sebagai berikut: ht = 0.1 + 0.80ht-1 + 0.27(ut-1 - 0.25√ )2 Tabel 9 menunjukkan nilai dugaan k, α, β, dan γ pada model NAGARCH (1,1) tidak nyata.

20 Peubah ar1 ar2 ar12 arch0 arch1 garch1 gamma

Tabel 9. Pendugaan Parameter Model NAGARCH (1,1) Dugaan Standard Error Nilai t Nilai p 0.57 -0.52 -0.45 0.10 0.27 0.80 -0.25

0.08 0.08 0.06 0.00 0.09 0.04 0.21

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

Biased Biased Biased Biased Biased Biased Biased

Overfitting model NAGARCH (1,1), yaitu NAGARCH (1,2), NAGARCH (2,1), dan NAGARCH (2,2) juga menghasilkan nilai dugaan yang tidak nyata (Lampiran 8). Dengan demikian model NAGARCH merupakan model ragam yang tidak memadai. Pemeriksaan Model Ragam Setelah didapatkan model yang memadai, yaitu: model GARCH, model QGARCH, dan model TGARCH. Selanjutnya, pemeriksaan model dilakukan dengan melakukan pemeriksaan pada galat baku. Pemeriksaan model yang dilakukan adalah pemeriksaan kehomogenan galat baku.

Lag 1 2 3 4 5 6 7 8 9 10 11 12

Tabel 10. Uji kehomogenan ragam galat baku pada model GARCH, QGARCH, dan TGARCH GARCH (1,2) QGARCH (1,2) TGARCH (1,1) LM Nilai p LM Nilai p LM Nilai p 0.2399 2.2296 4.5851 4.9455 5.609 5.7394 6.3062 6.6568 6.6588 6.7729 9.5563 9.7184

0.6243 0.328 0.2048 0.2929 0.3461 0.453 0.5045 0.5741 0.6726 0.7467 0.5707 0.6406

0.1067 0.6388 3.2594 3.8497 4.3675 4.4286 5.8121 6.8398 7.3238 7.839 10.7533 11.1295

0.7440 0.7266 0.3533 0.4267 0.4978 0.6189 0.5619 0.554 0.6035 0.6446 0.4642 0.5179

28.8395 37.8264 40.3523 41.9357 42.3127 42.7313 42.7710 42.8969 42.9009 43.0070 43.0086 43.0907

<.0001 <.0001 <.0001 <.0001 <.0001 <.0001 <.0001 <.0001 <.0001 <.0001 <.0001 <.0001

Pengujian kehomogenan ragam galat baku pada model GARCH dan model QGARCH sudah menunjukkan adanya kehomogenan ragam galat baku dengan nilai p yang tidak signifikan pada α = 0.05 seperti yang terlihat pada Tabel 10. Sedangkan pengujian kehomogenan ragam galat baku pada model TGARCH masih menunjukkan adanya keheterogenan ragam galat baku dengan nilai p yang

21 signifikan pada α = 0.05 sehingga model TGARCH merupakan model ragam yang tidak memadai. Prediksi dan Validasi Prediksi seluruh periode model ARIMA (2,0,0)(1,1,0)12 ‒ GARCH(1,2) dan model ARIMA (2,0,0)(1,1,0)12 ‒ QGARCH(1,2) secara deskriptif dapat dilihat pada Gambar 6. Setelah periode Januari 2001, prediksi model ARIMA (2,0,0)(1,1,0)12 ‒ GARCH(1,2) dan ARIMA (2,0,0)(1,1,0)12 ‒ QGARCH(1,2) hingga periode Februari 2012, yang digunakan untuk pembangunan model, menunjukkan pola yang hampir sama dengan nilai aktual. Selanjutnya, periode Maret 2012 hingga Desember 2013, yang digunakan sebagai validasi model, juga menunjukkan pola hampir sama dengan nilai aktual. 3000000

Luas Panen (Ha)

2500000 2000000 1500000 1000000 500000 Jan-00 Mei-00 Sep-00 Jan-01 Mei-01 Sep-01 Jan-02 Mei-02 Sep-02 Jan-03 Mei-03 Sep-03 Jan-04 Mei-04 Sep-04 Jan-05 Mei-05 Sep-05 Jan-06 Mei-06 Sep-06 Jan-07 Mei-07 Sep-07 Jan-08 Mei-08 Sep-08 Jan-09 Mei-09 Sep-09 Jan-10 Mei-10 Sep-10 Jan-11 Mei-11 Sep-11 Jan-12 Mei-12 Sep-12 Jan-13 Mei-13 Sep-13

0

Actual

GARCH

QGARCH

Gambar 6. Prediksi dan validasi model GARCH dan QGARCH Pada periode validasi model, setelah April 2012 luas panen Mei - Agustus 2012 (subround II 2012) dan September - Desember 2012 (subround III 2012) mengalami penurunan. Sedangkan luas panen Januari - April 2013 (subround I 2013) mengalami peningkatan. Berikutnya, luas panen Mei - Agustus 2013 (subround II 2013) dan September - Desember 2013 (subround III 2013) kembali mengalami penurunan. Hal ini terjadi karena periode Januari - April (subround I) merupakan musim panen raya sedangkan periode Mei - Agustus (subround II) dan periode September - Desember (subround III) merupakan musim gadu dan musim paceklik. Berdasarkan nilai MAPE hingga dua puluh dua periode ke depan, model ARIMA (2,0,0)(1,1,0)12 ‒ QGARCH(1,2) lebih baik dibandingkan dengan model ARIMA (2,0,0)(1,1,0)12 ‒ GARCH(1,2) di mana nilai MAPE pada model ARIMA (2,0,0)(1,1,0)12 ‒ QGARCH(1,2) sebesar 25.18% lebih kecil dibandingkan nilai MAPE model ARIMA (2,0,0)(1,1,0)12 ‒ GARCH(1,2) sebesar 27.58% (Tabel 11).

22 Tabel 11. Ringkasan hasil validasi dua puluh dua periode ke depan Model GARCH (1,2) QGARCH (1,2) MAD 334785 274494 MSE 473191 325027 MAPE 27.58% 25.18% Nilai MAPE untuk kedua model tersebut terlihat cukup tinggi, nilai tersebut terjadi karena terdapat beberapa nilai prediksi yang menyimpang jauh dari nilai aktual. Namun nilai MAPE hingga dua belas periode kedepan bernilai cukup baik untuk model ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2), yaitu 16.88%. Disamping itu, berdasarkan nilai MAD dan MSE terlihat juga bahwa model ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2) lebih baik daripada model ARIMA(2,0,0)(1,1,0)12 ‒ GARCH(1,2) (Tabel 11). Dengan demikian, dapat disimpulkan bahwa model GARCH kuadratik merupakan model prediksi luas panen padi nasional dengan hasil prediksi yang cukup baik. Penerapan Model Penerapan model ARIMA (2,0,0)(1,1,0)12 – QGARCH (1,2) untuk memprediksi data luas panen padi nasional periode Januari 2014 sampai dengan Desember 2014 dapat dilihat pada Gambar 7, di mana model rataan ARIMA (2,0,0)(1,1,0)12 dirumuskan sebagai berikut: (1-0.37𝐵+0.37𝐵2)(1+0.30𝐵12)(1-B12)1Yt = ut dan model QGARCH (1,2) yang merupakan model ragam dirumuskan sebagai berikut:

3.000.000

Luas Panen (Ha)

2.500.000 2.000.000 1.500.000 1.000.000 500.000 -

Prediksi

Batas Bawah

Batas Atas

Gambar 7. Penerapan model ARIMA(2,0,0)(1,1,0)12 – QGARCH(1,2) untuk prediksi luas panen padi nasional (ha) periode Januari 2014 hingga Desember 2014

23 Prediksi luas panen padi nasional periode Januari 2014 hingga Desember 2014 mencapai 13.816.278 hektar mengalami peningkatan 41.576 hektar dibandingkan dengan periode Januari 2013 hingga Desember 2013 yang mencapai 13.744.702 hektar. Selain itu, pola luas panen yang ditunjukkan hampir sama, yaitu periode Januari hingga Maret mengalami peningkatan dan setelah periode April mengalami penurunan.

4 SIMPULAN DAN SARAN Simpulan Data luas panen padi nasional memiliki fluktuasi yang sangat besar. Hal tersebut berakibat pada ragam bersyarat yang dimiliki menjadi tidak homogen. Untuk mengatasi masalah tersebut, fungsi rataan dan fungsi ragam akan dimodelkan secara simultan. Selain itu, terdapat pula pengaruh ketidaksimetrikan setelah dilakukan pemeriksaan. Maka pemodelan yang sesuai untuk data luas panen padi nasional adalah model QGARCH. Model QGARCH yang diperoleh adalah ARIMA (2,0,0)(1,1,0)12 ‒ QGARCH(1,2). Validasi model QGARCH(1,2) hingga dua puluh dua periode ke depan menghasilkan nilai MAPE 25.18%. Nilai MAPE tersebut terlihat cukup besar, nilai tersebut terjadi karena terdapat beberapa nilai prediksi yang menyimpang cukup jauh dari nilai aktual. Namun nilai MAPE hingga dua belas periode ke depan bernilai cukup baik, yaitu sebesar 16.88%. Kemudian berdasarkan nilai MAD dan MSE terlihat bahwa model ARIMA(2,0,0)(1,1,0)12 ‒ QGARCH(1,2) sudah baik. Dengan demikian, dapat disimpulkan bahwa model GARCH kuadratik merupakan model prediksi luas panen padi nasional dengan hasil prediksi yang cukup baik. Saran Pemodelan ragam bersyarat dengan model-model GARCH asimetri seperti model QGARCH terus mengalami perkembangan. Penelitian selanjutnya dapat menggunakan modifikasi model-model GARCH asimetri lainnya dengan harapan mendapatkan hasil pemodelan dan prediksi yang lebih baik. Selain itu, perlu dipertimbangkan untuk memasukkan peubah eksogen yang berpengaruh terhadap luas panen padi nasional ke dalam fungsi rataan agar hasil peramalan menjadi lebih baik, seperti ketersediaan benih, serangan organisme penganggu tanaman, dan dampak perubahan iklim.

24

DAFTAR PUSTAKA Agus W. 2011. Ekonometrika Pengantar dan Aplikasinya. Ekonisia. Fakultas Ekonomi UII. Yogyakarta. Bollerslev T. 1986. Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics. 31: 307–327. BPS-RI. 2012. Konversi Gabah Kering Giling ke Beras Tahun 2012. Jakarta. Crayer JD, Chan KS. 2008. Time Series Analysis with Application in R. New York: Springer. Ditjen Tanaman Pangan, Kementerian Pertanian. 2012. Perkembangan Luas Panen, Produktivitas dan Produksi Tanaman Pangan. Jakarta. Enders W. 2004. Applied Econometric Time Series 2nd Edition. New York : John Willey & Sons, Inc. Engle RF. 1982. Autoregressive conditional heteroskedasticity with estimates of the variances of the United Kingdom inflation. Econometrica. 50: 987–1008. Engle RF, David ML, Russell PR. 1987. Estimating time varying risk premium in the term structure: the ARCH-M model. Econometrica. 55: 391–407. Engle RF dan Ng V. 1993. Measuring and testing the impact of news in volatility. Journal of Finance. 48: 1749–1778. Hamilton JD. 1994. Time Series Analysis. New Jersey : Princeton University Press. Makridakis S dan Wheelwright SC. 1989. Forecasting Methods for Management. 5 ed. John Wiley & Sons, New York. Mood G dan Boes DC. 1974. Introduction to the theory of statistics. New York. McGrawhill. Nelson D. 1991. Conditional heteroskedasticity in asset returns: a new approach. Econometrica. 59: 347–370. Ramires OA dan Shonkwiler JS. 2001. Autoregresive conditional heteroscedasticity under error-term non-normality. CASNR. Manu. No. 1-513. Texas Tech University. Lubbock. Rezitis A dan Stavropoulos KS. 2007a. Supply response in the Greek broiler industry: application of GARCH models under rational expectations. Paper presented at the Hellenic Operational Research Society Conference, Arta, Greece, 21–23, June 2007. Rezitis A dan Stavropoulos KS. 2007b. Modeling meat supply response under rational expectations and CAP reforms: application to the Greek sheep industry. Working Paper, University of Ioannina, Agrinio, Greece. Sentana E. 1995. Quadratic ARCH models. Review of Economic Studies. 62: 639661. Taylor S. 1986. Modeling Financial Time Series. Wiley, New York. Zakoian JM. 1994. Threshold heteroskedastic models. Journal of Economic Dynamics and Control. 18: 931-955. Zheng Y, Kinnucan HW, Thompson H. 2008. News and food price volatility. Applied Economics. 40: 1629–1635.

25 Lampiran 1. Plot data bulanan luas panen padi Nasional setelah dilakukan pembedaan musiman

Lampiran 2. Plot residual untuk model tentatif ARIMA (2,0,0)(1,1,0)12

26 Lampiran 3. Hasil pendugaan parameter model-model ARIMA overfitting No Model ARIMA Tipe p-value 1

(3,0,0)(1,1,0)12

2

(2,0,1)(1,1,0)12

3

(2,0,0)(2,1,0)12

4

(2,0,0)(1,1,1)12

AR(1) AR(2) AR(3) SAR(12) AR(1) AR(2) MA(1) SAR(12) AR(1) AR(2) SAR(12) SAR(24) AR(1) AR(2) SAR(12) SMA(12)

Lampiran 4. Pemilihan model GARCH GARCH GARCH GARCH Koefisien (1,1) (1,2) (2,1) arch 0 Biased* 0.0342 Biased* arch 1 Biased* 0.0443 Biased* arch 2 0.0314 <0.0001 garch 1 Biased* Biased* garch 2 Biased* * tidak signifikan pada α = 0.05 Lampiran 5. Pemilihan model EGARCH EGARCH EGARCH EGARCH Koefisien (1,1) (1,2) (2,1) <0.0001 earch 0 Biased* 0.0000** <0.0001 earch 1 Biased* 0.0000** earch 2 Biased* <0.0001 egarch 1 Biased* 0.0000** egarch 2 0.0000** Theta 0.8441 Biased* 0.2545* ** * tidak signifikan pada α = 0.05 ** Menggunakan E-Views 5.0

0.0000 0.0000 0.0986 0.0000 0.2787 0.0000 0.0000 0.0685 0.0000 0.0000 0.0000 0.0272 0.0000 0.0000 0.7484 0.0000

GARCH (2,2) Biased* Biased* Biased* Biased* Biased*

EGARCH (2,2) 0.0000** 0.0000** 0.0098** 0.6715* ** 0.3939* ** 0.2776* **

27 Lampiran 6. Pemilihan model QGARCH QGARCH QGARCH QGARCH Koefisien (1,1) (1,2) (2,1) arch 0 Biased* 0.0026 Biased* <0.0001 arch 1 Biased* Biased* <0.0001 arch 2 <0.0001 garch 1 Biased* Biased* garch 2 Biased* Phi Biased* 0.0064 Biased* * tidak signifikan pada α = 0.05 Lampiran 7. Pemilihan model TGARCH TGARCH TGARCH Koefisien (1,1) (1,2) arch 0 <0.0001 Biased* arch1_plus <0.0001 Biased* arch1_minus <0.0001 Biased* arch2_plus Biased* arch2_minus Biased* garch 1 <0.0001 Biased* garch 2 * tidak signifikan pada α = 0.05

QGARCH (2,2) Biased* Biased* Biased* Biased* Biased* Biased*

TGARCH (2,1) Biased* Biased* Biased* Biased* Biased*

Lampiran 8. Pemilihan model NAGARCH NAGARCH NAGARCH NAGARCH Koefisien (1,1) (1,2) (2,1) arch0 Biased* Biased* Biased* arch1 Biased* Biased* Biased* arch2 Biased* Biased* garch1 Biased* Biased* garch2 Biased* gamma Biased* Biased* Biased* * tidak signifikan pada α = 0.05

TGARCH (2,2) <0.0001 Biased* Biased* Biased*

NAGARCH (2,2) Biased* Biased* Biased* Biased* Biased* Biased*

28

RIWAYAT HIDUP Penulis dilahirkan di Jakarta pada tanggal 30 Agustus 1983 sebagai anak kedua dari pasangan Teuku Syarif dan Kemala Sari. Pendidikan sarjana ditempuh di Program Studi Meteorologi, Fakultas Matematika dan Ilmu Pengetahuan Alam IPB, lulus pada tahun 2006. Kesempatan untuk melanjutkan ke program magister pada program studi Statistika Terapan di perguruan tinggi yang sama diperoleh pada tahun 2010. Penulis bekerja di Balai Penelitian Tanaman Kelapa dan Palma Lain pada tahun 2008 sampai dengan tahun 2012, di Manado. Selanjutnya, tahun 2012 hingga saat ini penulis bekerja di Sekretariat Direktorat Jenderal Tanaman Pangan di Jakarta. Bidang pekerjaan yang menjadi tanggung jawab penulis ialah data dan informasi tanaman pangan.