1. Pendahuluan Industri kelapa sawit di Indonesia dewasa ini telah berkembang pesat sebagai penghasil kelapa sawit terbesar di dunia. Investasi kelapa sawit di Indonesia dewasa ini adalah sebagai hasil ekspor ke negara-negara di Eropa di bidang industri. Hal ini disebabkan, kelapa sawit merupakan tumbuhan industri penting sebagai penghasil minyak nabati dan bahan bakar (biodiesel)[1]. Kelapa sawit merupakan tanaman yang tumbuh liar dan dibudidayakan. Pembudidayaan kelapa sawit memerlukan keadaan lingkungan yang baik, dan dipengaruhi oleh keadaan iklim dan tanah sebagai faktor utama. Tanaman kelapa sawit merupakan jenis tanaman tropis yang pertumbuhannya dipengaruhi oleh iklim dan jenis tanah. Curah hujan yang dibutuhkan tanaman kelapa sawit rata-rata 1.500-4.000 mm per tahun, sedangkan curah hujan optimal untuk tanaman ini adalah 2.000-3.000 mm per tahun, dengan jumlah hari hujan tidak lebih dari 180 hari per tahun. Faktor iklim lain yang mempengaruhi pertumbuhan tanaman kelapa sawit adalah suhu, tinggi tempat, kelembaban dan penyinaran matahari[2]. Upaya pengelolaan sawit dalam rangka peningkatan hasil produksi di kabupaten Morowali, Sulawesi Tengah menjadi upaya mendasar yang dilakukan baik pihak swasta (perusahaan perkebunan sawit) dan pemerintah. Peningkatan hasil produksi tahunan kelapa sawit dipengaruhi oleh beberapa faktor umum dan mendasar seperti faktor iklim yang telah disebutkan diatas, serta perluasan wilayah yang telah menjadi permasalahan pada umumnya[3]. Perkembangan produksi kelapa sawit pada umumnya disebabkan oleh keadaan iklim dan alam yang dirasakan kurang menguntungkan. Produktifitas kelapa sawit di Morowali, Sulawesi Tengah yang berhubungan dengan intensitas curah hujan di wilayah tersebut menjadi ranah penelitian yang dipaparkan sebagai tujuan penelitian ini. Intensitas curah hujan khususnya di wilayah Indonesia di peroleh dengan metode peramalan (forecasting) berdasarkan deret waktu. Dalam hal ini, metode peramalan yang digunakan adalah pemodelan ARIMA (Autoregressive Integrated Moving Average). Pemodelan ini akan digunakan dalam memprediksi atau meramalkan jumlah produksi kelapa sawit secara tahunan. Tujuan dari penelitian ini adalah memprediksi produksi kelapa sawit dengan indikator curah hujan tahunan menggunakan Model ARIMA (Autoregressive Integrated Moving Average) dalam melakukan peramalan curah hujan dan produksi kelapa sawit tahunan. 2. Kajian Pustaka Curah Hujan Karakteristik curah hujan di Indonesia dibedakan menjadi tiga tipe [4], yakni (1)Tipe ekuatorial, berhubungan dengan pergerakan zona konvergensi yaitu pertemuan dua massa udara (angin) yang berasal dari dua belahan dunia dan kemudian bergerak ke atas, seperti ditunjukkan pada Gambar 1.
2
Gambar 1 Curah Hujan Tipe Ekuatorial[4]
Gambar 1 menunjukkan karakteristik curah hujan tipe ekuatorial dimana posisi daerah konvergensi relatif sempit karena hanya berada pada lintang rendah (InterTropical Convergence Zone-ITCZ). (2)Tipe monsun, yaitu dipengaruhi oleh adanya tekanan tinggi dan tekanan rendah di benua Asia dan Australia secara bergantian. Di Indonesia tipe monsun dicirikan dengan perbedaan antara musim hujan dan musim kemarau. Curah hujan tipe monsun disajikan dalam Gambar 2.
Gambar 2 Curah Hujan Tipe Monsun[4]
Gambar 2 (kiri) menunjukan pergerakan udara dari benua Asia menuju benua Australia menyebabkan musim hujan di Indonesia. Gambar 2 (kanan) menunjukan pergerakan udara dari benua Australia menuju benua Asia yang menyebabkan musim kemarau. (3)Tipe lokal, terbentuk oleh naiknya udara yang menuju ke dataran tinggi atau pegunungan karena pemanasan lokal yang intensif. ARIMA (Autoregressive Integrated Moving Average) Memprediksikan kejadian dan kondisi di masa depan disebut dengan prakiraan, dan melakukan suatu tindakan prediksi disebut sebagai peramalan. Tujuan peramalan adalah guna mengurangi resiko pada proses pengambilan keputusan serta mengurangi biaya tak terduga. Time series adalah jenis data yang dikumpulkan menurut urutan waktu dalam suatu rentang waktu tertentu. Jika waktu dipandang bersifat diskrit (waktu dapat dimodelkan bersifat continue), frekuensi pengumpulan selalu sama (equidistant). Dalam kasus diskrit, frrekuensi dapat berupa misalnya detik, menit, jam, hari, minggu, bulanan
3
atau tahunan[5]. Secara umum model ARIMA dirumuskan dengan notasi sebagai berikut[6] : ARIMA(p,d,q) p menunjukkan orde atau derajat Autoregressive (AR), d menunjukkan orde atau derajat differencing (pembedaan), dan q menunjukkan orde atau derajat Moving Average (MA). Klasifikasi model ARIMA (Autoregressive Integrated Moving Average) 1. Autoregressive Model (AR) Model Autoregresive adalah model yang menggambarkan bahwa variabel dependent (terikat) dipengaruhi oleh variabel dependent itu sendiri pada periode-periode dan waktu-waktu sebelumnya. Variabel dependent merupakan variabel terikat yang besarannya tergantung dari besaran variabel independent (bebas). Banyaknya nilai lampau yang digunakan oleh model, yaitu sebanyak p, menentukan tingkat model ini. Apabila hanya digunakan satu lag dependent, maka model ini dinamakan model autoregressive tingkat satu (first-order autoregressive) atau AR(1). Apabila nilai yang digunakan sebanyak p lag dependent, maka model ini dinamakan model autoregressive tingkat p atau AR(p). 2. Moving Average Model (MA) Perbedaan model moving average dengan model autoregressive terletak pada jenis variabel independent. Bila variabel independent pada model autoregressive adalah nilai sebelumnya dari model dependent itu sendiri, maka pada model moving average sebagai variabel independent adalah nilai residual pada periode sebelumnya. Orde dari nilai MA (yang diberi notasi q) ditentukan oleh jumlah periode variabel independent yang masuk dalam model. 3. Model campuran a) Proses ARMA Model yang menggabungkan antara Autoregressive Model (AR) dan Moving Average Model (MA). b) Proses ARIMA Apabila nonstasioneritas ditambahkan pada campuran proses ARMA, maka model umum ARIMA (p,d,q) terpenuhi. Persamaan proses arima adalah Yt = Yt-1 + + + …...………………………………(1) Dimana : = nilai actual pada saat t periode Yt Yt-1 = nilai data lampau = Parameter yang akan diestimasi = Nilai kesalahan pada saat t periode ARIMA (Autoregressive Integrated Moving Average) sering juga disebut metode runtun waktu Box-Jenkins. ARIMA mempunyai tingkat ketepatan yang baik untuk peramalan jangka pendek menggunakan nilai standard error estimate yang paling kecil, sedangkan untuk peramalan jangka panjang ketepatan peramalannya kurang baik. Biasanya akan cenderung flat (mendatar/konstan) untuk periode yang cukup panjang. ARIMA adalah model yang secara penuh mengabaikan independensi variabel dalam membuat peramalan. ARIMA menggunakan nilai masa lalu dan sekarang dari variabel dependen untuk menghasilkan peramalan jangka pendek yang akurat. ARIMA cocok
4
jika observasi dari deret waktu (time series) secara statistik berhubungan satu sama lain (dependen). Model ARIMA secara umum dapat dinotasikan pada persamaan 2 [8]: ARIMA (p,d,q) = dimana : p = orde/derajat autoregressive(AR) d = orde/derajat differencing (pembedaan) q = orde/derajat moving average (MA) = kesalahan peramalan periode t = konstanta = = = Tahapan Metode ARIMA Urutan tahapan penerapan metode ARIMA adalah sebagai berikut: (1) Identifikasi. (2) Penaksiran parameter dan pengujian. (3) Penerapan model untuk peramalan. Gambar 3 memaparkan tahapan penerapan metode ARIMA dengan pendekatan Bob Jenkins.
Gambar 3 Tahap Pendekatan Box Jenkins
5
Mengidentifikasi model ARIMA yang sesuai yaitu merubah data, apabila diperlukan dapat menjadi bentuk yang stasioner dan memutuskan model sementara dengan menganalisis ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function). Data stasioner adalah data yang tidak mengandung trend, nilainya berfluktuasi di sekitar nilai rataan yang konstan, hal ini dapat dilihat melalui nilai autokorelasi (plot ACF) atau merupakan suatu sekumpulan nilai yang tidak mempunyai hubungan dengan nilai lainnya. Nilai terdahulu tidak dapat memprediksi keluarnya nilai yang berikutnya. Atau dengan kata lain data stasioner adalah data acak[9]. Pembedaan atau differencing dari deret berkala adalah : Zt = Yt – Yt-1 Dimana : Yt = nilai actual pada saat t periode Zt = nilai data yang stasioner pada saat t periode Yt-1 = nilai data lampau ACF (Autocorrelation Function) adalah korelasi diantara variabel itu sendiri dengan selang satu atau beberapa periode kebelakang. Sedangkan PACF (Partial Autocorrelation Function) adalah suatu ukuran dari korelasi dua variabel time series stasioner setelah pengaruh dari variabel lainnya dihilangkan. Data yang tidak stasioner biasanya ditransformasi menjadi data stasioner dengan melakukan differencing (menghitung perubahan atau selisih nilai observasi). Nilai selisih yang diperoleh dicek lagi apakah sudah stasioner atau belum. Jika belum stasioner maka dilakukan differencing lagi. Karakteristik dan persamaan yang umum untuk model ARIMA [ARIMA (p,d,q)(P,D,Q)L], yaitu sebagai berikut[10] : 1. Jika ACF terpotong (cuts off) setelah lag 1 atau lag 2; lag musiman tidak signifikan dan PACF perlahan-lahan menghilang (dies down), maka diperoleh model nonseasonal MA (q=1 atau 2). 2. Jika ACF terpotong (cuts off) setelah lag musiman L; lag nonmusiman tidak signifikan dan PACF perlahan-lahan menghilang (dies down), maka diperoleh model seasonal MA (Q=1). 3. Jika ACF terpotong (cuts off) setelah lag musiman L; lag non musiman terpotong (cuts off) setelah lag 1 atau 2, maka diperoleh model nonseasonal-seasonal MA (q=1 atau 2;Q=1). 4. Jika ACF perlahan-lahan menghilang (dies down) dan PACF terpotong (cuts off) setelah lag 1 atau 2; lag musiman tidak signifikan, maka diperoleh model nonseasonal AR (p=1 atau 2). 5. Jika ACF perlahan-lahan menghilang (dies down) dan PACF terpotong (cuts off) setelah lag musiman L; lag nonmusiman tidak signifikan, maka diperoleh model seasonal AR (P=1). 6. Jika ACF perlahan-lahan menghilang (dies down) dan PACF terpotong (cuts off) setelah lag musiman L; lag nonmusiman terpotong (cuts off) setelah lag 1 atau 2, maka diperoleh model nonseasonal-seasonal AR(p=1 atau 2;P=1).
6
7. Jika ACF dan PACF perlahan-lahan menghilang (dies down) maka diperoleh mixed model (ARMA atau ARIMA). Setelah model ditemukan, maka parameter dari model harus diestimasi. Terdapat dua cara mendasar yang dapat digunakan untuk pendugaan terhadap parameterparameter tersebut, yaitu : 1. Trial and Error yaitu dengan menguji beberapa nilai yang berbeda dan memilih diantaranya dengan syarat yang meminimumkan jumlah kuadrat nilai galat (sum square of residuals). 2. Perbaikan secara iteratif yaitu dengan cara memilih taksiran awal dan kemudian membiarkan program computer untuk memperhalus penaksiran tersebut secara iteratif. Metode ini banyak digunakan dan telah tersedia suatu logaritma (proses komputer). Secara mendasar, model sudah memadai apabila residualnya tidak dapat dipergunakan untuk memperbaiki ramalan atau dengan kata lain residualnya bersifat acak. Jika tidak, hal ini mengindikasikan bahwa model yang digunakan belum sesuai dengan data. Tahap berikutnya setalah model ditemukan, maka peramalan satu atau beberapa periode ke depan dapat dilakukan menggunakan R-tool. R-tool merupakan sebuah Perangkat lunak yang digunakan untuk menganalisa data statistik dan presentasi grafik. Berikut adalah contoh kasus yang memberikan diagnosa dalam menentukan jika sebuah deret waktu stasioner.
Gambar 4a Plot Time Series [5]
Berikut adalah contoh plot time series pada Gambar 4a dapat dijelaskan bahwa data relatif stasioner dan tidak mengandung tren.
7
Gambar 4b Plot ACF dan PACF [5]
Hasil identifikasi bentuk plot ACF dan PACF dari data recruit menunjukkan bahwa ACF cenderung dies down (turun cepat) dan PACF cenderung terputus setelah lag 2. Apabila hasil identifikasi plot ACF diatas cenderung tidak dies down atau plot PACF cenderung tidak terputus pada lag 2, maka perlu dilakukan pengidentifikasian lanjutan [11]. 3. Metode Penelitian Penelitian ini dibagi menjadi tiga tahapan yaitu: (1) Tahap penyusunan data awal. (2) Desain dan arsitektural simulasi. (3) Pemodelan dan visualisasi. Tahap penyusunan data bertujuan untuk menentukan data, lokasi dan studi pustaka yang digunakan dalam proses penelitian. Tahap penyusunan data awal berupa data yang diperoleh dari Dinas Perkebunan Sulawesi Tengah. Tahap desain dan arsitektural simulasi terdiri dari proses input data, peramalan curah hujan menggunakan metode Autoregressive Integrated Moving Average (ARIMA). Data dan variabel yang digunakan dalam penelitian ini meliputi: 1.) Data curah hujan tahunan di Kabupaten Morowali dimulai dari periode tahun 1962 sampai dengan periode tahun 2011, 2.) Data produksi kelapa sawit di Kabupaten Morowali dimulai dari periode tahun 1962 sampai dengan periode tahun 2011. Perangkat yang digunakan untuk melakukan peramalan adalah R-tool. Tahap desain dan arsitektural dapat dilihat pada Gambar 5.
8
Gambar 5 Desain Simulasi
Gambar 5 menunjukkan desain simulasi yang dijelaskan sebagai berikut. Pada bagian Data Layer, terdiri dari data curah hujan di Kabupaten Morowali, Sulawesi Tengah periode 1970-2010. Alasan utamanya, seperti yang telah dipaparkan pada latar belakang masalah adalah curah curah merupakan salah satu faktor yang mempengaruhi pertumbuhan tanaman kelapa sawit. Kedua data tersebut sebagai data inputan pada proses Application Layer. Pada bagian Application Layer, dilakukan proses peramalan produksi kelapa sawit menggunakan metode ARIMA, dengan pemrosesan data menggunakan tool R untuk mendapatkan hasil peramalan curah hujan dan produksi kelapa sawit berupa grafik. Penggunaan tool R diharapkan dapat memberikan hasil analisa yang mudah dalam diinterpretasikan berdasarkan data yang diperoleh. Hasil ini dipaparkan melalui proses Visualisation Layer. 4. Hasil dan Pembahasan Sumber data yang digunakan adalah data curah hujan tahunan dan produksi kelapa sawit di kabupaten morowali, sulawesi tengah[12]. Data curah hujan tahunan dan produksi kelapa sawit yang dipakai adalah data sejak tahun 1962-2011. Berdasarkan data curah hujan tahunan di kabupaten morowali, tahun 1991 merupakan tahun dimana jumlah curah hujan paling tinggi, dengan curah hujan total mencapai 5220 mm, sedang curah hujan terendah terjadi pada tahun 2003 dengan total curah hujan mencapai 2115 mm. Produksi kelapa sawit tertinggi adalah pada tahun 2008 dengan total jumlah produksi sebesaar 279.540 kg, sedang yang terendah pada tahun 1990 sebesar 440.328 kg. Produksi kelapa sawit mengalami peningkatan seiring dengan pertumbuhan atau umur kelapa sawit serta perluasan wilayah perkebunan.
9
Berdasarkan data temuan ini Dapat disimpulkan secara cepat bahwa hubungan antara curah hujan tahunan terhadap produksi kelapa sawit tidak memiliki keterkaitan. Dengan demikian, diperlukan suatu metode aplikatif dalam antisipasi perubahan curah hujan tahunan dan peramalan produksi kelapa sawit. Pengolahan data jumlah curah hujan tahunan ini akan dibuat dalam grafik dengan menggunakan Program R. Visualisasi disajikan dalam bentuk grafik atau plot. Data curah hujan merupakan data awal input selain data produksi kelapa sawit yang akan diproses. Data curah hujan yang ada kemudian diproses menjadi grafik atau plot dengan menggunakan R-tool. Berikut adalah kode program untuk membuat plot Jumlah curah hujan di Kabupaten Morowali > Jumlah = ts(Jumlah[,2], start=1962, frequency=1) # Menampilkan data jumlah curah hujan dari kolom 2 dimulai dari tahun 1962 > plot(Jumlah, type="o") # Menampilan grafik atau plot data jumlah curah hujan
Gambar 6 Data Awal Curah Hujan Tahunan
Gambar 6 menunjukkan grafik atau plot besaran Jumlah curah hujan di Kabupaten Morowali Dari tahun 1962-2011. Besaran curah hujan dalam bentuk grafik tersebut menunjukkan bahwa data tidak membentuk pola garis secara acak tetapi data hanya berada disekitar garis lurus atau rata-rata konstan yang disebut bersifat tidak stasioner atau membentuk pola naik atau turun secara teratur (bersifat trend). Karena data tidak stasioner maka diperlukan tahap identifikasi berikutnya yaitu dengan melakukan proses differencing plot Autocorrelation Function (ACF) dan plot Partial Autocorrelation Function (PACF) dari data jumlah curah hujan. Berikut adalah kode program differencing ACF dan PACF pada R-tool > par(mfrow=c(2,1)) # Menampilkan grafik 2 x 1 dalam 1 jendela > acf(diff(diff(Jumlah),1), 48) # Memanggil fungsi dan menampilkan hasil differencing pada acf data Jumlah > pacf(diff(diff(Jumlah),1), 48) # Memanggil fungsi dan menampilkan hasil differencing pada pacf data Jumlah
10
Gambar 7 Data ACF dan PACF
Gambar 7 menunjukan Gambar ACF dan PACF data yang sudah stasioner. Pola ACF cenderung cut off pada lag 0 dan 1 sedangkan pola PACF cenderung dies down lag 1, 2, 3. Berdasarkan Karakteristik dan persamaan yang umum untuk model ARIMA, Penentuan parameter pada model ARIMA dilakukan dengan cara trial and error, maka dapat dilakukan pendugaan terhadap 2 model menggunakan ARIMA (p,d,q) (P,D,Q) yaitu (2,1,1) (1,1,1), (1,1,1) (2,1,0), dan (1,1,1) (0,1,1). Model dugaan ARIMA yang telah ditetapkan diperlukan estimasi untuk menentukan lag yang akan digunakan. Pada tahap estimasi, ketiga model dugaan tersebut kemudian dibandingkan berdasarkan perbandingan kriteria nilai (Akaike Information Criteria) AIC dan nilai likelihoodnya. Model ramalan yang baik adalah jika nilai likelihood yang lebih besar dan nilai AIC yang lebih kecil. AIC dan Log Likelihood adalah indikator untuk memutuskan lag yang digunakan[13]. Tahap penaksiran dan pengujian dari data Jumlah curah hujan di Kabupaten Morowali dengan ARIMA (2,1,1) (1,1,1), (1,1,1) (2,1,0), dan (1,1,1) (0,1,1) yang di terapkan di R-tool adalah sebagai berikut: >Jumlah.fit1=arima(Jumlah,order=c(2,1,1),seasonal=list(order=c(1,1,1),period=1) ) # Memberi nama model pertama yang ditaksir yaitu (2,1,1) (1,1,1) >Jumlah.fit1 # Memanggil fungsi dari model pertama Call: arima(x = Jumlah, order = c(2, 1, 1), seasonal = list(order = c(1, 1, 1), period = 1)) Coefficients: ar1 ar2 ma1 sar1 sma1 -0.0545 -0.1522 -1.0000 0.1796 -1.0000 s.e. 0.7482 0.2176 0.1322 0.7731 0.1322
11
sigma^2 estimated as 466871: log likelihood = -387.98, aic = 787.96 >Jumlah.fit2=arima(Jumlah,order=c(1,1,1),seasonal=list(order=c(2,1,0),period=1) ) # Memberi nama model kedua yang ditaksir yaitu (1,1,1) (2,1,0) >Jumlah.fit2 # Memanggil fungsi dari model kedua Call: arima(x = Jumlah, order = c(1, 1, 1), seasonal = list(order = c(2, 1, 0), period = 1)) Coefficients: ar1 ma1 sar1 sar2 -0.3336 -1.0000 -0.2021 -0.3685 s.e. 0.3281 0.0564 0.3033 0.1436 sigma^2 estimated as 639772: log likelihood = -391.87, aic = 793.74 >Jumlah.fit3=arima(Jumlah,order=c(1,1,1),seasonal=list(order=c(0,1,1),period=1) ) # Memberi nama model ketiga yang ditaksir yaitu (1,1,1) (0,1,1) >Jumlah.fit3 # Memanggil fungsi dari model ketiga Call: arima(x = Jumlah, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 1)) Coefficients: ar1 ma1 sma1 0.1124 -1.000 -1.000 s.e. 0.1476 0.115 0.115 sigma^2 estimated as 480631: log likelihood = -388.46, aic = 784.91 Pada tahap estimasi model dugaan yang baik adalah model ARIMA (1,1,1)(0,1,1), karena memiliki nilai AIC yang lebih kecil, nilai AIC yang lebih kecil dianggap sebagai hasil yang lebih baik. Apabila ingin menggunakan lag dari variabel dalam model, maka panjang distribusi lag yang digunakan adalah yang meminimumkan nilai AIC. Informasi Akaike Kriteria (AIC) adalah cara untuk memilih model dari satu set model yang paling cocok dengan kebenaran dari parameter. Proses diagnosa dari model ARIMA (1,1,1)(0,1,1) yang telah di pilih yaitu dengan menguji distribusi estimasi residual-nya menggunakan uji statistik Ljung-Box adalah tahap selanjutnya, setelah terlebih dahulu menentukan model yang akan digunakan. Kebenaran model dugaan diatas dibuktikan dengan uji Hasil diagnostics check, diagnosa model yang didapat apakah sudah memenuhi asumsi atau belum dapat dilakukan dengan menulis perintah pada R-tool sebagai berikut > tsdiag(Jumlah.fit3, gof.lag=48) # Perintah untuk mendiagnosa model adalah tsdiag. Jumlah.fit3 model arima tepilih dan lag sebesar 48.
12
Gambar 8 Diagnostics Check Curah Hujan
Kebenaran model dugaan diatas dibuktikan dengan uji Hasil diagnostics check pada Gambar 8. Dan disimpulkan bahwa residual model ARIMA (1,1,1)(0,1,1) telah terdistribusi secara random (white noise). Ini ditunjukkan oleh p-value dari uji LjungBox yang semuanya lebih besar dari 5% atau 0,05 (alpha atau tingkat signifikansi pengujian). Setelah tahap diagnosa diatas, tahap berikutnya adalah peramalan dengan menggunakan model yang terpilih, yaitu dengan model ARIMA (1,1,1)(0,1,1). Berikut adalah tahap peramalan pada data Jumlah curah hujan menggunakan R-tool > Jumlah.pr = predict(Jumlah.fit3, n.ahead=4) # Memberi nama peramalan data Jumlah. Proses peramalannya dengan memasukkan model kedua dan lama waktu yang diinginkan. > U = Jumlah.pr$pred + 2*Jumlah.pr$se # Menentukan nilai batas atas dari peramalan. > L = Jumlah.pr$pred - 2*Jumlah.pr$se # Menentukan nilai batas bawah dari peramalan. >ts.plot(Jumlah,Jumlah.pr$pred,col=1:2,type="o",ylim=c(0,5000), xlim=c(2005, 2015)) # Menampilkan grafik peramalan area dengan memberikan warna merah pada garis peramalannya. > lines(U, col="green", lty="dashed") # Menampilkan garis batas atas dari peramalan dengan memberikan warna hijau. > lines(L, col="yellow", lty="dashed") # Menampilkan garis batas atas dari peramalan dengan memberikan warna kuning
13
Gambar 9 Grafik Prediksi Curah Hujan
Gambar 9 menunjukan hasil prediksi curah hujan selama 4 tahun, garis hitam menunjukkan curah hujan pada tahun 2000-2011. Sedangkan, garis merah adalah hasil peramalan selama 4 tahun yaitu 2012-2015. Garis hijau dan kuning menunjukkan batas atas dan batas bawah peramalan. Dapat disimpulkan bahwa curah hujan di Kabupaten Morowali akan mengalami penurunan, ini dapat dilihat dari plot peramalan Curah Hujan pada gambar 9. Adapun hasil peramalan jumlah curah hujan tahun 2012-2014 disajikan pada Tabel 1. Tabel 1 Hasil Peramalan Curah Hujan Tahun 2012 2013 2014 2015
Jumlah 2738 2714 2699 2686
Data awal yang akan digunakan sebagai input selanjutnya adalah data produksi kelapa sawit, Berikut adalah kode program untuk membuat plot Jumlah curah hujan di Kabupaten Morowali > Jumlah = ts(Jumlah[,2], start=1962, frequency=1) # Menampilkan data jumlah produksi dari kolom 2 dimulai dari tahun 1962 > plot(Jumlah, type="o") # Menampilan grafik atau plot data jumlah produksi
Gambar 10 Data Awal Produksi Kelapa Sawit
Gambar 10 adalah grafik atau plot besaran jumlah produksi kelapa sawit di Kabupaten Morowali dari tahun 1962-2011. Data tersebut berguna untuk memprediksi
14
data peramalan yang akan datang. Grafik gambar 10 diatas memperlihatkan pola data produksi kelapa sawit pada tahun 1962-2011 mengandung model ARIMA yang tidak stasioner karena data tidak membentuk pola garis secara acak tetapi data hanya berada disekitar garis lurus atau rata-rata konstan dan bersifat trend. Karena data tidak stasioner maka diperlukan tahap identifikasi berikutnya yaitu dengan melakukan proses differencing plot Autocorrelation Function (ACF) dan plot Partial Autocorrelation Function (PACF) dari data jumlah curah hujan. Ini dilakukan agar dapat mendapatkan model terbaik yang nanti akan di jadikan untuk peramalan. Berikut adalah kode program differencing ACF dan PACF pada R-tool > par(mfrow=c(2,1)) # Menampilkan grafik 2 x 1 dalam 1 jendela > acf(diff(diff(Jumlah),1), 48) # Memanggil fungsi dan menampilkan hasil differencing pada acf data Jumlah > pacf(diff(diff(Jumlah),1), 48) # Memanggil fungsi dan menampilkan hasil differencing pada pacf data Jumlah
Gambar 11 Data ACF dan PACF
Gambar 11 menunjukan data ACF dan PACF yang telah stationer, Pola ACF cenderung cut off lag 0 dan 1 melewati garis moving average dan pola PACF cenderung dies down (lag 1-9), garis pada lag tersebut tidak terpotong melainkan turun secara bertahap dengan berbentuk gelombang sinus. Berdasarkan Karakteristik dan persamaan yang umum untuk model ARIMA, Penentuan parameter pada model ARIMA dilakukan dengan cara trial and error, maka dapat dilakukan pendugaan terhadap 3 model menggunakan ARIMA (p,d,q) (P,D,Q) yaitu (1,1,1) (1,1,0), (2,1,1) (1,1,1), dan (2,1,0) (0,1,1). Tahap selanjutnya merupakan tahap estimasi, pada tahap estimasi ketiga model dugaan tersebut kemudian dibandingkan berdasarkan perbandingan kriteria nilai (Akaike Information Criteria) AIC dan nilai likelihood-nya. Model ramalan yang baik adalah jika nilai likelihood yang lebih besar dan nilai AIC yang lebih kecil. AIC dan Log
15
Likelihood adalah indikator untuk memutuskan lag yang digunakan. Tahap penaksiran dan pengujian dari data Jumlah produksi kelapa sawit di Kabupaten Morowali dengan ARIMA (1,1,1) (1,1,0), (2,1,1) (1,1,1), dan (2,1,0) (0,1,1) yang di terapkan di R-tool adalah sebagai berikut: >Jumlah.fit1=arima(Jumlah,order=c(1,1,1),seasonal=list(order=c(1,1,0),period=1) ) # Memberi nama model pertama yang ditaksir yaitu (1,1,1) (1,1,0) > Jumlah.fit1 # Memanggil fungsi dari model pertama Call: arima(x = Jumlah, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 0), period = 1)) Coefficients: ar1 ma1 sar1 -0.0210 -1.0000 -0.0210 s.e. 0.3749 0.0582 0.3749 sigma^2 estimated as 1.415e+11: log likelihood = -686.31, aic = 1380.62 >Jumlah.fit2=arima(Jumlah,order=c(2,1,1),seasonal=list(order=c(1,1,1),period=1) ) # Memberi nama model kedua yang ditaksir yaitu (2,1,1) (1,1,1) > Jumlah.fit2 # Memanggil fungsi dari model kedua Call: arima(x = Jumlah, order = c(2, 1, 1), seasonal = list(order = c(1, 1, 1), period = 1)) Coefficients: ar1 ar2 ma1 sar1 sma1 0.7674 0.0094 -1.0000 0.0693 -1.0000 s.e. 2.2651 1.6722 0.1602 2.2534 0.1602 sigma^2 estimated as 1.244e+11: log likelihood = -685.45, aic = 1382.9 >Jumlah.fit3=arima(Jumlah,order=c(2,1,1),seasonal=list(order=c(0,1,1),period=1) ) # Memberi nama model ketiga yang ditaksir yaitu (2,1,1) (0,1,1) > Jumlah.fit3 # Memanggil fungsi dari model ketiga Call: arima(x = Jumlah, order = c(2, 1, 0), seasonal = list(order = c(0, 1, 1), period = 1)) Coefficients: ar1 ar2 sma1 -0.0483 -0.0819 -1.0000 s.e. 0.1477 0.1474 0.0598
16
sigma^2 estimated as 1.401e+11: log likelihood = -686.16, aic = 1380.32 Pada tahap estimasi ini model dugaan yang baik adalah model ARIMA (2,1,1)(0,1,1) karena memiliki nilai AIC yang lebih kecil, Nilai AIC yang lebih kecil dianggap sebagai hasil yang lebih baik. Apabila ingin menggunakan lag dari variable dalam model, maka panjang distribusi lag yang digunakan adalah yang meminimumkan nilai AIC. Informasi Akaike Kriteria (AIC) adalah cara untuk memilih model dari satu set model yang paling cocok dengan kebenaran dari parameter. Tahap selanjutnya merupakan cek diagnosa dari model yang telah di pilih, yaitu dengan menguji distribusi estimasi residual-nya menggunakan uji statistik Ljung-Box. Kebenaran model dugaan terpilih diatas dibuktikan dengan uji Hasil diagnostics check, diagnosa model yang didapat apakah sudah memenuhi asumsi atau belum dapat dilakukan dengan menulis perintah pada R-tool sebagai berikut > tsdiag(Jumlah.fit3, gof.lag=48) # Perintah untuk mendiagnosa model adalah tsdiag, Jumlah.fit3 model arima tepilih dan lag sebesar 48.
Gambar 12 Data Diagnostics Check Kelapa Sawit
Berdasarkan Hasil diagnostics check pada Gambar 12 disimpulkan bahwa residual model ARIMA (1,1,1)(0,1,1) telah terdistribusi secara random (white noise), ini ditunjukkan oleh p-value dari uji Ljung-Box yang semuanya lebih besar dari 5% atau 0,05 (alpha atau tingkat signifikansi pengujian). Hasil diagnostics check dikatakan perlu dikarenakan membuktikan kebenaran model dugaan diatas. Setelah tahap diagnosa diatas, tahap berikutnya adalah peramalan dengan menggunakan model yang terpilih,
17
yaitu dengan model ARIMA (2,1,1)(0,1,1). Berikut adalah tahap peramalan pada data Jumlah curah hujan menggunakan R-tool > Jumlah.pr = predict(Jumlah.fit3, n.ahead=4) # Memberi nama peramalan data Jumlah. Proses peramalannya dengan memasukkan model kedua dan lama waktu yang diinginkan. > U = Jumlah.pr$pred + 2*Jumlah.pr$se # Menentukan nilai batas atas dari peramalan. > L = Jumlah.pr$pred - 2*Jumlah.pr$se # Menentukan nilai batas bawah dari peramalan. >ts.plot(Jumlah,Jumlah.pr$pred,col=1:2,type="o",ylim=c(0,300000),xlim=c(200 5, 2015)) # Menampilkan grafik peramalan area dengan memberikan warna merah pada garis peramalannya. > lines(U, col="green", lty="dashed") # Menampilkan garis batas atas dari peramalan dengan memberikan warna hijau. > lines(L, col="yellow", lty="dashed") # Menampilkan garis batas atas dari peramalan dengan memberikan warna kuning
Gambar 13 Grafik Prediksi Curah Hujan
Gambar 13 menjelaskan mengenai hasil prediksi produksi kelapa sawit selama empat tahun. Garis hitam menunjukkan hasil produksi kelapa sawit sebelumnya yaitu pada tahun 2000-2011. Sedangkan Garis merah adalah hasil peramalan selama empat tahun kedepan yaitu tahun 2012-2015. Garis hijau dan kuning menunjukkan batas atas dan batas bawah peramalan. Dari plot peramalan produksi kelapa sawit pada gambar diatas juga dapat kita simpulkan bahwa produksi kelapa sawit di kabupaten Morowali akan mengalami kenaikan. Hasil prediksi jumlah produksi kelapa sawit tahun 2012-2015 disajikan dalam Tabel 2. Tabel 2 Hasil Peramalan Produksi Kelapa Sawit Dalam Kilogram Tahun 2012 2013 2014 2015
Jumlah Produksi 147,349,0 152,804,8 153,350,6 154,026,4
18
Rata-rata hasil peramalan produksi kelapa sawit dan curah hujan tahunan dari tahun 2012 – 2015 disajikan dalam Tabel 3. Tabel 3 Peramalan Produksi dan Curah Hujan 2012 – 2015 Tahun 2012 2013 2014 2015 Rata - rata
Produksi Kelapa Sawit 174.349.0 152.804.8 153.350,6 154.026,4 152.805,0
Curah Hujan 2738 2714 2699 2686 2709,25
5. Simpulan Berdasarkan pembahasan dan hasil peramalan menggunakan pemodelan ARIMA. Dapat disimpulkan bahwa curah hujan tahunan di kabupaten Morowali dari tahun 19622011, dapat memberikan hasil yang terbaik. Dengan hasil peramalan curah hujan antara tahun 2012-2015 mengalami penurunan, sedangkan produksi kelapa sawit di tahun 2012-2015 mengalami kenaikan. Hasil ini memberikan deskripsi bahwa, produksi kelapa sawit juga dipengaruhi oleh curah hujan tahunan. Penurunan curah hujan yang stabil dengan rata-rata mencapai 2079,25 mm per tahun akan memberikan hasil produktivitas kelapa sawit tahunan yang dinamis meningkat. Dengan pencapaian produksi dengan rata-rata produksi mencapai 152,805 kg per tahun. Dikatakan memberikan produktifitas yang dinamis dikarenakan, curah hujan optimal untuk tanaman ini adalah 2.000-3.000 mm. Dengan melihat hasil yang terjadi maka dapat dikatakan bahwa tingkat produktifitas produksi kelapa sawit sangat bergantung pada kedua faktor diantaranya curah hujan dan disertai keadaan iklim dalam hal ini keadaan tanah, dari faktor tersebut inilah dapat dikatakan berpengaruh untuk pertumbuhan kelapa sawit tersebut. Karena dalam proses untuk meningkatkan tingkat produsi tersebut tanpa disertai dari faktor -faktor yang ada, maka tentu produksi dari kelapa sawit tersebut akan mengalami hambatan dalam peningkatan produksi. 6. DAFTAR PUSTAKA [1]
[2]
Djulin, Adimesra. Darwis, Valeriana. Dkk. Prospek Pengembangan Sumber Energy Alternatif (Biofuel) :Fokus Pada Jarak Pagar. Makalah Seminar Hasil Penelitian T.A. 2006. Badan Penelitian Dan Pengembangan Pertanian Departemen Pertanian. 2006. Profil Kelapa Sawit Final. Dirjen PPPHP-Kementerian Pertanian Republik Indonesia. http://pphp.deptan.go.id/xplore/files/PENGOLAHANHASIL/PENGOLAHAN HASIL/8-Profil Usaha/PROFIL INVESTASI
19
[3]
[4]
[5]
[6]
[7] [8]
[9]
[10] [11] [12] [13]
BIOENERGI/Profil Kelapa Sawit Final.pdf. Di unduh pada tanggal 20 – 03 – 2012. Zakaria, Junaiddin. Subsektor Perkebunan Di Kabupaten MOROWALI Provinsi SULAWESI TENGAH. Jurnal Economic Resources, ISSN. 0852-1158, Vol.11 No.31, Juni 2010. Pusat Penelitian dan Publikasi Ilmiah FE-UMI. 2010. Tukidin. Karakter Curah Hujan Di Indonesia. Jurnal Geografi Fakultas Ilmu Sosial Universitas Negeri Semarang. Volume 7 nomor 2 tahun 2010. Fakultas Ilmu Sosial Universitas Negeri Semarang. 2010. Shen, Gang. Load Forecasting, Using Time Series Models. 2003. http://www.ee.iastate.edu/~jdm/ee653/Load_Forecasting.doc. Diakses tanggal 22 -04 - 2012 Makridakis, Spyros G., Wheelwright, Steven C., Hyndman, Rob J, 1998, Forecasting Method and Applications 3rd Edition, New York: John Wiley & Sons. Gaynor, PE and Kirkpatrick RC. 1994. Introduction to Time Series Modelling and Forecasting in Business and Economics. Mc Grow Hill, Singapore. Ramdani ,Ahmad Luky.2011.Penggunaan Model Arima dalam peramalan suhu udara di sekitar Palangkaraya [Skripsi] Departemen Ilmu komputer Fakultas Matematika dan ilmu Pengetahuan alam Institut Pertanian Bogor. http://repository.ipb.ac.id/handle/123456789/48231. Diakses tanggal 02 -05 -2012 Petra Christian University Library - jiunkpe s1 tmi 2007 jiunkpe-ns-s1-200725403020-5535-arima-chapter4.pdf http://digilib.petra.ac.id/viewer.php?page=7&submit.x=14&submit.y=8&submit= next&qual=high&submitval=next&fname=%2Fjiunkpe%2Fs1%2Ftmi%2F2007% 2Fjiunkpe-ns-s1-2007-25403020-5535-arima-chapter4.pdf. Diakses tanggal 20 04 - 2012 http://daps.bps.go.id/file_artikel/77/arima.pdf . Diakses tanggal 17 – 05 -2012 Suhartono, Analisis Data Statistik Dengan R, Bahan Ajar Jurusan Statistika,ITS Surabaya, 2008 CV. Ramayana Rancang Bangun Konsultan. Laporan Akhir. RTSP dan RTJ Lokasi Umpanga dan Bente, Kab. Morowali. 2010. Maravall, Agustín. A Class Of Diagnostics In The ARIMA-Model-Based Decomposition Of A Time Series. Bank of Spain. 2003.
20
