Pendahuluan
Terdapat tiga jenis data menurut waktu, yaitu:
-
Cross-section Data
Terdiri dari beberapa objek data pada satu waktu tertentu. Misalnya data penduduk dan pendapatan perkapita tingkat kabupaten pada tahun 2021.
-
Time-series Data
Merupakan data yang terdiri atas satu objek tetapi meliputi beberapa periode waktu yaitu harian, bulanan, mingguan, tahunan, dll. Misalnya data jumlah penduduk kabupaten A pada tiga tahun terakhir.
-
Panel Data
Merupakan data yang menggabungkan antara data time-series dan data cross-section. Sehingga data panel akan memiliki beberapa objek dan beberapa periode waktu.
| Jenis Data | Objek | Waktu |
|---|---|---|
| Cross-section Data | Banyak | Satu |
| Time-series Data | Satu | Banyak |
| Panel Data | Banyak | Banyak |
Analisis deret waktu (time series analysis) adalah metode statistik yang digunakan untuk menganalisis data yang terkumpul dalam interval waktu tertentu. Data deret waktu terdiri dari observasi yang diambil secara berkala atau berurutan, seperti data harian harga saham, data bulanan penjualan, data suhu harian, atau data tahunan produksi industri. Analisis deret waktu bertujuan untuk memahami pola, tren, dan fluktuasi dalam data seiring waktu.
Analisis deret waktu memiliki banyak aplikasi, termasuk dalam ekonomi, keuangan, meteorologi, ilmu sosial, dan banyak bidang lainnya di mana data berkumpul seiring waktu. Dengan memahami pola dan tren dalam data deret waktu, kita dapat membuat perkiraan yang lebih baik dan mengambil keputusan yang lebih informasi dalam berbagai konteks.
Konsep Penting
Berikut adalah beberapa konsep penting dalam analisis deret waktu:
Stasioneritas: Analisis deret waktu biasanya dimulai dengan memeriksa stasioneritas data. Stasioneritas berarti bahwa statistik dasar seperti rata-rata, varians, dan kovarians tidak berubah seiring waktu. Jika data tidak stasioner, langkah-langkah seperti diferensiasi dapat digunakan untuk membuatnya stasioner.
Pola dan Tren: Analisis deret waktu mencari pola dan tren dalam data. Pola adalah fluktuasi teratur yang dapat berulang dalam data, sedangkan tren adalah perubahan panjang waktu yang menunjukkan peningkatan atau penurunan secara keseluruhan.
Komponen Deret Waktu: Deret waktu dapat dibagi menjadi tiga komponen utama: komponen musiman (seasonal), komponen tren (trend), dan komponen residual (error). Komponen musiman adalah fluktuasi berulang dalam data yang terkait dengan faktor-faktor musiman seperti musim panas dan musim dingin. Komponen tren adalah perubahan panjang waktu dalam data. Komponen residual adalah bagian dari data yang tidak dapat dijelaskan oleh musim atau tren.
Metode Peramalan (Forecasting): Salah satu tujuan utama analisis deret waktu adalah melakukan peramalan, yaitu memprediksi nilai-nilai di masa depan berdasarkan pola dan tren dalam data historis. Metode peramalan melibatkan penggunaan model statistik atau teknik matematika seperti Moving Average, Exponential Smoothing, atau model ARIMA (AutoRegressive Integrated Moving Average).
Evaluasi Model: Setelah membangun model peramalan, perlu dievaluasi untuk memastikan kinerjanya yang baik. Metrik evaluasi seperti Mean Absolute Error (MAE), Mean Squared Error (MSE), dan Root Mean Squared Error (RMSE) digunakan untuk mengukur sejauh mana model cocok dengan data sebenarnya.
Analisis Spektral: Analisis deret waktu juga dapat mencakup analisis spektral, yang digunakan untuk mengidentifikasi frekuensi dominan atau siklus dalam data deret waktu. Ini sering digunakan dalam aplikasi seperti analisis gelombang dan pengolahan sinyal.
Kompleksitas Model: Dalam analisis deret waktu, perlu mempertimbangkan kompleksitas model. Model yang terlalu sederhana mungkin tidak mampu menangkap pola yang kompleks, sedangkan model yang terlalu rumit dapat mengalami overfitting.
Metode Peramalan
Metode Smoothing
Metode smoothing (pemulusan) dalam analisis deret waktu adalah teknik yang digunakan untuk menghaluskan fluktuasi acak atau variabilitas dalam data deret waktu, sehingga pola dan tren yang mendasarinya menjadi lebih jelas. Tujuannya adalah membuat data yang lebih mudah diinterpretasikan dan digunakan untuk peramalan atau analisis lainnya. Berikut adalah beberapa metode smoothing umum dalam analisis deret waktu:
-
Simple Moving Average (SMA):
- SMA adalah metode pemulusan yang paling dasar. Ini melibatkan perhitungan rata-rata dari nilai-nilai dalam jangka waktu tertentu. Contoh termudah adalah rata-rata bulanan penjualan selama 12 bulan untuk menghaluskan fluktuasi musiman. SMA memberikan bobot yang sama pada semua data dalam jangka waktu yang dipilih.
-
Exponential Smoothing (Pemulusan Eksponensial):
- Metode ini memberikan bobot yang berbeda pada nilai dalam deret waktu, dengan bobot yang semakin berkurang seiring berjalannya waktu. Pemulusan eksponensial berguna ketika nilai yang lebih baru dianggap lebih penting daripada yang lebih lama. Ada beberapa variasi dari metode ini, termasuk pemulusan eksponensial sederhana (SES), Holt-Winters (yang mempertimbangkan komponen musiman dan tren), dan metode eksponensial beradaptasi.
-
Double Exponential Smoothing (Holt’s Linear Exponential Smoothing):
- Metode ini adalah ekstensi dari pemulusan eksponensial sederhana yang juga mempertimbangkan tren dalam deret waktu. Ini berguna ketika data memiliki tren yang signifikan yang harus diperhitungkan selain fluktuasi musiman.
-
Triple Exponential Smoothing (Holt-Winters Exponential Smoothing):
- Metode ini adalah ekstensi dari double exponential smoothing yang juga memasukkan komponen musiman. Ini cocok untuk data deret waktu yang memiliki pola musiman yang jelas selain tren.
-
Metode Keseluruhan (Holistic Methods):
- Pendekatan ini mencoba untuk mengkombinasikan beberapa metode smoothing, seperti SMA dan metode eksponensial, untuk memanfaatkan kekuatan masing-masing metode dalam mengatasi berbagai jenis fluktuasi dalam data.
-
Metode Jendela (Window Methods):
- Metode ini melibatkan pemilihan jangka waktu tertentu (jendela) dan mengambil nilai rata-rata atau bobot lainnya di dalam jendela tersebut. Metode ini digunakan untuk menghaluskan fluktuasi harian atau mingguan dalam data.
-
Lowess (Locally Weighted Scatterplot Smoothing):
- Metode ini merupakan pendekatan nonparametrik yang menghaluskan data dengan cara memberikan bobot berdasarkan jarak relatif antara titik data dalam deret waktu. Ini cocok untuk data yang memiliki fluktuasi yang tidak teratur.
Moving Average
Metode Moving Average (Rata-rata Bergerak) adalah salah satu metode smoothing yang umum digunakan dalam analisis deret waktu. Metode ini digunakan untuk menghaluskan fluktuasi acak dalam data deret waktu dan membantu dalam mengidentifikasi tren atau pola yang mendasarinya.
Simple Moving Average
SMA adalah metode pemulusan yang paling sederhana. Ini melibatkan penghitungan rata-rata dari sejumlah nilai dalam deret waktu yang sama.
Contoh sederhana dari SMA adalah jika kita ingin menghaluskan data penjualan bulanan selama setahun, kita dapat menggunakan SMA 12-bulan, di mana kita menjumlahkan penjualan selama 12 bulan terakhir dan kemudian membaginya dengan 12.
Cocok untuk pola data konstan/stasioner.
-
Prinsip dasar:
-
Data smoothing pada periode ke-t merupakan rata-rata dari m buah data dari data periode ke-t hingga ke-(t-m+1)
\[ S_t=\frac {1}{m} \sum_{i=t-m+1}^t X_i \]
-
Data smoothing pada periode ke-t berperan sebagai nilai forecasting pada periode ke-t+1.
\[ F_t=S_{t-1} \]
\[ F_{n,h}=S_n \]
MA dengan m yang lebih besar akan menghasilkan pola data yang lebih halus.
-
Double Moving Average
Double Moving Average (DMA), juga dikenal sebagai Moving Average Crossover Method, adalah metode peramalan yang melibatkan penggunaan dua moving average dengan periode yang berbeda untuk menghasilkan sinyal peramalan. Metode ini umumnya digunakan dalam analisis deret waktu untuk mengidentifikasi perubahan tren dalam data.
-
Proses penghalusan dengan rata-rata dilakukan dua kali.
-
Tahap I:
\[ S_{1,t}=\frac {1}{m} \sum_{i=t-m+1}^t X_i \]
-
Tahap II:
\[ S_{2,t}=\frac {1}{m} \sum_{i=t-m+1}^t S_{1,i} \]
-
Forecasting dilakukan dengan formula:
\[ F_{2,t,t+h}=A_t+B_t(h) \]
dengan
\[ A_t=2S_{1,t}-S_{2,t} \]
\[ B_t=\frac {2}{m-1}(S_{1,t}-S_{2,t}) \]
-
Metode Moving Average memiliki beberapa keuntungan, seperti:
Penghalusan Data: Metode ini membantu mengurangi fluktuasi acak dalam data, sehingga memudahkan pengamatan pola yang mendasari seperti tren atau musiman.
Ketahanan terhadap Outlier: lebih tahan terhadap outlier daripada beberapa metode lainnya karena mereka mengambil rata-rata sejumlah nilai dalam jangka waktu tertentu.
Sederhana dan Mudah Dimengerti: Metode ini relatif mudah dimengerti dan diimplementasikan.
Namun, metode Moving Average juga memiliki beberapa kelemahan, termasuk:
Lag: Salah satu kelemahan utama dari SMA adalah bahwa ia menghasilkan lag dalam perkiraan. Ini berarti bahwa prediksi SMA selalu tertinggal dibandingkan dengan tren aktual dalam data.
Kurang Responsif terhadap Perubahan: SMA tidak responsif terhadap perubahan yang cepat dalam data karena mendasarkan prediksi pada nilai rata-rata dalam jangka waktu tertentu.
Salah satu kelemahan utama DMA adalah bahwa ia dapat menghasilkan sinyal palsu ketika tren datar atau dalam kondisi pasar yang bergejolak. Selain itu, DMA cenderung memberikan sinyal dengan lag, yang berarti sinyal perubahan tren dapat terlambat.
Exponential Smoothing
Exponential Smoothing (Pemulusan Eksponensial) adalah metode peramalan yang digunakan dalam analisis deret waktu untuk menghaluskan fluktuasi data dan meramalkan nilai masa depan berdasarkan pola-pola dalam data historis. Metode ini memberikan bobot yang berkurang eksponensial seiring waktu kepada nilai dalam deret waktu, dengan memberikan lebih banyak bobot pada data yang lebih baru. Exponential Smoothing berguna untuk memperkirakan nilai masa depan berdasarkan perkiraan sebelumnya dan memberikan penekanan yang lebih besar pada data yang paling baru.
Berikut adalah komponen utama dalam metode Exponential Smoothing:
-
Level (Tingkat):
- Tingkat adalah perkiraan rata-rata atau nilai tengah dari data deret waktu. Ini adalah representasi dari “tingkat” yang mendasari data. Dalam metode Exponential Smoothing, tingkat diperbarui pada setiap periode waktu berdasarkan peramalan sebelumnya.
-
Faktor Penghalusan (Smoothing Factor atau Alpha):
- Faktor penghalusan (α, alpha) adalah angka antara 0 dan 1 yang menentukan sejauh mana bobot diberikan kepada data yang baru dibandingkan dengan data sebelumnya. Semakin besar nilai alpha, semakin besar bobotnya diberikan kepada data yang lebih baru.
-
Periode (Time Period atau P):
- Periode adalah jangka waktu antara dua pengamatan atau periode waktu di dalam data deret waktu. Periode ini digunakan untuk menghitung faktor penghalusan (alpha) yang optimal.
Single Exponential Smoothing
Metode Moving Average mengakomodir pengaruh data beberapa periode sebelumnya melalui pemberian bobot yang sama dalam proses merata-rata. Hal ini berarti bobot pengaruh sekian periode data tersebut dianggap sama. Dalam kenyataannya, bobot pengaruh data yang lebih baru mestinya lebih besar. Adanya perbedaan bobot pengaruh ini diakomodir metode Single Exponential Smoothing (SES) dengan menetapkan bobot secara eskponensial.
Ini adalah bentuk paling dasar dari metode Exponential Smoothing. Ini hanya mempertimbangkan tingkat (level) dalam peramalan.
Nilai smoothing pada periode ke-t:
\[ S_t= \alpha X_t +(1-\alpha)S_{t-1} \]
Nilai \(\alpha\) merupakan parameter pemulusan dengan nilai \(0<\alpha<1\). \(S_1\) biasanya diambil dari rataan beberapa data pertama. Nilai smoothing pada periode ke-t bertindak sebagai nilai forecast pada periode ke-(t+1).
Nilai forecasting diperoleh dengan formula \(F_t=S_{t-1}\) dan \(F_{n,h}=S_n\) .
Double Exponential Smoothing
Metode ini menggabungkan konsep tingkat (level) dan tren dalam peramalan. Ini cocok untuk data yang menunjukkan tren. Double Exponential Smoothing (DES) mirip dengan SES, hanya saja dilakukan dua kali yaitu pertama untuk tahapan ‘level’, dan kedua untuk tahapan ‘tren’.
Nilai smoothing data ke-t:
\[ S_t=L_{t-1}+T_{t-1} \]
\[ T_t=\gamma (L_t-L_{t-1})+(1-\gamma)T_{t-1} \]
\[ L_t=\alpha X_t + (1-\alpha)(L_{t-1}+T_{t-1}) \]
Bila \(Y_t=a+b*t+e\), maka \(L_0=a\) dan \(T_0=b\)
Nilai forecasting diperoleh dengan formula \(F_{t,h}=L_t+h*T_t\)
Di mana Lt adalah peramalan tingkat (level) pada periode saat ini, Tt adalah peramalan tren untuk periode saat ini, Xt adalah data aktual periode saat ini, \(\alpha\) adalah faktor penghalusan untuk tingkat (level), dan \(\gamma\) adalah faktor penghalusan untuk tren.
Metode Winters
Merupakan salah satu pendekatan smoothing untuk data yang berpola musiman (seasonal).
Memiliki dua prosedur penghitungan tergantung kondisi data:
Aditif: komponen musiman bersifat aditif dengan komponen level dan tren
Multiplikatif: komponen musiman bersifat multiplikatif dengan komponen level dan tren.
Winters Aditif
Komponen model:
\[ L_t=\alpha(Y_t-M_{t-p})+(1-\alpha)(L_{t-1}+T_{t-1}) \]
\[ T_t=\gamma(L_t-L_{t-1})+(1-\gamma)T_{t-1} \]
\[ M_t=\delta(Y_t-L_t)+(1-\delta)M_{t-p} \]
Nilai smoothing:
\[ S_t=L_t+T_t+M_{t-p} \]
Nilai Forecast:
\[ F_{t,h}=L_t+h*T_t+M_{t-p+h} \]
Misal, ambil 2q data pertama (q: ordo musiman), lalu hitung rata-rata masing-masing musim.
-
Musim I
\[ V_1=\frac{1}{q}\sum_{i=-2q+1}^{-q} X_i \]
-
Musim II
\[ V_2=\frac{1}{q}\sum_{i=-q+1}^{0} X_i \]
Winter Multiplikatif
Komponen model:
\[ L_t=\alpha(Y_t/M_{t-p})+(1-\alpha)(L_{t-1}+T_{t-1}) \]
\[ T_t=\gamma(L_t-L_{t-1})+(1-\gamma)T_{t-1} \]
\[ M_t=\delta(Y_t/L_t)+(1-\delta)M_{t-p} \]
Nilai smoothing:
\[ S_t=(L_t+T_t)M_{t-p} \]
Nilai Forecast:
\[ F_{t,h}=(L_t+h*T_t)M_{t-p+h} \]
Serupa dengan versi aditif, hanya berbeda dalam menghitung deseasonalized data, di mana:
Accuracy Measures
Beberapa ukuran yang dapat digunakan untuk mengevaluasi seberapa baik metode forecasting data:
Mean Absolute Deviation (MAD)
\[ MAD=\frac{1}{n} \sum_{t=1}^n|X_t-\hat X_t| \]
Mean Squared Deviation (MSD)
\[ MSD=\frac{1}{n} \sum_{t=1}^n(X_t-\hat X_t)^2 \]
Mean Absolute Percentage Error (MAPE)
\[ MAPE= \frac{1}{n} \sum_{t=1}^n |{\frac {X_t-\hat X_t}{X_t}}| \times 100 \]
Ilustrasi 1
Sebagai ilustrasi, tersedia data bagi hasil suatu bank syariah per bulan (File excel Moving Average.csv). Data ini dicatat setiap tanggal 1 di masing-masing bulan. Periodenya dari Januari 1989 hingga Desember 1992, sehingga terdapat 48 pengamatan.
#install.packages("forecast")
#install.packages("TTR")
library(forecast)## Warning: package 'forecast' was built under R version 4.2.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(TTR)## Warning: package 'TTR' was built under R version 4.2.3library(graphics)Data1<-read.csv("Moving Average.csv",header = TRUE,sep = ";")
Data1str(Data1)## 'data.frame': 48 obs. of 2 variables:
## $ t : int 1 2 3 4 5 6 7 8 9 10 ...
## $ xt: int 86 90 86 102 111 92 96 106 102 104 ...#membentuk objek time series
Data1.ts<-ts(Data1$xt)
head(Data1.ts)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] 86 90 86 102 111 92summary(Data1.ts)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 85.0 96.0 106.0 106.9 118.2 130.0#time series plot
ts.plot(Data1.ts, xlab="Waktu", ylab="Xt", main="Time Series Plot")
points(Data1.ts)
SMA m = 3
Smoothing data dilakukan dengan m = 3.
#melakukan Single Moving Average dengan n=3
Data1.sma3<-SMA(Data1.ts, n=3)
ramal.sma3<-c(NA,Data1.sma3)
Datagab.sma3<-cbind(bagihasil_aktual=c(Data1.ts,NA),pemulusan=c(Data1.sma3,NA),ramal.sma3)
Datagab.sma3## bagihasil_aktual pemulusan ramal.sma3
## [1,] 86 NA NA
## [2,] 90 NA NA
## [3,] 86 87.33333 NA
## [4,] 102 92.66667 87.33333
## [5,] 111 99.66667 92.66667
## [6,] 92 101.66667 99.66667
## [7,] 96 99.66667 101.66667
## [8,] 106 98.00000 99.66667
## [9,] 102 101.33333 98.00000
## [10,] 104 104.00000 101.33333
## [11,] 99 101.66667 104.00000
## [12,] 92 98.33333 101.66667
## [13,] 87 92.66667 98.33333
## [14,] 85 88.00000 92.66667
## [15,] 94 88.66667 88.00000
## [16,] 112 97.00000 88.66667
## [17,] 128 111.33333 97.00000
## [18,] 101 113.66667 111.33333
## [19,] 124 117.66667 113.66667
## [20,] 124 116.33333 117.66667
## [21,] 118 122.00000 116.33333
## [22,] 120 120.66667 122.00000
## [23,] 124 120.66667 120.66667
## [24,] 130 124.66667 120.66667
## [25,] 108 120.66667 124.66667
## [26,] 109 115.66667 120.66667
## [27,] 128 115.00000 115.66667
## [28,] 124 120.33333 115.00000
## [29,] 106 119.33333 120.33333
## [30,] 122 117.33333 119.33333
## [31,] 119 115.66667 117.33333
## [32,] 125 122.00000 115.66667
## [33,] 127 123.66667 122.00000
## [34,] 105 119.00000 123.66667
## [35,] 103 111.66667 119.00000
## [36,] 105 104.33333 111.66667
## [37,] 93 100.33333 104.33333
## [38,] 107 101.66667 100.33333
## [39,] 103 101.00000 101.66667
## [40,] 110 106.66667 101.00000
## [41,] 94 102.33333 106.66667
## [42,] 102 102.00000 102.33333
## [43,] 94 96.66667 102.00000
## [44,] 96 97.33333 96.66667
## [45,] 106 98.66667 97.33333
## [46,] 106 102.66667 98.66667
## [47,] 118 110.00000 102.66667
## [48,] 109 111.00000 110.00000
## [49,] NA NA 111.00000#membuat plot
ts.plot (Data1.ts,xlab="periode waktu",ylab="Bagi hasil", col="blue",lty=3)
points(Data1.ts)
lines (Data1.sma3,col="red",lwd=2)
lines (ramal.sma3,col="black",lwd= 2)
title("Rataan bergerak Sederhana n=3",cex.main=1,font.main=4 ,col.main="black")
legend("topleft",legend=c("Data aktual","Pemulusan SMA","Ramalan SMA"),
lty=1:3,col=c ("blue","red","black"), cex=0.7)
Mencari nilai keakuratan
Mean Absolute Percentage Error (MAPE)
\[ MAPE= \frac{1}{n} \sum_{t=1}^n |{\frac {x_t-\hat x_t}{x_t}}| \times 100 \]
# Ukuran Keakuratan
error.sma3<-Data1.ts-ramal.sma3[1:length(Data1.ts)]
## SSE (Sum Square Error)
SSE.sma3 <- sum(error.sma3[4:length(Data1.ts)]^2)
## MSE (Mean Squared Error)
MSE.sma3 <- mean(error.sma3[4:length(Data1.ts)]^2)
## RMSE (Root Mean Square Error)
RMSE.sma3 <- sqrt(mean(error.sma3[4:length(Data1.ts)]^2))
## MAD (Mean Absolute Deviation)
MAD.sma3 <- mean(abs(error.sma3[4:length(Data1.ts)]))
## MAPE (Mean Absolute Percentage Error)
r.error.sma3 <- (error.sma3[4:length(Data1.ts)]/ramal.sma3[4:length(Data1.ts)])*100 # Relative Error
MAPE.sma3 <- mean(abs(r.error.sma3))
performa <- data.frame(
"Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Simple Moving Average N=3" = c(SSE.sma3, MSE.sma3, MAPE.sma3, RMSE.sma3, MAD.sma3))
performaSMA m = 6
Smoothing data dilakukan dengan m = 6.
#melakukan Single Moving Average dengan n=6
Data1.sma6<-SMA(Data1.ts, n=6)
ramal.sma6<-c(NA,Data1.sma6)
Datagab.sma6<-cbind(bagihasil_aktual=c(Data1.ts,NA),pemulusan=c(Data1.sma6,NA),ramal.sma6)
Datagab.sma6## bagihasil_aktual pemulusan ramal.sma6
## [1,] 86 NA NA
## [2,] 90 NA NA
## [3,] 86 NA NA
## [4,] 102 NA NA
## [5,] 111 NA NA
## [6,] 92 94.50000 NA
## [7,] 96 96.16667 94.50000
## [8,] 106 98.83333 96.16667
## [9,] 102 101.50000 98.83333
## [10,] 104 101.83333 101.50000
## [11,] 99 99.83333 101.83333
## [12,] 92 99.83333 99.83333
## [13,] 87 98.33333 99.83333
## [14,] 85 94.83333 98.33333
## [15,] 94 93.50000 94.83333
## [16,] 112 94.83333 93.50000
## [17,] 128 99.66667 94.83333
## [18,] 101 101.16667 99.66667
## [19,] 124 107.33333 101.16667
## [20,] 124 113.83333 107.33333
## [21,] 118 117.83333 113.83333
## [22,] 120 119.16667 117.83333
## [23,] 124 118.50000 119.16667
## [24,] 130 123.33333 118.50000
## [25,] 108 120.66667 123.33333
## [26,] 109 118.16667 120.66667
## [27,] 128 119.83333 118.16667
## [28,] 124 120.50000 119.83333
## [29,] 106 117.50000 120.50000
## [30,] 122 116.16667 117.50000
## [31,] 119 118.00000 116.16667
## [32,] 125 120.66667 118.00000
## [33,] 127 120.50000 120.66667
## [34,] 105 117.33333 120.50000
## [35,] 103 116.83333 117.33333
## [36,] 105 114.00000 116.83333
## [37,] 93 109.66667 114.00000
## [38,] 107 106.66667 109.66667
## [39,] 103 102.66667 106.66667
## [40,] 110 103.50000 102.66667
## [41,] 94 102.00000 103.50000
## [42,] 102 101.50000 102.00000
## [43,] 94 101.66667 101.50000
## [44,] 96 99.83333 101.66667
## [45,] 106 100.33333 99.83333
## [46,] 106 99.66667 100.33333
## [47,] 118 103.66667 99.66667
## [48,] 109 104.83333 103.66667
## [49,] NA NA 104.83333#membuat plot
ts.plot (Data1.ts,xlab="periode waktu",ylab="Bagi hasil", col="blue",lty=3)
points(Data1.ts)
lines (Data1.sma6,col="red",lwd=2)
lines (ramal.sma6,col="black",lwd= 2)
title("Rataan bergerak Sederhana n=6",cex.main=1,font.main=4 ,col.main="black")
legend("topleft",legend=c("Data aktual","Pemulusan SMA","Ramalan SMA"),
lty=1:3,col=c ("blue","red","black"), cex=0.7)
Mencari nilai keakuratan
# Ukuran Keakuratan
error.sma6<-Data1.ts-ramal.sma6[1:length(Data1.ts)]
## SSE (Sum Square Error)
SSE.sma6 <- sum(error.sma6[7:length(Data1.ts)]^2)
## MSE (Mean Squared Error)
MSE.sma6 <- mean(error.sma6[7:length(Data1.ts)]^2)
## RMSE (Root Mean Square Error)
RMSE.sma6 <- sqrt(mean(error.sma6[7:length(Data1.ts)]^2))
## MAD (Mean Absolute Deviation)
MAD.sma6 <- mean(abs(error.sma6[7:length(Data1.ts)]))
## MAPE (Mean Absolute Percentage Error)
r.error.sma6 <- (error.sma6[7:length(Data1.ts)]/ramal.sma6[7:length(Data1.ts)])*100 # Relative Error
MAPE.sma6 <- mean(abs(r.error.sma6))
performa <- data.frame(
"Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Simple Moving Average N=6" = c(SSE.sma6, MSE.sma6, MAPE.sma6, RMSE.sma6, MAD.sma6))
performaDMA m = 3
bagihasil.dma3<-SMA(Data1.sma3,n=3)
At<-2*Data1.sma3-bagihasil.dma3
Bt<-2/(3-1)*(Data1.sma3-bagihasil.dma3)
pemulusan.dma3<-At+Bt
ramal.dma3<-c(NA,pemulusan.dma3)
#forecasting 5 periode ke depan
t = 1:5
f = c()
for (i in t) {
f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}
Datagab.dma3<-cbind(bagihasil_aktual=c(Data1.ts,rep(NA,5)),
pemulusan.sma = c(Data1.sma3, rep(NA,5)),
pemulusan.dma =c(pemulusan.dma3,rep(NA,5)),
At=c(At, rep(NA,5)),
Bt=c(Bt,rep(NA,5)),
ramalan = c(ramal.dma3, f[-1]))
Datagab.dma3## bagihasil_aktual pemulusan.sma pemulusan.dma At Bt
## [1,] 86 NA NA NA NA
## [2,] 90 NA NA NA NA
## [3,] 86 87.33333 NA NA NA
## [4,] 102 92.66667 NA NA NA
## [5,] 111 99.66667 112.55556 106.11111 6.4444444
## [6,] 92 101.66667 109.00000 105.33333 3.6666667
## [7,] 96 99.66667 98.33333 99.00000 -0.6666667
## [8,] 106 98.00000 94.44444 96.22222 -1.7777778
## [9,] 102 101.33333 104.66667 103.00000 1.6666667
## [10,] 104 104.00000 109.77778 106.88889 2.8888889
## [11,] 99 101.66667 100.33333 101.00000 -0.6666667
## [12,] 92 98.33333 92.33333 95.33333 -3.0000000
## [13,] 87 92.66667 82.88889 87.77778 -4.8888889
## [14,] 85 88.00000 78.00000 83.00000 -5.0000000
## [15,] 94 88.66667 86.44444 87.55556 -1.1111111
## [16,] 112 97.00000 108.55556 102.77778 5.7777778
## [17,] 128 111.33333 136.00000 123.66667 12.3333333
## [18,] 101 113.66667 126.33333 120.00000 6.3333333
## [19,] 124 117.66667 124.55556 121.11111 3.4444444
## [20,] 124 116.33333 117.22222 116.77778 0.4444444
## [21,] 118 122.00000 128.66667 125.33333 3.3333333
## [22,] 120 120.66667 122.66667 121.66667 1.0000000
## [23,] 124 120.66667 119.77778 120.22222 -0.4444444
## [24,] 130 124.66667 130.00000 127.33333 2.6666667
## [25,] 108 120.66667 118.00000 119.33333 -1.3333333
## [26,] 109 115.66667 106.33333 111.00000 -4.6666667
## [27,] 128 115.00000 110.77778 112.88889 -2.1111111
## [28,] 124 120.33333 127.00000 123.66667 3.3333333
## [29,] 106 119.33333 121.55556 120.44444 1.1111111
## [30,] 122 117.33333 114.00000 115.66667 -1.6666667
## [31,] 119 115.66667 112.11111 113.88889 -1.7777778
## [32,] 125 122.00000 129.33333 125.66667 3.6666667
## [33,] 127 123.66667 130.11111 126.88889 3.2222222
## [34,] 105 119.00000 113.88889 116.44444 -2.5555556
## [35,] 103 111.66667 98.77778 105.22222 -6.4444444
## [36,] 105 104.33333 89.66667 97.00000 -7.3333333
## [37,] 93 100.33333 90.11111 95.22222 -5.1111111
## [38,] 107 101.66667 100.77778 101.22222 -0.4444444
## [39,] 103 101.00000 101.00000 101.00000 0.0000000
## [40,] 110 106.66667 113.77778 110.22222 3.5555556
## [41,] 94 102.33333 100.33333 101.33333 -1.0000000
## [42,] 102 102.00000 98.66667 100.33333 -1.6666667
## [43,] 94 96.66667 89.33333 93.00000 -3.6666667
## [44,] 96 97.33333 94.66667 96.00000 -1.3333333
## [45,] 106 98.66667 100.88889 99.77778 1.1111111
## [46,] 106 102.66667 108.88889 105.77778 3.1111111
## [47,] 118 110.00000 122.44444 116.22222 6.2222222
## [48,] 109 111.00000 117.22222 114.11111 3.1111111
## [49,] NA NA NA NA NA
## [50,] NA NA NA NA NA
## [51,] NA NA NA NA NA
## [52,] NA NA NA NA NA
## [53,] NA NA NA NA NA
## ramalan
## [1,] NA
## [2,] NA
## [3,] NA
## [4,] NA
## [5,] NA
## [6,] 112.55556
## [7,] 109.00000
## [8,] 98.33333
## [9,] 94.44444
## [10,] 104.66667
## [11,] 109.77778
## [12,] 100.33333
## [13,] 92.33333
## [14,] 82.88889
## [15,] 78.00000
## [16,] 86.44444
## [17,] 108.55556
## [18,] 136.00000
## [19,] 126.33333
## [20,] 124.55556
## [21,] 117.22222
## [22,] 128.66667
## [23,] 122.66667
## [24,] 119.77778
## [25,] 130.00000
## [26,] 118.00000
## [27,] 106.33333
## [28,] 110.77778
## [29,] 127.00000
## [30,] 121.55556
## [31,] 114.00000
## [32,] 112.11111
## [33,] 129.33333
## [34,] 130.11111
## [35,] 113.88889
## [36,] 98.77778
## [37,] 89.66667
## [38,] 90.11111
## [39,] 100.77778
## [40,] 101.00000
## [41,] 113.77778
## [42,] 100.33333
## [43,] 98.66667
## [44,] 89.33333
## [45,] 94.66667
## [46,] 100.88889
## [47,] 108.88889
## [48,] 122.44444
## [49,] 117.22222
## [50,] 120.33333
## [51,] 123.44444
## [52,] 126.55556
## [53,] 129.66667#membuat plot
ts.plot (Data1.ts,xlab="periode waktu",ylab="Bagi hasil", col="blue",lty=3,ylim=c(80,140))
points(Datagab.dma3[,1])
lines (Datagab.dma3[,3],col="red",lwd=2)
lines (Datagab.dma3[,6],col="black",lwd= 2)
title("Rataan bergerak Bergkita n=3",cex.main=1,font.main=4 ,col.main="black")
legend(0,140,legend=c ("Data aktual","Pemulusan","Ramalan"),lty=1:3,col=c ("blue","red","black"),cex=0.7)
# Ukuran Keakuratan
error.dma3<-Data1.ts-ramal.dma3[1:length(Data1.ts)]
## SSE (Sum Square Error)
SSE.dma3 <- sum(error.dma3[6:length(Data1.ts)]^2)
## MSE (Mean Squared Error)
MSE.dma3 <- mean(error.dma3[6:length(Data1.ts)]^2)
## RMSE (Root Mean Square Error)
RMSE.dma3 <- sqrt(mean(error.dma3[6:length(Data1.ts)]^2))
## MAD (Mean Absolute Deviation)
MAD.dma3 <- mean(abs(error.dma3[6:length(Data1.ts)]))
## MAPE (Mean Absolute Percentage Error)
r.error.dma3 <- (error.dma3[6:length(Data1.ts)]/ramal.dma3[6:length(Data1.ts)])*100 # Relative Error
MAPE.dma3 <- mean(abs(r.error.dma3))
performa <- data.frame(
"Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Double Moving Average N=3" = c(SSE.dma3, MSE.dma3, MAPE.dma3, RMSE.dma3, MAD.dma3))
performaPerbandingan SMA dan DMA
ts.plot (Data1.ts,xlab="periode waktu",ylab="Bagi hasil", col="blue",lty=3,ylim=c(80,140))
points(Data1.ts)
lines (pemulusan.dma3,col="red",lwd=2)
lines (Data1.sma3,col="black",lwd= 2)
lines (Data1.sma6,col="green",lwd= 2)
title("Perbandingan SMA dan DMA",cex.main=1,font.main=4 ,col.main="black")
legend(0,140,legend=c ("Data aktual","Pemulusan SMA 3","Pemulusan SMA 6", "Pemulusan DMA 3"),lty=1:3,col=c("blue","black","green","red"),cex=0.7)
performa <- data.frame(
"Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
"Simple Moving Average N=3" = c(SSE.sma3, MSE.sma3, MAPE.sma3, RMSE.sma3, MAD.sma3),
"Simple Moving Average N=6" = c(SSE.sma6, MSE.sma6, MAPE.sma6, RMSE.sma6, MAD.sma6),
"Double Moving Average N=3" = c(SSE.dma3, MSE.dma3, MAPE.dma3, RMSE.dma3, MAD.dma3))
performaIlustrasi 2
library(forecast)
library(TTR)
library(graphics)
obligasi<-read.csv("Exponential Smoothing.csv", header=TRUE, sep=";")
obligasi.ts<-ts(obligasi$yt)
#SES dengan alpha=0.2
obligasi.ses<- HoltWinters(obligasi.ts,alpha = 0.2, beta=F, gamma=F)
obligasi.ses$SSE #nilai keakuratan## [1] 115.5729obligasi.ses$fitted## Time Series:
## Start = 2
## End = 36
## Frequency = 1
## xhat level
## 2 9.640000 9.640000
## 3 9.231200 9.231200
## 4 9.321760 9.321760
## 5 8.802008 8.802008
## 6 8.543206 8.543206
## 7 8.461565 8.461565
## 8 8.423052 8.423052
## 9 8.636242 8.636242
## 10 8.906193 8.906193
## 11 8.215955 8.215955
## 12 7.741364 7.741364
## 13 7.780891 7.780891
## 14 7.990313 7.990313
## 15 8.107050 8.107050
## 16 7.943840 7.943840
## 17 7.842672 7.842672
## 18 7.653338 7.653338
## 19 6.879470 6.879470
## 20 6.949976 6.949976
## 21 6.616781 6.616781
## 22 6.670225 6.670225
## 23 6.060380 6.060380
## 24 5.857904 5.857904
## 25 6.140323 6.140323
## 26 6.016458 6.016458
## 27 5.992367 5.992367
## 28 5.371093 5.371093
## 29 5.075875 5.075875
## 30 4.451100 4.451100
## 31 3.937480 3.937480
## 32 3.395984 3.395984
## 33 3.621787 3.621787
## 34 3.798830 3.798830
## 35 3.518464 3.518464
## 36 3.305971 3.305971datases<-data.frame(obligasi.ts,c(NA,obligasi.ses$fitted[,1]))
colnames(datases) = c("y","yhat")
datasesobligasi.ramal<-forecast(obligasi.ses,h=5)
obligasi.ramal## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 37 2.760377 0.7729830 4.747771 -0.2790797 5.799833
## 38 2.760377 0.7336249 4.787129 -0.3392728 5.860026
## 39 2.760377 0.6950166 4.825737 -0.3983190 5.919073
## 40 2.760377 0.6571169 4.863637 -0.4562816 5.977035
## 41 2.760377 0.6198882 4.900865 -0.5132180 6.033972SES dengan R
#membuat plot
ts.plot (obligasi.ts,xlab="periode waktu",ylab="obligasi",col="blue",lty=3)
points(obligasi.ts)
lines (datases[,2],col="red",lwd=2)
title("Single Exponential Smoothing alpha=0.2",cex.main=1,font.main=4 ,col.main="black")
legend(0,2,legend=c ("Data aktual","Fitted SES"),lty=1:3,col=c("blue","red"),cex=0.7)
plot(obligasi.ramal,xlab="periode waktu",ylab="obligasi")
DES dengan R
#DES alpha=0.2 dan beta=0.3
obligasi.des<- HoltWinters(obligasi.ts,alpha = 0.2, beta=0.3, gamma=F)
obligasi.des$SSE #nilai keakuratan## [1] 267.0453obligasi.des$fitted## Time Series:
## Start = 3
## End = 36
## Frequency = 1
## xhat level trend
## 3 5.552000 7.596000 -2.044000000
## 4 4.582320 6.378400 -1.796080000
## 5 3.342817 5.010456 -1.667639200
## 6 2.758125 4.175853 -1.417728208
## 7 2.738384 3.833500 -1.095115722
## 8 3.081229 3.844508 -0.763278790
## 9 3.983971 4.362783 -0.378812517
## 10 5.165686 5.184376 -0.018690747
## 11 5.222217 5.223549 -0.001331887
## 12 5.382288 5.346373 0.035915115
## 13 6.082949 5.893631 0.189317809
## 14 6.985980 6.631959 0.354020896
## 15 7.752886 7.303584 0.449302112
## 16 8.082098 7.660509 0.421588958
## 17 8.336221 7.953278 0.382943097
## 18 8.344707 8.048177 0.296529823
## 19 7.455453 7.432565 0.022887414
## 20 7.420243 7.410762 0.009480242
## 21 6.874300 6.992994 -0.118694310
## 22 6.758127 6.876240 -0.118112293
## 23 5.824362 6.130702 -0.306339942
## 24 5.316168 5.669090 -0.352921665
## 25 5.471243 5.706934 -0.235691743
## 26 5.248488 5.481194 -0.232706301
## 27 5.184135 5.377990 -0.193855569
## 28 4.392764 4.724508 -0.331743650
## 29 3.931602 4.293211 -0.361609495
## 30 3.055296 3.535681 -0.480385602
## 31 2.270113 2.820837 -0.550723351
## 32 1.448960 2.062091 -0.613130149
## 33 1.635601 2.064168 -0.428567779
## 34 1.953597 2.209880 -0.256283816
## 35 1.812598 2.042277 -0.229679617
## 36 1.750203 1.941278 -0.191075480datades<-data.frame(obligasi.ts,c(NA,NA,obligasi.des$fitted[,1]))
colnames(datades) = c("y","yhat")
obligasides.ramal<-forecast(obligasi.des,h=5)
obligasides.ramal## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 37 1.2543545 -2.209531 4.718240 -4.043201 6.551910
## 38 0.9929469 -2.586103 4.571997 -4.480738 6.466632
## 39 0.7315392 -3.015226 4.478304 -4.998643 6.461722
## 40 0.4701316 -3.501119 4.441382 -5.603372 6.543635
## 41 0.2087239 -4.044950 4.462398 -6.296709 6.714157#membuat plot
ts.plot (obligasi.ts,xlab="periode waktu",ylab="obligasi",col="blue",lty=3)
points(obligasi.ts)
lines (datades[,2],col="red",lwd=2)
title("Double Exponential Smoothing alpha=0.2 beta 0.3",cex.main=1,font.main=4 ,col.main="black")
legend(0,2,legend=c ("Data aktual","Fitted DES"),lty=1:3,col=c("blue","red"),cex=0.7)
plot(obligasides.ramal,xlab="periode waktu",ylab="obligasi")
Perbandingan SES dengan DES
ts.plot (obligasi.ts,xlab="periode waktu",ylab="obligasi",col="blue",lty=3)
points(obligasi.ts)
lines (datases[,2],col="red",lwd=2)
lines (datades[,2],col="black",lwd= 2)
title("Perbandingan SES dan DES",cex.main=1,font.main=4 ,col.main="black")
legend(0,2.3,legend=c ("Data aktual","Fitted SES","Fitted DES"),lty=1:3,col=c("blue","red","black"),cex=0.7)
Ilustrasi 3
Winter’s Aditif dengan R
Laba1<-read.csv("Winter.csv", header=TRUE, sep=";")
str(Laba1)## 'data.frame': 48 obs. of 2 variables:
## $ t : int 1 2 3 4 5 6 7 8 9 10 ...
## $ yt: int 11354 11273 17050 18007 20288 21619 18749 20429 19689 17904 ...#membentuk objek time series
Laba1.ts<-ts(Laba1$yt,frequency=12)
#melakukan Winters aditif method
Laba1aditif.hw<-HoltWinters(Laba1.ts)## Warning in HoltWinters(Laba1.ts): optimization difficulties: ERROR:
## ABNORMAL_TERMINATION_IN_LNSRCHLaba1aditif.hw$SSE #nilai keakuratan## [1] 13160801Laba1aditif.hw$fitted## xhat level trend season
## Jan 2 12421.06 17515.14 68.54283 -5162.62500
## Feb 2 13119.20 17602.26 70.31546 -4553.37500
## Mar 2 18524.14 17690.26 72.00398 761.87500
## Apr 2 17635.79 17778.14 73.51910 -215.87500
## May 2 22265.15 17865.18 74.80907 4325.16667
## Jun 2 21491.14 17951.67 75.92434 3463.54167
## Jul 2 19360.77 18036.69 76.79171 1247.29167
## Aug 2 20937.81 18102.67 75.76036 2759.37500
## Sep 2 20116.85 18200.98 77.91188 1837.95833
## Oct 2 18296.89 18252.23 75.36782 -30.70833
## Nov 2 17323.31 18340.66 76.61408 -1093.95833
## Dec 2 15140.59 18403.91 75.33930 -3338.66667
## Jan 3 13534.49 18475.01 74.93457 -5015.45307
## Feb 3 14260.91 18594.88 79.22248 -4413.18656
## Mar 3 19651.77 18683.94 80.16184 887.66733
## Apr 3 18730.87 18759.89 79.75964 -108.77555
## May 3 23333.26 18836.08 79.41911 4417.76174
## Jun 3 22504.00 18891.33 77.11242 3535.55499
## Jul 3 20242.40 19000.56 80.17764 1161.66371
## Aug 3 22101.75 19083.32 80.42434 2938.00397
## Sep 3 20900.55 19190.81 83.00683 1626.73919
## Oct 3 19418.78 19263.96 82.06603 72.76210
## Nov 3 18184.33 19305.89 78.23596 -1199.79686
## Dec 3 16141.38 19430.94 82.70414 -3372.26898
## Jan 4 14930.91 19508.18 82.18228 -4659.45085
## Feb 4 15310.71 19566.04 79.86117 -4335.19638
## Mar 4 20580.37 19646.20 79.88998 854.27486
## Apr 4 19681.23 19737.32 80.96119 -137.04801
## May 4 24144.48 19835.61 82.61550 4226.24925
## Jun 4 23765.33 19894.90 80.38872 3790.04404
## Jul 4 21240.90 19978.10 80.65733 1182.14589
## Aug 4 23299.80 20066.04 81.35188 3152.41458
## Sep 4 21783.10 20152.62 81.85148 1548.62927
## Oct 4 20050.60 20215.76 80.06596 -245.22840
## Nov 4 19548.97 20297.56 80.23146 -828.82772
## Dec 4 17020.20 20357.50 78.29449 -3415.59615datawintersaditif<- data.frame(Laba1.ts[13:length(Laba1.ts)],Laba1aditif.hw$fitted[,1])
colnames(datawintersaditif) = c("y","yhat")
datawintersaditifLaba1aditif.ramal<-forecast(Laba1aditif.hw,h=12)
Laba1aditif.ramal## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 5 15657.12 14878.82 16435.42 14466.81 16847.43
## Feb 5 16254.35 15475.56 17033.14 15063.29 17445.41
## Mar 5 21608.24 20828.87 22387.62 20416.29 22800.19
## Apr 5 20743.21 19963.15 21523.26 19550.22 21936.19
## May 5 24862.15 24081.32 25642.99 23667.97 26056.33
## Jun 5 24711.00 23929.28 25492.73 23515.46 25906.55
## Jul 5 22216.34 21433.61 22999.08 21019.25 23413.43
## Aug 5 24248.30 23464.43 25032.17 23049.48 25447.13
## Sep 5 22532.67 21747.53 23317.81 21331.91 23733.43
## Oct 5 20978.67 20192.13 21765.21 19775.76 22181.58
## Nov 5 20298.39 19510.30 21086.48 19093.11 21503.67
## Dec 5 17915.52 17125.73 18705.32 16707.64 19123.41#membuat plot
plot(Laba1aditif.hw,xlab="periode waktu",ylab="Laba1")
points(Laba1.ts)
legend(4.25,14000,legend=c ("Data aktual","Fitted Winter's"),lty=1:3,col=c("blue","red"),cex=0.7)
plot(Laba1aditif.ramal,xlab="periode waktu",ylab="Laba1")
Winter’s Multiplikatif dengan R
#melakukan Winters multiplikatif method
Laba1multiplikatif.hw<-HoltWinters(Laba1.ts, seasonal="mult")
Laba1multiplikatif.hw$SSE #nilai keakuratan## [1] 13910813Laba1multiplikatif.hw$fitted## xhat level trend season
## Jan 2 12577.93 17515.14 68.54283 0.7153183
## Feb 2 13249.39 17603.91 69.08577 0.7496969
## Mar 2 18494.38 17692.28 69.60324 1.0412394
## Apr 2 17629.77 17779.26 70.06928 0.9876995
## May 2 22162.53 17864.23 70.46928 1.2357345
## Jun 2 21416.75 17947.75 70.81945 1.1885929
## Jul 2 19384.79 18028.96 71.09833 1.0709794
## Aug 2 20988.38 18088.44 70.78661 1.1557968
## Sep 2 20121.51 18178.65 71.30766 1.1025513
## Oct 2 18255.91 18223.85 70.60708 0.9978931
## Nov 2 17258.64 18309.94 71.02255 0.9389410
## Dec 2 15037.86 18368.09 70.67711 0.8155568
## Jan 3 13366.93 18437.57 70.64515 0.7222156
## Feb 3 14114.06 18583.11 72.65468 0.7565521
## Mar 3 19682.93 18676.50 73.21116 1.0497725
## Apr 3 18716.08 18744.37 73.06772 0.9946136
## May 3 23481.78 18814.10 72.97827 1.2432723
## Jun 3 22614.06 18862.04 72.30637 1.1943409
## Jul 3 20274.13 18960.01 72.99496 1.0652093
## Aug 3 22282.41 19034.57 73.03700 1.1661538
## Sep 3 20915.30 19127.12 73.56062 1.0893002
## Oct 3 19359.56 19190.49 73.28721 1.0049721
## Nov 3 18010.48 19222.95 72.19177 0.9334205
## Dec 3 15837.12 19355.50 73.81120 0.8151148
## Jan 4 14566.11 19435.12 73.96718 0.7466321
## Feb 4 14939.55 19491.15 73.48581 0.7635995
## Mar 4 20584.92 19582.01 73.95200 1.0472609
## Apr 4 19606.05 19667.32 74.25674 0.9931353
## May 4 24389.09 19762.95 74.83046 1.2294263
## Jun 4 24018.60 19810.47 74.09748 1.2079018
## Jul 4 21269.57 19879.76 73.96842 1.0659450
## Aug 4 23562.98 19960.12 74.13997 1.1761345
## Sep 4 21803.15 20031.23 74.05883 1.0844481
## Oct 4 19898.96 20086.12 73.54446 0.9870676
## Nov 4 19387.24 20166.91 73.73883 0.9578365
## Dec 4 16585.96 20223.79 73.28631 0.8171601datawintersmultiplikatif<- data.frame(Laba1.ts[13:length(Laba1.ts)],Laba1multiplikatif.hw$fitted[,1])
colnames(datawintersmultiplikatif) = c("y","yhat")
datawintersmultiplikatifLaba1multiplikatif.ramal<-forecast(Laba1multiplikatif.hw,h=12)
Laba1multiplikatif.ramal## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 5 15101.44 14841.53 15361.35 14703.95 15498.93
## Feb 5 15737.29 15475.72 15998.87 15337.24 16137.34
## Mar 5 21605.18 21339.13 21871.22 21198.30 22012.06
## Apr 5 20648.63 20381.48 20915.77 20240.07 21057.18
## May 5 25127.98 24854.79 25401.17 24710.17 25545.79
## Jun 5 25015.03 24739.95 25290.11 24594.33 25435.72
## Jul 5 22258.20 21984.55 22531.85 21839.69 22676.71
## Aug 5 24548.75 24269.80 24827.71 24122.13 24975.38
## Sep 5 22561.50 22283.55 22839.44 22136.42 22986.57
## Oct 5 20836.54 20559.34 21113.75 20412.59 21260.49
## Nov 5 20087.76 19809.77 20365.74 19662.62 20512.89
## Dec 5 17409.67 17319.12 17500.22 17271.18 17548.15#membuat plot
plot(Laba1multiplikatif.hw,xlab="periode waktu",ylab="Laba1")
points(Laba1.ts)
legend(4.25,14000,legend=c ("Data aktual","Fitted Winter's"),lty=1:3,col=c("blue","red"),cex=0.7)
plot(Laba1multiplikatif.ramal,xlab="periode waktu",ylab="Laba1")
Perbandingan Winter’s Aditif dengan Winter’s Multiplikatif
ts.plot (Laba1.ts,xlab="periode waktu",ylab="Laba1", col="blue",lty=3)
points(Laba1.ts)
lines (datawintersaditif[,2],col="red",lwd=2)
lines (datawintersmultiplikatif[,2],col="black",lwd= 2)
title("Perbandingan Winter's Aditif dan Multiplikatif",cex.main=1,font.main=4 ,col.main="black")
legend(3.75,14000,legend=c ("Data aktual","Ramalan aditif","Ramalan Multiplikatif"),lty=1:3,col=c("blue","red","black"),cex=0.7)
ARIMA
Sekumpulan data deret waktu dinyatakan stasioner jika nilai rata-rata (mean) dan ragam (variance) dari data tersebut tidak mengalami perubahan secara sistematik sepanjang waktu. Pemodelan deret waktu mengharuskan data yang stasioner. Karena kestasioneran data berkaitan dengan ketepatan estimasi.
Pemeriksaan Kestasioneran Data
Pemeriksaan kestasioneran data dapat dilakukan dengan uji formal yaitu menggunakan uji Augmented Dickey Fuller (ADF).
Hipotesis pada uji ADF:
\(H_0\): Data mempunyai unit root (tidak stasioner)
\(H_1\): Data tidak mempunyai unit root (stasioner)
Selain dengan uji formal ADF, ketidakstasioneran juga dapat dilihat dari plot Fungsi autokorelasi pada lag(k), ACF(k).
Differencing
Jika data tidak stasioner dalam rata-rata (mean), maka perlu dilakukan differencing.
Fungsi Autokorelasi (ACF)
Fungsi autokorelasi pada lag k, ACF(k) merupakan korelasi data deret waktu dengan dirinya sendiri k unit waktu sebelumnya. ACF merupakan salah satu alat utama untuk memeriksa trend, pola musiman, dan kandidat model. Bentuk kurva ACF digunakan untuk memperkirakan model deret waktu yang tepat.
Fungsi Autokorelas Parsial (PACF)
Fungsi autokorelasi parsial pada lag k, PACF(k) merupakan korelasi data deret waktu dengan dirinya sendiri pada lag k unit waktu sebelumnya, setelah menghilangkan efek yang terletak diantara kedua pengamatan tersebut.
Kurva ACF dan PACF biasa disebut dengan korelogram.
Plot ACF dan PACF sangat berguna dalam memprediksikan orde p, q, P dan Q dalam model, dengan melihat dua kemungkinan bentuk dasar: cuts off (terpotong) atau dies down (menurun dengan cepat) dengan pola sinus/eksponensial.
Model Time Series Stasioner
-
Model Autoregressive (AR)
Berdasarkan deviasi acak saat ini dan data masa lalu.
Model AR mempunyai ordo yang besarnya dinotasikan dengan huruf “p”, sehingga dinotasikan dengan AR(p)
Model AR mengasumsikan tiap observasi dibentuk oleh rata-rata tertimbang pengamatan masa lalu, p periode ke belakang dan deviasi periode sekarang.
Model umum AR(p) ditulis:
-
Model Moving Average (MA)
Berdasarkan data masa lalu.
Modek MA mempunyai ordo yang besarnya dinotasikan dengan hurut “q”, sehingga dinotasikan dengan MA(q). Model MA mengasumsikan tiapobservasi dibentuk oleh rata-rata q periode ke belakang.
Model umum MA(q) ditulis:
-
Model campuran antaran AR dan MA (ARMA)
Adakalanya proses acak yang stasioner mempunyai dua karakteristik AR dan MA. Proses acak seperti ini perlu didekati dengan model campuran antara autoregressive dan moving average, disebut ARMA(p,q).
Modul umum ARMA(p,q):
Acuan ACF dan PACF dalam memperkirakan model deret waktu
Pendugaan Parameter dan Diagnostik Model
Metode pendugaan parameter yang umum digunakan adalam maximum likelihood estimation (MLE). Pemilihan model terbaik dari beberapa model yang telah diidentifikasi sebelumnya dilakukan dengan menggunakan Akaike’s Information Criterion (AIC) dan Bayesian Information Criterion (BIC).
Pemeriksaan Kenormalan Sisaan
Dilihat dari Normal Quantile Plot, jika plotnya menyerupai garis lurus, maka asumsi kenormalan terpenuhi.
Dapat dilihat juga dari histogram data apakah memiliki sebaran simetris.
Pemeriksaan Kebebasan Sisaan
Dilihat dari plot pencaran residu(sisaan), jika plotnya tidak membentuk pola tertentu, maka asumsi kebebasan terpenuhi.
Dilihat dari ACF dan PACF sisaan, jika ACF dan PACF tidak ada yang signifikan, maka sisaan saling bebas.
Dilihat dari uji Ljung-Box (Uji Q), jika nilai p-value cukup besar, berarti sisaan sudah saling bebas (tidak ada autokorelasi).
Overfitting Model
Jika model tentatif yang memenuhi ketentuan-ketentuan tadi telah berhasil diperoleh, maka usaha berikutnya adalah melakukan overfitting, yaitu menambah orde model. Kadang-kadang usaha ini dapat dmemberikan alternatif model yang lebih baik, dalam artian penyimpangannya lebih kecil, namun tetap memenuhi syarat signifikansi dan stasioneritas parameter secara kebebasan dan kenormalan sisaan.
Setelah semua tahapan ini dikerjakan barulah dapat ditentukan model finalnya, yaitu model terbaik yang (idealnya)“lolos” dari seluruh batasan asumsi. Dengan model akhir ini dapat dilakukan peramalan nilai-nilai variabel dimaksud untuk beberapa periode yang akan datang.
Ilustrasi 1
dt.arima <- read.csv("Bagi Hasil.csv",header=T,sep=";")
str(dt.arima)## 'data.frame': 100 obs. of 2 variables:
## $ t : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Xt: num 121 120 122 122 123 ...Bagi.hasil <- ts(dt.arima$Xt)
Bagi.hasil## Time Series:
## Start = 1
## End = 100
## Frequency = 1
## [1] 120.8399 120.4160 122.3064 121.5694 122.8942 123.3013 124.1932 124.8930
## [9] 125.8147 125.8902 126.6755 124.1835 125.5775 124.7087 124.9495 125.8825
## [17] 123.9903 124.6912 125.2508 125.9463 123.3904 123.3333 122.6446 123.6364
## [25] 123.5462 123.9536 123.9251 122.7099 122.1322 122.5747 123.7589 122.6723
## [33] 124.0033 123.2013 123.8163 125.2421 126.3553 125.5886 125.1816 123.7251
## [41] 123.4024 123.9143 124.6083 125.0931 124.3820 124.6854 124.7064 124.7050
## [49] 126.2886 127.5864 128.0930 127.2421 124.8627 123.8595 124.4472 123.3074
## [57] 123.8689 124.1007 126.1777 126.9949 125.9848 125.1947 123.4171 124.2836
## [65] 123.4681 125.9491 126.4359 127.6295 128.2976 127.4389 126.4682 125.3810
## [73] 125.5677 124.6353 124.2816 125.8262 126.1171 125.0234 125.2856 124.6220
## [81] 126.8500 127.5602 127.6669 126.7697 126.4150 123.8775 123.6915 125.7046
## [89] 123.9605 124.9237 124.7858 126.4416 126.8498 125.6157 124.8311 124.0431
## [97] 123.6912 124.9995 124.3380 123.7311ts.plot(Bagi.hasil, type="l", ylab="Bagi Hasil")
Partisi Data Data Bagi Hasil akan dibagi menjadi dua bagian yaitu data training (80%) dan data testing (20%). Data training digunakan untuk membangun model, sedangkan data testing digunakan untuk validasi model.
dataTrain <- subset(Bagi.hasil, start = 1, end = 80)
head(dataTrain)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] 120.8399 120.4160 122.3064 121.5694 122.8942 123.3013dataTest <- subset(Bagi.hasil, start = 81, end = 100)
head(dataTest)## Time Series:
## Start = 81
## End = 86
## Frequency = 1
## [1] 126.8500 127.5602 127.6669 126.7697 126.4150 123.8775Step 1 Periksa Kestastioneran
#install.packages("fpp")
library(fpp)## Warning: package 'fpp' was built under R version 4.2.3## Loading required package: fma## Warning: package 'fma' was built under R version 4.2.3## Loading required package: expsmooth## Warning: package 'expsmooth' was built under R version 4.2.3## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: tseries## Warning: package 'tseries' was built under R version 4.2.3#Cek Stasioner data
#Plot TS
plot(dataTrain, main = "Plot Bagi Hasil", ylab = "Bagi hasil")
#Melalui Plot ACF
acf(dataTrain, main = "Plot ACF Bagi Hasil")
pacf(dataTrain, main="Plot PACF Bagi Hasil")
#Melaui Uji Augmented Dickey-Fuller (ADF)
adf.test(Bagi.hasil, alternative = "stationary")##
## Augmented Dickey-Fuller Test
##
## data: Bagi.hasil
## Dickey-Fuller = -3.8475, Lag order = 4, p-value = 0.0194
## alternative hypothesis: stationaryHipotesis:
\(H_0\): Data tidak stasioner
\(H_1\): Data stasioner
Jika \(p-value < \alpha\) maka data stastioner.
Dapat dilihat bahwa p-value = 0.0194 \(< \alpha = 0.05\) maka data Bagi Hasil
merupakan data yang stasioner.
Step 2 Identifikasi Model
#Identifikasi model melalui Plot ACF dan PACF
par(mfrow = c(1,2))
acf(dataTrain, main = "Plot ACF bagi hasil", lag.max=20)
pacf(dataTrain, main = "Plot PACF bagi hasil", lag.max=20)
Berdasarkan plot, dapat dilihat bahwa ACF menurun eksponensial. Sedangkan PACF signifikan pada lag ke-1 dan cut off pada lag ke-2.
Step 3 Pendugaan Parameter
#Fit Model
fit.ar1 <- Arima(dataTrain, order = c(1,0,0))
summary(fit.ar1)## Series: dataTrain
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.7849 124.5010
## s.e. 0.0725 0.5041
##
## sigma^2 = 1.036: log likelihood = -114.4
## AIC=234.8 AICc=235.12 BIC=241.95
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05333909 1.005112 0.8281074 0.03662076 0.6644191 0.9399806
## ACF1
## Training set 0.02153699fit.ar1$coef## ar1 intercept
## 0.7848972 124.5009592#t hitung untuk ar1
thit.ar1 <- fit.ar1$coef[1]/0.0725
thit.ar1## ar1
## 10.82617p.value.ar1 = 2*pt(-abs(thit.ar1), df = length(dataTrain)-1)
p.value.ar1## ar1
## 2.845985e-17#t hitung untuk intersep
thit.int <- fit.ar1$coef[2]/0.5041
thit.int## intercept
## 246.9767p.value.int = 2*pt(-abs(thit.int), df = length(dataTrain)-1)
p.value.int## intercept
## 7.349694e-116Model: \(y_t = 124.50 + 0.78y_{t-1}+e_t\)
auto.arima(dataTrain, trace = TRUE)##
## ARIMA(2,1,2) with drift : 241.9522
## ARIMA(0,1,0) with drift : 237.1099
## ARIMA(1,1,0) with drift : 238.9546
## ARIMA(0,1,1) with drift : 239.0098
## ARIMA(0,1,0) : 235.1659
## ARIMA(1,1,1) with drift : 239.9066
##
## Best model: ARIMA(0,1,0)## Series: dataTrain
## ARIMA(0,1,0)
##
## sigma^2 = 1.12: log likelihood = -116.56
## AIC=235.11 AICc=235.17 BIC=237.48lmtest::coeftest(fit.ar1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.784897 0.072506 10.825 < 2.2e-16 ***
## intercept 124.500959 0.504143 246.956 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Step 4 Diagnostik Model
Diagnostik model dilakukan untuk mengetahui apakah sisaan dari model yang dibangun telah memenuhi asumsi pemodelan. Asumsi sisaan yang harus terpenuhi pada model ARIMA(p,d,q) adalah sebagai berikut:
Sisaan saling bebas (tidak ada autokorelasi) dan identik
Sisaan mengikuti sebaran Normal \((0,\sigma^2)\)
Eksploratif
#Plot Histogram sisaan
sisaan.model <- fit.ar1$residuals
#Eksplorasi
h <- hist(sisaan.model, col="lightblue", xlab="Sisaan", main="Histogram Sisaan")
xfit <- seq(min(sisaan.model), max(sisaan.model), length=100)
yfit <- dnorm(xfit,mean=mean(sisaan.model), sd=sd(sisaan.model))
yfit <- yfit*diff(h$mids[1:2]*length(sisaan.model))
lines(xfit, yfit, col="red", lwd=2)
par(mfrow=c(1,2))
#Plot Q-Q
qqnorm(sisaan.model, ylab="Plot QQ Sisaan Model")
qqline(sisaan.model, col="red", lwd=2)
#Plot sebaran sisaan
plot(sisaan.model, type="p", main="Plot Sisaan")
#Plot ACF dan PACF sisaan Model
par(mfrow=c(1,2))
acf(sisaan.model, main = "Plot ACF sisaan", lag.max=20)
pacf(sisaan.model, main = "Plot PACF sisaan", lag.max=20)
checkresiduals(fit.ar1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 10.527, df = 9, p-value = 0.3095
##
## Model df: 1. Total lags used: 10Berdasarkan plot di atas dapat dilihat pada Normal Q-Q Plot menunjukkan bahwa sisaan mengikuti sebaran normal karena banyak titik-titik yang berada di sekitar garis. Jika dilihat dari plot ACF bahwa ada lag yang signifikan pada lag ke-1. Hal ini menunjukkan bahwa ada gejala autokorelasi pada sisaan. Selanjutnya, untuk memastikan kembali akan dilakukan uji asumsi secara formal.
Uji Non-Autokorelasi dengan
Ljung-Box test
Hipotesis:
\(H_0\): Sisaan saling bebas
\(H_1\): Sisaan tidak saling bebas
Keputusan: Tolak \(H_0\) jika \(p-value < \alpha\)
#L-Jung Box test
test <- Box.test(sisaan.model, lag = 24, type = "Ljung")
test##
## Box-Ljung test
##
## data: sisaan.model
## X-squared = 28.224, df = 24, p-value = 0.2507Berdasarkan hasil uji Ljung-Box di atas dapat dilihat bahwa gagal
tolak \(H_0\) karena diperoleh
p-value = 0.2507 yang berarti bahwa tidak terdapat gejala
autokorelasi pada sisaan dari model ARIMA (1,0,0) atau antar \(e_t\) tidak berkorelasi.
Uji Normalitas dengan
Jarque Bera Test
Hipotesis:
\(H_0\): Sisaan berdistribusi Normal
\(H_1\): Sisaan tidak berdistribusi Normal
Keputusan: Tolak \(H_0\) jika \(p-value < \alpha\)
library(tseries)
jarque.bera.test(residuals(fit.ar1))##
## Jarque Bera Test
##
## data: residuals(fit.ar1)
## X-squared = 1.3808, df = 2, p-value = 0.5014Berdasarkan hasil uji normalitas di atas dapat disimpulkan bahwa \(H_0\) gagal ditolak sehingga sisaan dari model ARIMA (1,0,0) mengikuti sebaran Normal.
Step 5 Overfitting
Pada kasus data bagi hasil, model yang diperoleh adalah AR(1). Sehingga overfitting yang dilakukan adalah pada model AR(2).
Pendugaan Parameter: Model AR(2)
#Overfitting
fit.ar2 <- Arima(dataTrain, order = c(2,0,0))
summary(fit.ar2)## Series: dataTrain
## ARIMA(2,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 mean
## 0.8381 -0.0739 124.5299
## s.e. 0.1108 0.1163 0.4637
##
## sigma^2 = 1.044: log likelihood = -114.2
## AIC=236.4 AICc=236.93 BIC=245.93
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04619852 1.002627 0.8200764 0.03085011 0.6581858 0.9308646
## ACF1
## Training set -0.04964409fit.ar2$coef[1:3]## ar1 ar2 intercept
## 0.8381357 -0.0738928 124.5298562#t hitung untuk koef ar1
thit.ar1 <- fit.ar2$coef[1]/0.1108
thit.ar1## ar1
## 7.564402p.value.ar1 = 2*pt(-abs(thit.ar1), df=length(dataTrain)-1)
p.value.ar1## ar1
## 6.119474e-11#t hitung untuk koef ar1
thit.ar2 <- fit.ar2$coef[2]/0.1163
thit.ar2## ar2
## -0.6353637p.value.ar2 = 2*pt(-abs(thit.ar2), df=length(dataTrain)-1)
p.value.ar2## ar2
## 0.5270274#t hitung untuk intersep
thit.int <- fit.ar2$coef[3]/0.4637
thit.int## intercept
## 268.5569p.value.int <- 2*pt(-abs(thit.int), df=length(dataTrain)-1)
p.value.int## intercept
## 9.899179e-119#Fit Model
fit.ar1 <- Arima(dataTrain, order = c(1,0,0))
summary(fit.ar1)## Series: dataTrain
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.7849 124.5010
## s.e. 0.0725 0.5041
##
## sigma^2 = 1.036: log likelihood = -114.4
## AIC=234.8 AICc=235.12 BIC=241.95
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05333909 1.005112 0.8281074 0.03662076 0.6644191 0.9399806
## ACF1
## Training set 0.02153699lmtest::coeftest(fit.ar1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.784897 0.072506 10.825 < 2.2e-16 ***
## intercept 124.500959 0.504143 246.956 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Overfitting
fit.ar2 <- Arima(dataTrain, order = c(2,0,0))
summary(fit.ar2)## Series: dataTrain
## ARIMA(2,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 mean
## 0.8381 -0.0739 124.5299
## s.e. 0.1108 0.1163 0.4637
##
## sigma^2 = 1.044: log likelihood = -114.2
## AIC=236.4 AICc=236.93 BIC=245.93
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04619852 1.002627 0.8200764 0.03085011 0.6581858 0.9308646
## ACF1
## Training set -0.04964409lmtest::coeftest(fit.ar2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.838136 0.110754 7.5676 3.803e-14 ***
## ar2 -0.073893 0.116286 -0.6354 0.5251
## intercept 124.529856 0.463738 268.5347 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Berdasarkan hasil pendugaan parameter dan overfitting, diperoleh
model AR(1) sebagai model terbaik karena AR(1) memiliki nilai AIC
terkecil dan jika dilihat dari hasil uji signifikansi pendugaan
parameter dengan fungsi coeftest terlihat bahwa seluruh
penduga paramater pada model AR(1) telah signifikan terhadap model.
Sehingga model yang didapatkan adalah: Model AR(1): \(y_t=124.50+0.78y_{t-1}+e_t\)
Validasi Model
#plot dugaan dengan data asli
dugaan.ar1 <- fitted(fit.ar1)
dugaan.ar2 <- fitted(fit.ar2)
cbind(dataTrain,dugaan.ar1,dugaan.ar2)## Time Series:
## Start = 1
## End = 80
## Frequency = 1
## dataTrain dugaan.ar1 dugaan.ar2
## 1 120.8399 123.1084 123.1406
## 2 120.4160 121.6274 121.6466
## 3 122.3064 121.2947 121.3545
## 4 121.5694 122.7785 122.9703
## 5 122.8942 122.2000 122.2129
## 6 123.3013 123.2398 123.3777
## 7 124.1932 123.5594 123.6210
## 8 124.8930 124.2594 124.3385
## 9 125.8147 124.8087 124.8591
## 10 125.8902 125.5321 125.5799
## 11 126.6755 125.5914 125.5751
## 12 124.1835 126.2078 126.2277
## 13 125.5775 124.2518 124.0810
## 14 124.7087 125.3459 125.4335
## 15 124.9495 124.6640 124.6023
## 16 125.8825 124.8530 124.8684
## 17 123.9903 125.5853 125.6325
## 18 124.6912 124.1001 123.9777
## 19 125.2508 124.6503 124.7050
## 20 125.9463 125.0895 125.1222
## 21 123.3904 125.6354 125.6638
## 22 123.3333 123.6293 123.4702
## 23 122.6446 123.5845 123.6112
## 24 123.6364 123.0439 123.0382
## 25 123.5462 123.8224 123.9203
## 26 123.9536 123.7516 123.7714
## 27 123.9251 124.0713 124.1196
## 28 122.7099 124.0490 124.0656
## 29 122.1322 123.0952 123.0492
## 30 122.5747 122.6417 122.6548
## 31 123.7589 122.9890 123.0683
## 32 122.6723 123.9185 124.0282
## 33 124.0033 123.0656 123.0299
## 34 123.2013 124.1103 124.2258
## 35 123.8163 123.4809 123.4553
## 36 125.2421 123.9636 124.0300
## 37 126.3553 125.0827 125.1795
## 38 125.5886 125.9564 126.0072
## 39 125.1816 125.3546 125.2823
## 40 123.7251 125.0352 124.9979
## 41 123.4024 123.8920 123.8072
## 42 123.9143 123.6387 123.6444
## 43 124.6083 124.0405 124.0972
## 44 125.0931 124.5852 124.6411
## 45 124.3820 124.9657 124.9961
## 46 124.6854 124.4076 124.3643
## 47 124.7064 124.6457 124.6711
## 48 124.7050 124.6622 124.6663
## 49 126.2886 124.6611 124.6636
## 50 127.5864 125.9041 125.9910
## 51 128.0930 126.9227 126.9617
## 52 127.2421 127.3203 127.2904
## 53 124.8627 126.6525 126.5398
## 54 123.8595 124.7849 124.6084
## 55 124.4472 123.9975 123.9434
## 56 123.3074 124.4588 124.5101
## 57 123.8689 123.5641 123.5114
## 58 124.1007 124.0049 124.0662
## 59 126.1777 124.1868 124.2190
## 60 126.9949 125.8170 125.9427
## 61 125.9848 126.4584 126.4741
## 62 125.1947 125.6656 125.5671
## 63 123.4171 125.0455 124.9796
## 64 124.2836 123.6502 123.5481
## 65 123.4681 124.3304 124.4057
## 66 125.9491 123.6903 123.6582
## 67 126.4359 125.6376 125.7978
## 68 127.6295 126.0197 126.0225
## 69 128.2976 126.9565 126.9869
## 70 127.4389 127.4809 127.4587
## 71 126.4682 126.8069 126.6896
## 72 125.3810 126.0450 125.9395
## 73 125.5677 125.1917 125.1000
## 74 124.6353 125.3382 125.3368
## 75 124.2816 124.6064 124.5415
## 76 125.8262 124.3288 124.3140
## 77 126.1171 125.5411 125.6347
## 78 125.0234 125.7695 125.7644
## 79 125.2856 124.9110 124.8262
## 80 124.6220 125.1168 125.1268ts.plot(dataTrain,xlab="periode waktu",ylab="Bagi hasil", col="blue",lty=3)
points(dataTrain)
lines (dugaan.ar1,col="red",lwd=2)
lines (dugaan.ar2,col="black",lwd= 2)
legend("topleft",legend=c("Data aktual","AR(1)","AR(2)"),
lty=1:3,col=c ("blue","red","black"), cex=0.7)
Ukuran keakuratan ramalan
#Forecast
forecast.ar1 <- forecast(dataTrain, model = fit.ar1, h=20)
forecast.ar1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 81 124.5960 123.2915 125.9005 122.6009 126.5910
## 82 124.5755 122.9172 126.2339 122.0393 127.1118
## 83 124.5595 122.7167 126.4023 121.7411 127.3778
## 84 124.5469 122.5991 126.4947 121.5680 127.5258
## 85 124.5370 122.5273 126.5468 121.4634 127.6106
## 86 124.5293 122.4823 126.5762 121.3987 127.6598
## 87 124.5232 122.4536 126.5927 121.3580 127.6883
## 88 124.5184 122.4350 126.6018 121.3322 127.7046
## 89 124.5146 122.4228 126.6065 121.3155 127.7138
## 90 124.5117 122.4147 126.6087 121.3046 127.7188
## 91 124.5094 122.4092 126.6096 121.2974 127.7214
## 92 124.5076 122.4054 126.6098 121.2926 127.7226
## 93 124.5062 122.4028 126.6095 121.2893 127.7230
## 94 124.5050 122.4009 126.6092 121.2871 127.7230
## 95 124.5042 122.3996 126.6087 121.2855 127.7228
## 96 124.5035 122.3986 126.6083 121.2844 127.7226
## 97 124.5029 122.3979 126.6080 121.2835 127.7223
## 98 124.5025 122.3974 126.6077 121.2830 127.7221
## 99 124.5022 122.3970 126.6074 121.2825 127.7218
## 100 124.5019 122.3967 126.6072 121.2822 127.7216plot(forecast.ar1)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniongabTestForecast<- cbind(dataTest,forecast.ar1)
df.gab<- as.data.frame(gabTestForecast, digits=3)
df.gabaccuracy(forecast.ar1, dataTest)## ME RMSE MAE MPE MAPE MASE
## Training set 0.05333909 1.005112 0.8281074 0.03662076 0.6644191 0.9399806
## Test set 0.81445620 1.541135 1.2389227 0.63887856 0.9816734 1.4062949
## ACF1 Theil's U
## Training set 0.02153699 NA
## Test set 0.57892128 1.309237Step 6 Peramalan
#Forecast
fit.ar1 <- arima(Bagi.hasil, order=c(1,0,0))
forecast.bagihasil <- forecast(Bagi.hasil, model = fit.ar1, h=10)
forecast.bagihasil## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 101 123.9471 122.6274 125.2668 121.9288 125.9654
## 102 124.1120 122.4515 125.7725 121.5725 126.6515
## 103 124.2380 122.4078 126.0681 121.4390 127.0370
## 104 124.3341 122.4119 126.2564 121.3944 127.2739
## 105 124.4076 122.4337 126.3815 121.3888 127.4264
## 106 124.4637 122.4602 126.4671 121.3997 127.5277
## 107 124.5065 122.4860 126.5270 121.4165 127.5966
## 108 124.5392 122.5089 126.5696 121.4341 127.6443
## 109 124.5642 122.5281 126.6003 121.4503 127.6781
## 110 124.5833 122.5439 126.6227 121.4643 127.7023plot(forecast.bagihasil)
Ilustrasi 2
Data Harga Emas Japan (1979-2021)
hargaEmas <- read.csv("1979-2021.csv", sep = ",")
head(hargaEmas)hargaEmasJPN <- hargaEmas[,4]
head(hargaEmasJPN)## [1] 45160.3 50209.1 50274.3 54144.6 61057.1 60161.4str(hargaEmasJPN)## num [1:511] 45160 50209 50274 54145 61057 ...emasJPN.ts <- ts(hargaEmasJPN)
emasJPN.ts## Time Series:
## Start = 1
## End = 511
## Frequency = 1
## [1] 45160.3 50209.1 50274.3 54144.6 61057.1 60161.4 63972.9 69957.5
## [9] 89088.5 90261.1 104134.1 125630.3 159889.7 161056.0 125800.0 124139.1
## [17] 118403.7 145558.9 140586.5 138824.3 141296.4 134300.0 134552.9 119823.1
## [25] 104952.0 102819.8 108590.3 103761.7 107187.2 95631.8 97366.9 97572.7
## [33] 99122.9 99984.8 88997.8 87899.5 88640.0 85939.7 79495.5 84104.9
## [41] 79261.7 79918.3 88186.9 107111.6 105866.8 116886.8 107083.0 105366.4
## [49] 120786.2 99730.8 98978.0 101942.2 104996.1 99204.6 102097.1 102099.2
## [57] 95265.3 89752.6 94229.8 88442.8 87369.3 92597.7 87261.1 85408.4
## [65] 88828.1 89063.3 83799.7 84035.3 84562.2 82186.8 81327.1 77716.7
## [73] 77640.3 75449.3 82616.6 81705.2 78609.6 78952.5 77544.4 79544.6
## [81] 70455.9 68753.8 65726.5 65490.6 67594.7 60984.2 61718.4 57680.7
## [89] 59899.8 56735.8 54972.4 59682.6 65330.3 65478.2 63105.3 61832.1
## [97] 61447.1 61216.6 61445.0 63736.2 64952.4 65632.6 67909.6 64394.8
## [105] 66604.9 64960.2 65203.7 58714.4 59470.5 54719.1 56689.0 56080.3
## [113] 57140.4 58268.2 58820.2 58358.5 53617.9 51768.7 51497.6 51253.0
## [121] 51430.5 49082.2 50605.3 50121.1 51480.3 53698.5 50443.8 52020.6
## [129] 51173.9 53571.9 58335.1 57287.2 59972.3 60672.5 57987.1 58375.7
## [137] 55420.9 53602.3 54393.8 55780.7 56458.8 49292.1 51237.2 53028.2
## [145] 48101.6 48220.8 50208.2 48765.2 49902.4 50793.2 49851.5 47532.4
## [153] 47121.0 46668.5 47629.4 44151.4 44460.3 45619.6 45481.3 44875.1
## [161] 43087.1 43240.7 45570.0 41899.6 41874.1 41795.6 41568.0 41576.2
## [169] 41222.9 38676.6 38817.3 39404.3 40441.1 40417.0 42136.4 38903.5
## [177] 37737.0 40092.8 40475.4 43592.9 41194.5 39738.5 39991.6 38204.0
## [185] 40565.7 38004.2 38468.4 38642.4 39077.3 37215.4 37909.1 38236.9
## [193] 37096.3 36397.9 33868.8 32746.8 32529.1 32860.5 33750.1 37451.2
## [201] 37881.6 39124.0 39373.3 39921.0 43359.4 42136.4 42330.2 40990.6
## [209] 42218.4 41899.7 41121.1 41968.5 42211.1 43164.3 42253.9 42858.8
## [217] 41943.7 43279.4 43104.4 43168.4 40272.8 38282.5 38639.8 39136.3
## [225] 40094.4 37452.1 37879.1 37733.3 38655.0 37570.5 40139.9 41052.8
## [233] 40682.7 41122.0 41704.1 38601.3 39998.8 34061.7 36224.5 32463.8
## [241] 33186.3 34058.5 33095.3 34210.0 32714.1 31591.4 29333.9 27937.5
## [249] 31834.5 31214.1 29771.6 29708.5 30325.8 32258.9 28384.9 29726.0
## [257] 29314.5 30487.7 30326.3 29539.3 29570.6 28880.8 29809.5 31342.2
## [265] 30757.4 31282.6 32295.0 32508.2 31784.4 33749.2 33206.9 32488.4
## [273] 34917.0 34120.4 33919.6 36238.1 37764.7 39725.9 39946.0 39572.9
## [281] 40532.7 38175.4 36484.9 37087.1 39407.2 38823.4 39088.4 41202.2
## [289] 44078.0 41077.3 39706.5 40160.8 43225.2 41546.0 42763.3 43825.0
## [297] 43345.4 42462.4 43627.3 44609.5 42313.5 43250.6 44079.6 42872.9
## [305] 43464.0 43187.7 43623.5 44693.7 45810.9 45189.1 46607.2 44635.9
## [313] 43673.5 45410.9 45725.4 45711.5 44667.3 48435.0 48088.8 48149.2
## [321] 53638.2 54797.7 59307.0 60549.4 66580.7 64398.7 68670.2 73544.8
## [329] 73207.8 70132.2 72380.1 73127.2 70750.4 70732.3 74794.1 75306.0
## [337] 78684.5 78697.7 78136.1 80904.9 80288.2 80333.5 79237.8 77894.9
## [345] 85456.1 91017.5 86909.7 93142.4 98164.6 101138.0 92915.9 91028.2
## [353] 93504.2 98611.1 99222.0 90401.3 93902.9 71869.3 77577.1 78842.8
## [361] 82580.3 93148.4 90522.7 86898.6 93116.3 90165.2 89388.1 88637.0
## [369] 89154.5 94140.8 101285.0 101240.8 97771.4 98484.6 104232.3 110861.3
## [377] 109924.8 110081.6 101305.5 104626.6 109186.8 108480.7 115895.8 113993.1
## [385] 108721.1 115603.2 119264.3 124590.5 124825.3 121584.2 125703.9 138678.3
## [393] 124869.6 134272.9 135542.0 117795.1 132980.0 143263.8 136815.4 131844.0
## [401] 122170.6 127544.3 126678.2 129077.5 138172.8 137399.7 142343.2 143315.7
## [409] 151925.1 146555.0 150267.5 143095.3 140809.6 118407.3 129287.6 136845.9
## [417] 130162.8 129897.6 128263.3 126599.0 127583.2 135389.2 133030.9 131620.3
## [425] 127232.1 133216.1 132168.7 133570.1 133444.0 130518.2 140374.7 144593.4
## [433] 148035.3 145127.6 142351.0 141240.5 147846.8 143289.4 136086.2 137545.0
## [441] 133418.2 137853.1 130911.0 127512.7 134600.0 139389.3 139032.6 137558.1
## [449] 134427.9 135495.7 137534.9 135422.3 133923.0 133668.1 134238.6 133652.0
## [457] 136530.9 140482.8 138713.6 141171.2 140022.7 139579.2 140057.9 144371.2
## [465] 144432.1 144327.1 143260.8 145431.1 146818.9 140621.2 140791.4 143696.9
## [473] 141819.7 138506.1 136673.1 133333.7 134853.8 137113.2 138252.8 140325.5
## [481] 144009.3 146847.8 143381.3 142816.1 140664.3 151805.7 154996.2 162232.0
## [489] 160523.8 163341.2 159901.0 164615.4 171703.5 173646.5 173694.2 182083.6
## [497] 186241.5 190751.5 207748.9 207567.2 199124.5 196728.6 183789.9 194885.3
## [505] 195130.6 185683.3 186861.0 193213.0 207845.0 195692.0 200376.1ts.plot(emasJPN.ts)
acf(emasJPN.ts, lag.max = 100)
pacf(emasJPN.ts, lag.max = 100)
#UJI ADF
adf.test(emasJPN.ts, alternative = "stationary")##
## Augmented Dickey-Fuller Test
##
## data: emasJPN.ts
## Dickey-Fuller = -0.91431, Lag order = 7, p-value = 0.9514
## alternative hypothesis: stationary#Differencing
emasJPN.dif1 <- diff(emasJPN.ts, differences = 1)
ts.plot(emasJPN.dif1)
points(emasJPN.dif1)
adf.test(emasJPN.dif1, alternative = "stationary")## Warning in adf.test(emasJPN.dif1, alternative = "stationary"): p-value smaller
## than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: emasJPN.dif1
## Dickey-Fuller = -8.4508, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationaryauto.arima(emasJPN.dif1, d=1, trace = TRUE)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2) with drift : Inf
## ARIMA(0,1,0) with drift : 10528.65
## ARIMA(1,1,0) with drift : 10391.68
## ARIMA(0,1,1) with drift : 10178.25
## ARIMA(0,1,0) : 10526.63
## ARIMA(1,1,1) with drift : Inf
## ARIMA(0,1,2) with drift : 10179.28
## ARIMA(1,1,2) with drift : Inf
## ARIMA(0,1,1) : 10176.46
## ARIMA(1,1,1) : Inf
## ARIMA(0,1,2) : 10177.49
## ARIMA(1,1,0) : 10389.66
## ARIMA(1,1,2) : Inf
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(0,1,1) : Inf
## ARIMA(0,1,2) : Inf
## ARIMA(0,1,1) with drift : Inf
## ARIMA(0,1,2) with drift : Inf
## ARIMA(1,1,0) : 10407.23
##
## Best model: ARIMA(1,1,0)## Series: emasJPN.dif1
## ARIMA(1,1,0)
##
## Coefficients:
## ar1
## -0.4913
## s.e. 0.0388
##
## sigma^2 = 44104603: log likelihood = -5201.61
## AIC=10407.21 AICc=10407.23 BIC=10415.68Heteroscedasticity Model
Heteroskedastisitas dalam konteks model time series merujuk pada kondisi di mana varians dari residu atau error model time series tidak konstan sepanjang waktu. Dengan kata lain, heteroskedastisitas menunjukkan bahwa fluktuasi atau variasi dari nilai-nilai residual dalam model time series tidak stabil, dan bisa menjadi lebih besar atau lebih kecil pada beberapa periode waktu tertentu.
Penyebab heteroskedastisitas dalam model time series dapat bervariasi, dan hal ini dapat mempengaruhi keakuratan model serta interpretasi hasilnya. Beberapa penyebab umum heteroskedastisitas dalam model time series meliputi:
Perubahan Regime: Terkadang, dalam data time series, ada perubahan dalam kondisi atau regime yang dapat mempengaruhi volatilitas data. Misalnya, dalam data keuangan, kondisi ekonomi yang berubah bisa menghasilkan fluktuasi yang berbeda dalam volatilitas harga.
Musimanitas: Data time series seringkali menunjukkan pola musiman, yang dapat menyebabkan fluktuasi dalam varians residual. Misalnya, penjualan toko pakaian mungkin lebih volatil selama musim liburan daripada selama bulan-bulan lain.
Perubahan Variabilitas Tren: Jika ada perubahan dalam tren data seiring waktu, variabilitas residualnya juga dapat berubah. Sebagai contoh, jika sebuah produk baru sukses diluncurkan dan mengubah tren penjualan perusahaan, maka variabilitas dalam residu model penjualan mungkin meningkat.
Outlier: Kemunculan outlier atau nilai-nilai yang ekstrem dalam data time series dapat mengakibatkan heteroskedastisitas. Outlier dapat mempengaruhi varians residual pada periode di mana mereka terjadi.
Efek Ganda: Kadang-kadang, dua atau lebih faktor dapat berinteraksi secara kompleks dan menghasilkan fluktuasi yang tidak konstan dalam varians residual.
Heteroskedastisitas bisa menjadi masalah dalam analisis time series karena dapat mengganggu asumsi dasar model statistik, seperti model autoregresi (AR), model moving average (MA), atau model autoregresi terintegrasi moving average (ARIMA). Untuk mengatasi heteroskedastisitas, beberapa langkah yang dapat diambil adalah:
Transformasi Data: Kita dapat mencoba mentransformasi data time series, seperti mengambil logaritma atau diferensiasi data, untuk menjadikan varians residu lebih stabil.
Penggunaan Model yang Sesuai: Menggunakan model time series yang lebih canggih yang memungkinkan untuk mengatasi heteroskedastisitas, seperti model GARCH (Generalized Autoregressive Conditional Heteroskedasticity), yang dirancang khusus untuk menangani fluktuasi varians.
Penggunaan Weighted Least Squares: Jika heteroskedastisitas terdeteksi, kita dapat menggunakan metode Weighted Least Squares (WLS) untuk memberikan bobot berbeda pada observasi yang berbeda dalam upaya mengatasi masalah heteroskedastisitas.
Penelitian lebih lanjut: Penting untuk memahami penyebab heteroskedastisitas dalam data kita, sehingga kita dapat mengambil langkah-langkah yang sesuai untuk mengatasinya. Ini mungkin melibatkan penelitian lebih lanjut dan analisis data yang lebih mendalam.
ARCH Model
ARCH (Autoregressive Conditional Heteroskedasticity) dikembangkan oleh Robert F. Engle pada tahun 1982. Pendekatan ini mengasumsikan bahwa variabilitas dalam data time series bergantung pada nilai-nilai residual sebelumnya (autoregresif) dan bahwa variabilitas tersebut dapat dijelaskan melalui model autoregresi untuk varians.
Model ARCH(p) memodelkan varians residual pada waktu t sebagai fungsi dari nilai residual sebelumnya, yaitu variabel autoregresfid tingkat p (AR(p)). Model ini dirumuskan sebagai berikut:
\[ \sigma_t^2= \alpha_0+ \sum_{i=1}^p \alpha_i .\varepsilon_{t-i}^2 \]
di mana \(\sigma_t^2\) adalah varians residual pada waktu ke-t, \(\alpha_0\) adalah konstanta, \(\alpha_i\) adalah koefisien model ARCH, dan \(\varepsilon_{t-i}^2\) adalah nilai residual kuadrat pada waktu ke-(t-i).
Estimasi parameter model ARCH dapat dilakukan menggunakan metode kuadrat terkecil (least square) atau metode lainnya.
Salah satu kelemahan model ARCH adalah bahwa ia mengasumsikan bahwa variabilitas kondisional hanya bergantung pada nilai residual sebelumnya, sementara dalam situasi nyata, faktor-faktor lain seperti informasi eksternal juga dapat mempengaruhi variabilitas. Oleh karena itu, model GARCH yang lebih kompleks sering digunakan karena memungkinkan untuk memasukkan faktor-faktor tambahan dalam model.
GARCH Model
GARCH (Generalized Autoregressive Conditional Heteroskedasticity) adalah pengembangan dari model ARCH yang lebih fleksibel. Model GARCH mengasumsikan bahwa variabilitas kondisional dalam data time series dapat dijelaskan oleh nilai-nilai residual sebelumnya (autoregresif), nilai-nilai varians sebelumnya (moving average), dan komponen heteroskedastisitas konstan.
Model GARCH(p, q) adalah model yang paling umum digunakan dalam praktik. Model ini dirumuskan sebagai berikut:
\[ \sigma_t^2=\omega+\sum_{i=1}^p\alpha_i.\varepsilon_{t-1}^2+\sum_{j=1}^q\beta_j.\sigma_{t-j}^2 \]
di mana \(\sigma_t^2\) adalah varians residual pada waktu ke-t, \(\omega\) adalah komponen heteroskedastisitas konstan, \(\alpha_i\) adalah koefisein model GARCH untuk nilai residual, \(\varepsilon_{t-i}^2\) adalah nilairesidula kuadrat pada waktu ke \(t-i\) , \(\beta_j\) adalah koefisien model GARCH untuk nilai varians sebelumnya, dan \(\sigma_{t-j}^2\) adalah varians residual pada waktu ke \(t-j\).
Model GARCH memiliki keunggulan dibandingkan dengan model ARCH karena lebih fleksibel dan memungkinkan perubahan dalam variabilitas seiring waktu. Dengan kata lain, model GARCH dapat mengatasi lebih baik fluktuasi volatilitas yang kompleks.
Model GARCH tetap merupakan model statistik yang berdasarkan pada asumsi tertentu, sehingga tidak selalu cocok untuk semua jenis data time series. Selain itu, dalam aplikasi praktis, estimasi parameter GARCH bisa menjadi rumit, terutama jika ada banyak parameter yang harus diestimasi.
GARCH-in-Mean
Biasanya, dalam model GARCH (Generalized Autoregressive Conditional Heteroskedasticity), asumsi yang umum adalah bahwa nilai rata-rata dari residu adalah nol (zero mean). Ini berarti, secara statistik, nilai-nilai residu seharusnya memiliki rata-rata yang mendekati nol.
Asumsi rata-rata nol (zero mean) pada model GARCH penting karena GARCH digunakan untuk memodelkan variabilitas atau fluktuasi volatilitas dalam data time series, bukan perubahan dalam nilai rata-rata. Dengan asumsi nilai rata-rata nol, model GARCH fokus pada menggambarkan bagaimana volatilitas data berubah seiring waktu, bukan menggambarkan perubahan nilai rata-rata data.
Namun, ada situasi di mana kita mungkin ingin mempertimbangkan model yang mengizinkan rata-rata non-zero. Ini mungkin terjadi jika data kita memiliki tren yang signifikan atau nilai rata-rata yang tidak nol secara alami. Dalam kasus seperti itu, kita dapat mempertimbangkan model GARCH dengan komponen yang menggambarkan rata-rata non-zero. Salah satu model yang dapat kita pertimbangkan adalah model GARCH-in-Mean (GARCH-M), yang memungkinkan kita memodelkan baik volatilitas kondisional maupun perubahan dalam nilai rata-rata data.
Model GARCH-M umumnya memiliki bentuk persamaan berikut:
\[ y_t=\mu+\varepsilon_t \]
\[ \varepsilon_t=\sigma_t.z_t \]
\[ \sigma_t^2=\omega+\sum_{i=1}^p\alpha_i.\varepsilon_{t-1}^2+\sum_{j=1}^q\beta_j.\sigma_{t-j}^2 \]
Di sini, \(\mu\) adalah rata-rata non-zero, dan \(\varepsilon_t\) adalah residula yang tergantung pada volatilitas kondisional. Model ini memungkinkan kita untuk memodelkan bagaimana volatilitas kondisional bergantung pada nilai-nilai residual sebelumnya dan juga memungkinkan kita untuk menggambarkan perubahan dalam rata-rata data.
Identifikasi Model
Identifikasi Model ARCH:
-
Periksa Plot ACF Residu:
Pertama, periksa plot ACF dari residu model. Plot ACF digunakan untuk mengidentifikasi pola ketergantungan serial dalam residu.
Jika terdapat pola signifikan pada lag-lag yang berdekatan (lag dengan nilai t yang lebih kecil), maka ini bisa menjadi indikasi adanya heteroskedastisitas kondisional dalam data. Pola ini mungkin berarti bahwa varians residual bergantung pada nilai-nilai residual sebelumnya.
-
Periksa Plot PACF Residu:
Selanjutnya, periksa plot PACF dari residu model. Plot PACF digunakan untuk mengidentifikasi jumlah lag dalam model AR(p) (autoregresi) yang mungkin diperlukan untuk mengatasi heteroskedastisitas.
Jika terdapat pola signifikan pada lag-lag yang berdekatan pada plot PACF, ini bisa menunjukkan adanya ketergantungan serial dalam residu dan indikasi bahwa model AR(p) mungkin diperlukan untuk mengatasi heteroskedastisitas.
-
Model ARCH yang Mungkin Cocok:
- Jika plot ACF dan PACF menunjukkan pola ketergantungan serial yang signifikan dalam residu, maka model ARCH(p) mungkin cocok, di mana \(p\)adalah tingkat autoregresi yang sesuai dengan pola pada PACF.
Identifikasi Model GARCH:
-
Periksa Plot ACF Residu Kuadrat:
Pertama, periksa plot ACF dari residu kuadrat. Plot ini akan memberikan gambaran tentang apakah varians residual bergantung pada nilai-nilai residual sebelumnya.
Jika terdapat pola signifikan pada lag-lag yang berdekatan pada plot ACF dari residu kuadrat, ini bisa menjadi indikasi adanya heteroskedastisitas kondisional dalam data, yang sesuai dengan model GARCH.
-
Periksa Plot PACF Residu Kuadrat:
Selanjutnya, periksa plot PACF dari residu kuadrat. Plot PACF digunakan untuk mengidentifikasi jumlah lag dalam model MA(q) (moving average) yang mungkin diperlukan untuk mengatasi heteroskedastisitas.
Jika terdapat pola signifikan pada lag-lag yang berdekatan pada plot PACF dari residu kuadrat, ini bisa menunjukkan adanya ketergantungan serial dalam varians residual dan indikasi bahwa model MA(q) mungkin diperlukan untuk mengatasi heteroskedastisitas.
-
Model GARCH yang Mungkin Cocok:
- Jika plot ACF dan PACF dari residu kuadrat menunjukkan pola ketergantungan serial yang signifikan, maka model GARCH(p, q) mungkin sesuai, di mana \(p\) adalah tingkat autoregresi untuk varians (\(\sigma^2\)) dan \(q\) adalah tingkat moving average untuk varians.
Penting untuk dicatat bahwa interpretasi plot ACF dan PACF dapat menjadi subjektif, dan hasilnya dapat bervariasi tergantung pada kompleksitas data time series kita. Oleh karena itu, langkah ini sering disertai dengan uji statistik formal dan evaluasi goodness-of-fit model untuk memastikan kecocokan model dengan data.
Ilustrasi
###installing and loading multiple packages
list.packages<-c("fGarch", "PerformanceAnalytics","rugarch","tseries","xts","FinTS", "urca")
new.packages <- list.packages[!(list.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
#Loading Packages
invisible(lapply(list.packages, require, character.only = TRUE))BMW<-read.csv("BMW.csv", sep = ";")
head(BMW)Calculate Financial Returns and remove first NA value
#set timeseries element
BMW$Years <- as.Date(BMW$Years, "%d/%m/%Y")
class(BMW$Years)## [1] "Date"BMW.z = zoo(x=BMW$Prices, order.by=BMW$Years)
#calculate log returns and remove first NA value
Return.BMW<-Return.calculate(BMW.z, method = "log")[-1]
head(Return.BMW)## 1995-11-07 1995-11-08 1995-11-09 1995-11-10 1995-11-13 1995-11-14
## 0.001334223 0.013080740 0.012911841 -0.007743447 0.000000000 0.014782595ts.plot(BMW$Prices, type="l", ylab="Price of BMW", col="#FFC6AC")
ts.plot(Return.BMW, type="l", ylab="Return BMW", col="#FFC6AC")
Periksa Kestasioneran
Menggunakan package urca
#apply ADF test with drift
ADF_Returns = ur.df(Return.BMW, type = "drift",selectlags = "AIC" )
#summary of he test
summary(ADF_Returns)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.156291 -0.010786 -0.000023 0.010886 0.138128
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0004118 0.0003170 1.299 0.193978
## z.lag.1 -1.0096970 0.0201650 -50.072 < 2e-16 ***
## z.diff.lag 0.0552591 0.0145795 3.790 0.000152 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02171 on 4692 degrees of freedom
## Multiple R-squared: 0.48, Adjusted R-squared: 0.4798
## F-statistic: 2166 on 2 and 4692 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -50.0717 1253.588
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78adf.test(Return.BMW, alternative = "stationary")## Warning in adf.test(Return.BMW, alternative = "stationary"): p-value smaller
## than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Return.BMW
## Dickey-Fuller = -17.905, Lag order = 16, p-value = 0.01
## alternative hypothesis: stationaryDapat dilihat bahwa nilai statistik uji -50.0717 lebih
kecil dari tiga critical values (tau2 pada tabel), dan p-value sebesar
0.01 maka dapat dikatakan bahwa nilai return BMW adalah
stasioner.
Periksa Volatility
# plot returns with squared and absolute returns
dataToPlot = cbind(Return.BMW, Return.BMW^2, abs(Return.BMW))
colnames(dataToPlot) = c("Returns", "Returns^2", "abs(Returns)")
plot.zoo(dataToPlot, main="BMW Daily Returns", col="#FFC6AC")
Periksa Normalitas
#histograme with normal density curve
hist(Return.BMW,prob=T,breaks=50,xlab="Daily Returns",main = "Histogram of Daily S&P 500 Returns and the Normal Distribution",
ylab="Probabillty Distribution",col="#FFC6AC")
mu<-mean(Return.BMW)
sigma<-sd(Return.BMW)
x<-seq(min(Return.BMW),max(Return.BMW),length=80)
y<-dnorm(x,mu,sigma)
lines(x,y,lwd=2,col="#999B84") 
##QQ-plot
qqnorm(Return.BMW, main = "BMW Daily Returns -QQ Plot", col = "#999B84")
qqline(Return.BMW)
#summary table
table.Stats(cbind(Return.BMW))#conduct Jarque-Bera test for normality
jarque.bera.test(Return.BMW)##
## Jarque Bera Test
##
## data: Return.BMW
## X-squared = 3741.4, df = 2, p-value < 2.2e-16Berdasarkan Jarque Bera Test, karena p-value < \(\alpha\) maka tolak \(H_0\) yang artiya data tidak berdistribusi normal.
Periksa ARCH Effect
## Convert to xts for time series features
Return.BMW<-as.xts(Return.BMW)
# plot autocorrelations of returns, returns^2 and abs(returns)
options(repr.plot.width=15, repr.plot.height=5)
par(mfrow=c(1,3))
acf(Return.BMW, main="BMW Returns",cex.main=10)
acf(Return.BMW^2, main="BMW Returns^2",cex.main=10)
acf(abs(Return.BMW), main="BMW abs(Returns)",cex.main=10)
par(mfrow=c(1,1))Plot nilai Returns^2 dan abs(Returns) menunjukkan bahwa terdapat
autokorelasi yang tinggi. Kita akan periksa juga dengan
Ljung-Box test.
# use Ljung Box.test from stats package to check auto correlation in squre retruns
Box.test(coredata(Return.BMW^2), type="Ljung-Box", lag = 12)##
## Box-Ljung test
##
## data: coredata(Return.BMW^2)
## X-squared = 1222.9, df = 12, p-value < 2.2e-16Nilai p-value < \(\alpha\) maka tolak \(H_0\) yang artinya terdapat autokorelasi pada data.
#ARCH LM Test
ArchTest(Return.BMW)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: Return.BMW
## Chi-squared = 500.01, df = 12, p-value < 2.2e-16Berdasarkan ARCH-LM Test dapat dilihat p-value < \(\alpha\) sehingga mengindikasikan tolak \(H_0\) : Tidak ada ARCH Effect. Maka data BMW Returns menunjukkan adanya ARCH Effect.
ARCH (1) Norm
ARCH Model estimation with normal innovation
#Specify the model
spec = ugarchspec(variance.model=list(garchOrder=c(1,0)),
mean.model=list(armaOrder=c(0,0)),distribution.model="norm",)
#Fit ARCH Model
arch11.fit=ugarchfit(data=Return.BMW,spec=spec, out.sample = 100)
#fitted model outcome
arch11.fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,0)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.035813 0.000007 5023.56091 0.00000
## omega 0.000000 0.000000 0.22024 0.82568
## alpha1 0.937775 0.000394 2380.30724 0.00000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.035813 0.068103 0.525865 0.59898
## omega 0.000000 0.003829 0.000006 0.99999
## alpha1 0.937775 7.011710 0.133744 0.89360
##
## LogLikelihood : -13197.03
##
## Information Criteria
## ------------------------------------
##
## Akaike 5.7429
## Bayes 5.7471
## Shibata 5.7429
## Hannan-Quinn 5.7444
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.864 0.1722
## Lag[2*(p+q)+(p+q)-1][2] 1.897 0.2804
## Lag[4*(p+q)+(p+q)-1][5] 2.238 0.5634
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.04251 0.8366
## Lag[2*(p+q)+(p+q)-1][2] 0.04869 0.9573
## Lag[4*(p+q)+(p+q)-1][5] 0.09444 0.9983
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 0.01234 0.500 2.000 0.9116
## ARCH Lag[4] 0.05724 1.397 1.611 0.9925
## ARCH Lag[6] 0.08157 2.222 1.500 0.9994
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 15.6575
## Individual Statistics:
## mu 0.7137
## omega 14.3888
## alpha1 0.7143
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 0.846 1.01 1.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 8.395 6.112e-17 ***
## Negative Sign Bias 4.157 3.283e-05 ***
## Positive Sign Bias 6.972 3.560e-12 ***
## Joint Effect 122.344 2.413e-26 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 7435 0
## 2 30 8221 0
## 3 40 8886 0
## 4 50 9601 0
##
##
## Elapsed time : 2.221207coef(arch11.fit)## mu omega alpha1
## 3.581296e-02 2.275579e-08 9.377750e-01#Conditional Volatility Plot for ARCH Model
plot.ts(sigma(arch11.fit), ylab="sigma(t)", col="#FFC6AC", main = "Conditional Varaince of ARCH Model")
GARCH (1,1) Norm
GARCH Model estimation with normal innovation
#Specify GARCH models:
spec = ugarchspec(variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(0,0)),
distribution.model="norm",)
#fit GARCH model
garch.fit=ugarchfit(data=Return.BMW,spec=spec, out.sample = 100)
##summary of GARCH fit
garch.fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000750 0.000234 3.2025 0.001363
## omega 0.000002 0.000002 1.1134 0.265526
## alpha1 0.071842 0.013647 5.2644 0.000000
## beta1 0.926656 0.013767 67.3076 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000750 0.000293 2.55737 0.010547
## omega 0.000002 0.000013 0.16352 0.870107
## alpha1 0.071842 0.086677 0.82885 0.407192
## beta1 0.926656 0.089064 10.40438 0.000000
##
## LogLikelihood : 11722.32
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.0982
## Bayes -5.0927
## Shibata -5.0982
## Hannan-Quinn -5.0963
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 11.26 0.0007933
## Lag[2*(p+q)+(p+q)-1][2] 11.95 0.0005841
## Lag[4*(p+q)+(p+q)-1][5] 15.02 0.0004381
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 3.707 0.05417
## Lag[2*(p+q)+(p+q)-1][5] 4.286 0.22047
## Lag[4*(p+q)+(p+q)-1][9] 6.784 0.21797
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.3644 0.500 2.000 0.5461
## ARCH Lag[5] 0.3885 1.440 1.667 0.9159
## ARCH Lag[7] 3.0703 2.315 1.543 0.5004
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 23.2576
## Individual Statistics:
## mu 0.11020
## omega 7.54575
## alpha1 0.07893
## beta1 0.06899
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 3.219 1.295e-03 ***
## Negative Sign Bias 3.727 1.964e-04 ***
## Positive Sign Bias 2.904 3.701e-03 ***
## Joint Effect 22.326 5.578e-05 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 188.1 7.744e-30
## 2 30 222.7 9.243e-32
## 3 40 241.6 4.745e-31
## 4 50 263.2 4.259e-31
##
##
## Elapsed time : 1.370477Berdasarkan uji Ljung-Box dapat dilihat bahwa tidak terdapat autokorelasi dalam residual kuadrat dari GARCH model.
#Conditional Volatility Plot for GARCH Model
plot.ts(sigma(garch.fit), ylab="sigma(t)", col="#FFC6AC", main = "Conditional Varaince of GARCH Model")
ARCH(1) dan GARCH(1,1) dengan distribusi normal masih menjadi model yang kurang sesuai.
QQ-plot of GARCH-normal fitted values
##QQ-plot for fitted GARCH Model
plot(garch.fit, which=9)
Selanjutnya akan dilakukan pendugaan model ARCH dan GARCH dengan
Non-Normal Innovations menggunakan rugarch() Package
ARCH(1) t-dist
ARCH Model with t-distribution innovation
#Specify the model
spec = ugarchspec(variance.model=list(garchOrder=c(1,0)),
mean.model=list(armaOrder=c(0,0)),distribution.model="std")
#Fit ARCH Model
arch11.fit.t=ugarchfit(data=Return.BMW,spec=spec, out.sample=100)
arch11.fit.t##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,0)
## Mean Model : ARFIMA(0,0,0)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000569 0.000251 2.2663 0.023434
## omega 0.000366 0.000021 17.1357 0.000000
## alpha1 0.374801 0.047250 7.9323 0.000000
## shape 3.744069 0.238705 15.6849 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000569 0.000225 2.5254 0.011557
## omega 0.000366 0.000028 13.0066 0.000000
## alpha1 0.374801 0.067346 5.5653 0.000000
## shape 3.744069 0.275759 13.5773 0.000000
##
## LogLikelihood : 11481.02
##
## Information Criteria
## ------------------------------------
##
## Akaike -4.9933
## Bayes -4.9877
## Shibata -4.9933
## Hannan-Quinn -4.9913
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 10.01 0.0015608
## Lag[2*(p+q)+(p+q)-1][2] 11.51 0.0007611
## Lag[4*(p+q)+(p+q)-1][5] 16.28 0.0002009
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 6.794 9.145e-03
## Lag[2*(p+q)+(p+q)-1][2] 55.155 2.776e-15
## Lag[4*(p+q)+(p+q)-1][5] 131.922 0.000e+00
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 96.64 0.500 2.000 0
## ARCH Lag[4] 145.73 1.397 1.611 0
## ARCH Lag[6] 195.64 2.222 1.500 0
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 3.6398
## Individual Statistics:
## mu 0.1172
## omega 2.7804
## alpha1 0.8572
## shape 2.1391
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 2.0340 0.04201 **
## Negative Sign Bias 0.9717 0.33127
## Positive Sign Bias 2.5525 0.01073 **
## Joint Effect 7.5827 0.05547 *
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 127.3 4.696e-18
## 2 30 177.8 2.544e-23
## 3 40 278.0 7.908e-38
## 4 50 398.3 3.169e-56
##
##
## Elapsed time : 5.776982#Conditional Volatility Plot for ARCH Model
plot.ts(sigma(arch11.fit.t), ylab="sigma(t)", col="#999B84", main = "Conditional Variance of ARCH Model")
GARCH(1,1) t-dist
GARCH Model with t-distribution Innovation
#Specify GARCH models:
spec = ugarchspec(variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(0,0)),
distribution.model="std",)
#fit GARCH model
garch.fit.t=ugarchfit(data=Return.BMW,spec=spec, out.sample = 100)
garch.fit.t##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000587 0.000218 2.6976 0.006984
## omega 0.000001 0.000001 1.1809 0.237652
## alpha1 0.062468 0.009141 6.8338 0.000000
## beta1 0.936532 0.009528 98.2932 0.000000
## shape 6.669334 0.608116 10.9672 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000587 0.000209 2.8083 0.004981
## omega 0.000001 0.000005 0.3161 0.751925
## alpha1 0.062468 0.041465 1.5065 0.131934
## beta1 0.936532 0.041605 22.5100 0.000000
## shape 6.669334 0.827899 8.0557 0.000000
##
## LogLikelihood : 11820.85
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.1407
## Bayes -5.1337
## Shibata -5.1407
## Hannan-Quinn -5.1382
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 10.91 0.0009574
## Lag[2*(p+q)+(p+q)-1][2] 11.59 0.0007246
## Lag[4*(p+q)+(p+q)-1][5] 14.61 0.0005638
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 6.652 0.009902
## Lag[2*(p+q)+(p+q)-1][5] 7.594 0.037133
## Lag[4*(p+q)+(p+q)-1][9] 10.651 0.036284
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.08793 0.500 2.000 0.7668
## ARCH Lag[5] 0.12511 1.440 1.667 0.9820
## ARCH Lag[7] 3.41505 2.315 1.543 0.4388
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 71.7215
## Individual Statistics:
## mu 0.07115
## omega 19.05071
## alpha1 0.20396
## beta1 0.11069
## shape 0.13538
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 3.488 4.906e-04 ***
## Negative Sign Bias 4.197 2.755e-05 ***
## Positive Sign Bias 2.856 4.311e-03 ***
## Joint Effect 25.772 1.064e-05 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 123.6 2.350e-17
## 2 30 145.0 2.225e-17
## 3 40 196.1 8.168e-23
## 4 50 181.8 3.829e-17
##
##
## Elapsed time : 2.800889##QQ-plot for fitted GARCH Model with t distribution
plot(garch.fit.t, which=9)
We can observe that the QQ-plot for Garch (1,1) does give us better results.
Various plot of fitted Model
#various plots for fitted values
plot(garch.fit.t, which= "all")##
## please wait...calculating quantiles...
Model Selection
Model Selection using Information Criterion
#Model selection usinh information critarion
model.list = list("arch(1,1)" = arch11.fit,
"arch(1,1).t" = arch11.fit.t,
"garch(1,1)" = garch.fit,
"garch(1,1).t" = garch.fit.t)
info.mat = sapply(model.list, infocriteria)
rownames(info.mat) = rownames(infocriteria(garch.fit))
info.mat## arch(1,1) arch(1,1).t garch(1,1) garch(1,1).t
## Akaike 5.742888 -4.993266 -5.098248 -5.140679
## Bayes 5.747086 -4.987668 -5.092650 -5.133682
## Shibata 5.742887 -4.993267 -5.098249 -5.140682
## Hannan-Quinn 5.744365 -4.991296 -5.096277 -5.138217Model GARCH(1,1) t-dist dipilih menjadi model terbaik karena memiliki nilai information criterion terkecil.
Forecasting
forc = ugarchforecast(garch.fit.t, data = Return.BMW, n.ahead = 10, n.roll =10)
print(forc)##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: sGARCH
## Horizon: 10
## Roll Steps: 10
## Out of Sample: 10
##
## 0-roll forecast [T0=2013-06-19]:
## Series Sigma
## T+1 0.000587 0.01224
## T+2 0.000587 0.01229
## T+3 0.000587 0.01234
## T+4 0.000587 0.01240
## T+5 0.000587 0.01245
## T+6 0.000587 0.01250
## T+7 0.000587 0.01255
## T+8 0.000587 0.01261
## T+9 0.000587 0.01266
## T+10 0.000587 0.01271plot(forc, which= "all")