PEMODELAN ARIMAX PADA DATA IHSG PERIODE APRIL 2020-DESEMBER 2022 MENGGUNAKAN PEUBAH PREDIKTOR KURS DOLLAR

Packages

knitr::opts_chunk$set(dev = "png",
                      dpi = 300,
                      echo = TRUE,
                      cache = TRUE)
library(readxl)
library(TTR)

## Warning: package 'TTR' was built under R version 4.3.1

library(forecast)

## Warning: package 'forecast' was built under R version 4.3.1

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(tseries)

## Warning: package 'tseries' was built under R version 4.3.1

library(TSA)

## Warning: package 'TSA' was built under R version 4.3.1

## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

library(dynlm)

## Warning: package 'dynlm' was built under R version 4.3.1

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(lmtest)
library(imputeTS)

## Warning: package 'imputeTS' was built under R version 4.3.1

## 
## Attaching package: 'imputeTS'

## The following object is masked from 'package:zoo':
## 
##     na.locf

## The following object is masked from 'package:tseries':
## 
##     na.remove

library(stats)
library(MASS)

## Warning: package 'MASS' was built under R version 4.3.1

library(padr)

## Warning: package 'padr' was built under R version 4.3.2

library(astsa)

## Warning: package 'astsa' was built under R version 4.3.1

## 
## Attaching package: 'astsa'

## The following object is masked from 'package:forecast':
## 
##     gas

library(tfarima)

## Warning: package 'tfarima' was built under R version 4.3.2

## 
## Attaching package: 'tfarima'

## The following object is masked from 'package:TSA':
## 
##     spec

## The following objects are masked from 'package:forecast':
## 
##     easter, seasadj

library(FinTS)

## Warning: package 'FinTS' was built under R version 4.3.2

## 
## Attaching package: 'FinTS'

## The following object is masked from 'package:forecast':
## 
##     Acf

library(rio)

## Warning: package 'rio' was built under R version 4.3.1

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.2

library(tsibble)

## Warning: package 'tsibble' was built under R version 4.3.1

## 
## Attaching package: 'tsibble'

## The following object is masked from 'package:zoo':
## 
##     index

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

library(tseries)
library(MASS)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

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

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

library(TTR)
library(forecast)
library(lmtest)
library(orcutt)

## Warning: package 'orcutt' was built under R version 4.3.1

library(HoRM)

## Warning: package 'HoRM' was built under R version 4.3.1

library(dLagM)

## Warning: package 'dLagM' was built under R version 4.3.1

## Loading required package: nardl

## Warning: package 'nardl' was built under R version 4.3.1

## 
## Attaching package: 'nardl'

## The following object is masked from 'package:FinTS':
## 
##     ArchTest

## 
## Attaching package: 'dLagM'

## The following object is masked from 'package:forecast':
## 
##     forecast

library(dynlm)
library(MLmetrics)

## Warning: package 'MLmetrics' was built under R version 4.3.1

## 
## Attaching package: 'MLmetrics'

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

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

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

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

## The following object is masked from 'package:tfarima':
## 
##     S

library(TSA)
library(aTSA)

## 
## Attaching package: 'aTSA'

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

## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test

## The following object is masked from 'package:forecast':
## 
##     forecast

## The following object is masked from 'package:graphics':
## 
##     identify

Data

data <- import("https://raw.githubusercontent.com/divanm/mpdw/main/Data/DATA%20IDX%20COMPOSITE%20JKSE.csv")
data

##           Date     Open     High      Low    Close Adj Close     Volume
## 1   2020-04-06 4623.429 4975.536 4562.902 4649.079  4649.079  255758500
## 2   2020-04-13 4649.079 4747.725 4480.607 4634.821  4634.821  218419900
## 3   2020-04-20 4634.821 4669.542 4441.090 4496.064  4496.064  238397200
## 4   2020-04-27 4496.064 4726.774 4474.893 4716.403  4716.403  183906700
## 5   2020-05-04 4716.403 4716.403 4576.228 4597.430  4597.430  179791700
## 6   2020-05-11 4597.430 4659.862 4460.273 4507.607  4507.607  231060000
## 7   2020-05-18 4507.607 4609.042 4487.179 4545.952  4545.952  139287000
## 8   2020-05-25 4545.952 4755.957 4541.834 4753.612  4753.612  288083800
## 9   2020-06-01 4753.612 5014.764 4753.612 4947.782  4947.782  372565600
## 10  2020-06-08 4947.782 5139.406 4712.073 4880.359  4880.359  490084700
## 11  2020-06-15 4880.359 5018.985 4809.619 4942.275  4942.275  336775500
## 12  2020-06-22 4942.275 4977.651 4879.133 4904.088  4904.088  260206800
## 13  2020-06-29 4904.088 4997.482 4862.041 4973.794  4973.794  249816600
## 14  2020-07-06 4973.794 5111.564 4973.490 5031.256  5031.256  358094500
## 15  2020-07-13 5031.256 5116.464 5024.934 5079.585  5079.585  349862400
## 16  2020-07-20 5079.585 5162.980 5031.888 5082.991  5082.991  427056800
## 17  2020-07-27 5082.991 5149.630 5080.121 5149.627  5149.627  356913600
## 18  2020-08-03 5149.627 5187.964 4928.468 5143.893  5143.893  423927900
## 19  2020-08-10 5157.834 5279.349 5157.834 5247.690  5247.690  340206700
## 20  2020-08-17 5247.690 5327.316 5243.995 5272.810  5272.810  188318800
## 21  2020-08-24 5272.810 5381.948 5261.408 5346.659  5346.659  665155800
## 22  2020-08-31 5346.659 5369.447 5188.609 5239.851  5239.851  568186400
## 23  2020-09-07 5235.010 5256.305 4754.799 5016.712  5016.712  492882200
## 24  2020-09-14 5060.024 5187.279 5013.186 5059.223  5059.223  503336800
## 25  2020-09-21 5058.990 5075.819 4820.331 4945.791  4945.791  406370800
## 26  2020-09-28 4962.953 4991.762 4841.362 4926.734  4926.734  440809400
## 27  2020-10-05 4947.040 5057.836 4915.685 5053.663  5053.663  421591900
## 28  2020-10-12 5078.129 5182.531 5064.171 5103.414  5103.414  525143900
## 29  2020-10-19 5116.758 5135.078 5063.697 5112.188  5112.188  522912300
## 30  2020-10-26 5113.149 5158.265 5110.618 5128.225  5128.225  210436500
## 31  2020-11-02 5108.029 5335.529 5073.500 5335.529  5335.529  516535700
## 32  2020-11-09 5371.971 5520.908 5319.420 5461.058  5461.058  669558400
## 33  2020-11-16 5500.034 5628.442 5462.460 5571.656  5571.656  789934400
## 34  2020-11-23 5583.330 5795.836 5583.330 5783.335  5783.335 1255563200
## 35  2020-11-30 5779.671 5853.156 5563.862 5810.483  5810.483 1160743700
## 36  2020-12-07 5854.303 6004.423 5854.303 5938.329  5938.329 1003237500
## 37  2020-12-14 5959.275 6160.977 5959.275 6104.324  6104.324 1292048600
## 38  2020-12-21 6131.623 6195.154 5853.261 6008.709  6008.709  692226400
## 39  2020-12-28 6067.000 6143.870 5962.010 5979.073  5979.073  526030000
## 40  2021-01-04 5997.830 6275.742 5929.048 6257.835  6257.835 1007037000
## 41  2021-01-11 6278.413 6472.312 6278.374 6373.412  6373.412 1527600500
## 42  2021-01-18 6365.026 6504.992 6283.313 6307.127  6307.127 1081229300
## 43  2021-01-25 6322.519 6322.727 5825.292 5862.352  5862.352  813369600
## 44  2021-02-01 5856.777 6179.368 5735.469 6151.729  6151.729  872980000
## 45  2021-02-08 6193.551 6286.291 6157.135 6222.521  6222.521  563483000
## 46  2021-02-15 6244.385 6314.555 6173.594 6231.932  6231.932  739117700
## 47  2021-02-22 6267.025 6312.873 6184.517 6241.796  6241.796  880483600
## 48  2021-03-01 6281.857 6394.452 6245.308 6258.749  6258.749  987879800
## 49  2021-03-08 6304.000 6364.357 6167.718 6358.209  6358.209  652024500
## 50  2021-03-15 6374.762 6387.736 6268.839 6356.160  6356.160  742809300
## 51  2021-03-22 6346.795 6354.939 6058.843 6195.562  6195.562  712370600
## 52  2021-03-29 6205.569 6230.990 5892.645 6011.456  6011.456  462735200
## 53  2021-04-05 6040.055 6113.156 5944.421 6070.209  6070.209  720785500
## 54  2021-04-12 6080.792 6115.620 5883.524 6086.258  6086.258  669695000
## 55  2021-04-19 6084.785 6096.995 5973.247 6016.864  6016.864  675194300
## 56  2021-04-26 6018.340 6033.900 5950.868 5995.616  5995.616  641211400
## 57  2021-05-03 5999.661 6005.088 5922.493 5928.309  5928.309  667078400
## 58  2021-05-10 5942.106 5985.372 5911.376 5938.351  5938.351  240955800
## 59  2021-05-17 5950.137 5958.793 5742.038 5773.120  5773.120  696962000
## 60  2021-05-24 5793.067 5904.838 5759.340 5848.616  5848.616  679387300
## 61  2021-05-31 5869.218 6103.859 5860.536 6065.166  6065.166  687508700
## 62  2021-06-07 6075.507 6134.882 5972.449 6095.497  6095.497 1116799300
## 63  2021-06-14 6110.797 6124.872 5944.053 6007.120  6007.120  955804200
## 64  2021-06-21 5959.965 6130.096 5884.918 6022.399  6022.399  851733600
## 65  2021-06-28 6015.279 6043.432 5913.591 6023.008  6023.008  825550700
## 66  2021-07-05 6024.165 6080.222 5985.354 6039.844  6039.844  754205300
## 67  2021-07-12 6060.054 6114.286 5947.618 6072.510  6072.510  716924100
## 68  2021-07-19 6062.966 6166.305 6015.149 6101.690  6101.690  597748400
## 69  2021-07-26 6109.031 6160.463 6068.512 6070.039  6070.039  869133600
## 70  2021-08-02 6098.008 6263.539 6048.101 6203.431  6203.431 1228093800
## 71  2021-08-09 6203.747 6239.015 6042.475 6139.492  6139.492  869347300
## 72  2021-08-16 6144.943 6147.304 5938.407 6030.772  6030.772  796085900
## 73  2021-08-23 6037.700 6138.497 6021.954 6041.366  6041.366  977033500
## 74  2021-08-30 6062.113 6169.056 6054.281 6126.921  6126.921  893890700
## 75  2021-09-06 6137.126 6150.358 5982.766 6094.873  6094.873  958513400
## 76  2021-09-13 6080.062 6137.669 6052.970 6133.246  6133.246  992844400
## 77  2021-09-20 6132.034 6163.884 5996.407 6144.815  6144.815  950159400
## 78  2021-09-27 6141.470 6286.943 6086.263 6228.845  6228.845 1102928200
## 79  2021-10-04 6234.799 6497.656 6234.799 6481.769  6481.769 1172175100
## 80  2021-10-11 6482.957 6680.009 6451.678 6633.338  6633.338  941178200
## 81  2021-10-18 6649.281 6687.134 6585.557 6643.738  6643.738  757917300
## 82  2021-10-25 6626.367 6680.117 6509.878 6591.346  6591.346  899190300
## 83  2021-11-01 6618.122 6627.851 6480.010 6581.785  6581.785  828841100
## 84  2021-11-08 6599.405 6714.158 6592.055 6651.054  6651.054  921099300
## 85  2021-11-15 6661.849 6720.988 6592.228 6720.263  6720.263 1134263700
## 86  2021-11-22 6731.515 6754.464 6544.896 6561.553  6561.553 1087199900
## 87  2021-11-29 6552.803 6647.477 6484.578 6538.506  6538.506  795330100
## 88  2021-12-06 6552.816 6652.922 6525.985 6652.922  6652.922 1078695000
## 89  2021-12-13 6663.739 6688.379 6559.303 6601.932  6601.932 1011838100
## 90  2021-12-20 6571.345 6592.167 6529.593 6562.900  6562.900  941359400
## 91  2021-12-27 6570.558 6621.579 6562.554 6581.482  6581.482  763376600
## 92  2022-01-03 6586.260 6738.110 6586.133 6701.316  6701.316  884588500
## 93  2022-01-10 6697.376 6727.767 6625.757 6693.401  6693.401  854569300
## 94  2022-01-17 6711.406 6726.373 6534.270 6726.373  6726.373  805658500
## 95  2022-01-24 6700.822 6712.262 6523.929 6645.511  6645.511  986322200
## 96  2022-01-31 6656.707 6731.391 6626.656 6731.391  6731.391  777116900
## 97  2022-02-07 6751.349 6874.351 6748.909 6815.607  6815.607 1238259600
## 98  2022-02-14 6791.316 6899.410 6698.509 6892.818  6892.818 1092845600
## 99  2022-02-21 6896.752 6929.911 6758.861 6888.171  6888.171 1127214800
## 100 2022-02-28 6964.697 6996.935 6861.959 6928.328  6928.328  769293800
## 101 2022-03-07 6905.324 6929.863 6814.183 6922.602  6922.602 1085654100
## 102 2022-03-14 6949.206 7032.702 6894.921 6954.965  6954.965  949839000
## 103 2022-03-21 6945.893 7055.341 6929.934 7002.532  7002.532 1086012100
## 104 2022-03-28 6996.318 7099.497 6987.223 7078.760  7078.760  890142300
## 105 2022-04-04 7082.247 7216.494 7060.214 7210.835  7210.835  986772300
## 106 2022-04-11 7212.003 7355.300 7146.586 7235.532  7235.532 1073108900
## 107 2022-04-18 7245.504 7297.394 7174.810 7225.606  7225.606 1118782600
## 108 2022-04-25 7204.161 7267.112 7121.860 7228.914  7228.914  805992800
## 109 2022-05-02       NA       NA       NA       NA        NA         NA
## 110 2022-05-09 7154.917 7156.484 6509.879 6597.993  6597.993  778348200
## 111 2022-05-16 6603.893 6965.106 6574.136 6918.144  6918.144  816107000
## 112 2022-05-23 6930.939 7032.822 6802.713 7026.256  7026.256  750357000
## 113 2022-05-30 7052.173 7233.995 6974.616 7182.961  7182.961  949041500
## 114 2022-06-06 7163.174 7257.997 7051.223 7086.648  7086.648 1281571500
## 115 2022-06-13 6992.230 7138.496 6882.644 6936.967  6936.967 1322143900
## 116 2022-06-20 6937.008 7067.747 6859.599 7042.937  7042.937 1133599100
## 117 2022-06-27 7043.018 7070.519 6777.317 6794.328  6794.328  840889500
## 118 2022-07-04 6794.368 6796.519 6559.637 6740.219  6740.219  810131700
## 119 2022-07-11 6740.019 6757.087 6607.702 6651.905  6651.905  794121300
## 120 2022-07-18 6652.065 6901.720 6611.926 6886.962  6886.962  853012300
## 121 2022-07-25 6887.042 7032.353 6857.801 6951.123  6951.123 1387470600
## 122 2022-08-01 6951.123 7090.767 6902.630 7084.655  7084.655 1102623400
## 123 2022-08-08 7084.655 7181.100 7021.673 7129.277  7129.277 1169007600
## 124 2022-08-15 7129.277 7230.111 7080.753 7172.434  7172.434  931624200
## 125 2022-08-22 7172.434 7210.163 7064.501 7135.248  7135.248 1241361100
## 126 2022-08-29 7135.088 7223.126 7015.347 7177.179  7177.179 1359158600
## 127 2022-09-05 7177.139 7287.700 7147.976 7242.656  7242.656 1535050500
## 128 2022-09-12 7242.696 7377.495 7163.073 7168.870  7168.870 1462367300
## 129 2022-09-19 7168.950 7252.186 7127.457 7178.583  7178.583 1266925400
## 130 2022-09-26 7178.503 7178.503 6926.860 7040.798  7040.798 1036700200
## 131 2022-10-03 7040.798 7135.908 6995.062 7026.783  7026.783 1020558200
## 132 2022-10-10 7026.663 7026.663 6814.530 6814.530  6814.530 1046699200
## 133 2022-10-17 6814.530 7058.911 6747.380 7017.771  7017.771  995211200
## 134 2022-10-24 7017.811 7108.816 7016.702 7056.040  7056.040  989950000
## 135 2022-10-31 7056.040 7128.136 6962.853 7045.527  7045.527  907843500
## 136 2022-11-07 7045.527 7111.440 6956.288 7089.206  7089.206  989493100
## 137 2022-11-14 7089.206 7104.764 6955.529 7082.181  7082.181  831349300
## 138 2022-11-21 7082.181 7108.831 7013.221 7053.150  7053.150  816522900
## 139 2022-11-28 7053.110 7090.277 6967.950 7019.639  7019.639  932257000
## 140 2022-12-05 7019.639 7053.903 6683.629 6715.118  6715.118  763168800
## 141 2022-12-12 6715.118 6854.098 6641.810 6812.193  6812.193 1047258100
## 142 2022-12-19 6812.154 6844.121 6715.045 6800.673  6800.673  668236600
## 143 2022-12-26 6800.713 6858.147 6796.167 6835.808  6835.808  117543200

Mengambil data close

data_close <- as.numeric(data$Close)
data_close <- as.data.frame(data_close)

Mengecek keberadaan missing value

data_close[which(is.na(data_close$data)),]

## [1] NA

Menduga missing value

data_close <- na_interpolation(data_close, option = "spline")
data_close[109,]

## [1] 6840.217

Mengubah data menjadi data time series

data.ts <- ts(data_close)

Eksplorasi Data

Plot Data Penuh

plot.ts(data.ts, lty=1, xlab="waktu", ylab="IHSG yahoo finance", main="Plot Data IHSG YAHOO FINANCE")

Berdasarkan plot data deret waktu, terlihat bahwa data cenderung memiliki trend yang naik. Berdasarkan pola data, pembagian data latih dan data uji ditetapkan dengan proporsi 90%:10%.

Pembagian Data

Data kemudian dibagi menjadi data latih dan data uji. Pembagian kali ini dilakukan dengan proporsi / perbandingan, yaitu 90:10.

data.train <- data_close$data[1:128]
train.ts <- ts(data.train)
data.test <- data_close$data[129:143]
test.ts <- ts(data.test)

Plot Data Latih

train.ts<-ts(data.train)
plot.ts(train.ts, lty=1, xlab="waktu", ylab="IHSG", main="Plot IHSG Train")

Berdasarkan plot data deret waktu pada data latih, terlihat bahwa data cenderung memiliki trend yang naik dan cenderung tidak bergerak pada nilai tengah tertentu. Hal ini mengindikasikan bahwa data tidak stasioner dalam rataan.

Plot Data Uji

test.ts<-ts(data.test, start = 129)
plot.ts(test.ts, lty=1, xlab="waktu", ylab="IHSG", main="Plot IHSG")

Uji Stasioneritas Data

Plot ACF

acf(train.ts)

Berdasarkan plot ACF, terlihat bahwa plot ACF data menurun secara perlahan (tails of slowly). Hal ini juga menjadi indikasi bahwa data tidak stasioner dalam rataan

Uji ADF

tseries::adf.test(train.ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -2.9063, Lag order = 5, p-value = 0.1998
## alternative hypothesis: stationary

\(H_0\) : Data tidak stasioner dalam rataan

\(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.3898 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa data tidak stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga ketidakstasioneran model kedepannya harus ditangani

Plot Box-Cox

index <- seq(1:128)
bc = boxcox(train.ts~index, lambda = seq(-5,10,by=1))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda

## [1] 3.030303

#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]

##  [1] 2.272727 2.424242 2.575758 2.727273 2.878788 3.030303 3.181818 3.333333
##  [9] 3.484848 3.636364 3.787879

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 3.030303 dan pada selang kepercayaan 95% nilai memiliki batas bawah 2.272727 dan batas atas 3.787879. Selang tersebut tidak memuat nilai satu sehingga dapat dikatakan bahwa data tidak stasioner dalam ragam.

Penanganan Ketidakstasioneran Data

Ketidakstasioner dalam Ragam

data.trans <- data_close^3
train.trans <- (data.train)^3
test.trans <- (data.test)^3
train.trans.ts <- ts(train.trans, frequency = 1)

Plot Box-Cox

index <- seq(1:128)
bc = boxcox(train.trans.ts~index, lambda = seq(-5,10,by=1))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda

## [1] 1.060606

#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]

## [1] 0.7575758 0.9090909 1.0606061 1.2121212

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 1.060606 dan pada selang kepercayaan 95% nilai memiliki batas bawah 0.7575758 dan batas atas 1.2121212. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data stasioner dalam ragam. ## Ketidakstasioneran dalam Rataan

train.diff<-diff(train.trans.ts,differences = 1) 
plot.ts(train.diff, lty=1, xlab="waktu", ylab="Data Difference 1 IHSG", main="Plot Difference IHSG")

Berdasarkan plot data deret waktu, terlihat bahwa data sudah stasioner dalam rataan ditandai dengan data bergerak pada nilai tengah tertentu (tidak terdapat trend ataupun musiman pada data) ### Plot ACF

acf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cuts off pada lag ke 2. Hal ini menandakan data sudah stasioner dalam rataan dan ketidakstasioneran data telah berhasil tertangani.

Uji ADF

tseries::adf.test(train.diff)

## Warning in tseries::adf.test(train.diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.diff
## Dickey-Fuller = -5.8203, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

\(H_0\) : Data tidak stasioner dalam rataan

\(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.01 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) atau data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga dalam hal ini ketidakstasioneran data sudah berhasil ditangani dan dapat dilanjutkan ke pemodelan

Identifikasi Model

Plot ACF

acf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung cuts off pada lag ke 2, sehingga jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,1,2).

Plot PACF

pacf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot PACF cenderung cuts off pada lag ke 2, sehingga jika plot ACF dianggap tails of, maka model tentatifnya adalah ARIMA(2,1,0).

Jika baik plot ACF maupun plot PACF keduanya dianggap tails of, maka model yang terbentuk adalah ARIMA(2,1,2)

Plot EACF

eacf(train.diff)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o x o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 o x o o o o o o o o o  o  o  o 
## 4 x o o o o o o o o o o  o  o  o 
## 5 o o x x x o o o o o o  o  o  o 
## 6 o x x o o o o o o o o  o  o  o 
## 7 o x x o o x o o o o o  o  o  o

Identifikasi model menggunakan plot EACF dilakukan dengan melihat ujung segitiga pada pola segitiga nol. Dalam hal ini model tentatif yang terbentuk adalah ARIMA(0,1,2), ARIMA(2,1,0), ARIMA(2,1,2), ARIMA(1,1,2), ARIMA(1,1,1).

Pendugaan Parameter Model Tentatif

ARIMA(0,1,2)

model1.da=Arima(train.trans.ts, order=c(0,1,2),method="ML")
summary(model1.da) #AIC=1412.54

## Series: train.trans.ts 
## ARIMA(0,1,2) 
## 
## Coefficients:
##          ma1      ma2
##       0.0897  -0.1417
## s.e.  0.0907   0.0926
## 
## sigma^2 = 1.999e+20:  log likelihood = -3147.5
## AIC=6301.01   AICc=6301.2   BIC=6309.54
## 
## Training set error measures:
##                      ME        RMSE        MAE      MPE     MAPE      MASE
## Training set 2204339325 13973403058 9770310874 0.884391 4.298203 0.9677074
##                     ACF1
## Training set -0.01059438

lmtest::coeftest(model1.da) #ma2 tidak signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ma1  0.089688   0.090693  0.9889   0.3227
## ma2 -0.141712   0.092616 -1.5301   0.1260

ARIMA(2,1,0)

model2.da=Arima(train.trans.ts, order=c(2,1,0),method="ML")
summary(model2.da) #AIC=1445.9

## Series: train.trans.ts 
## ARIMA(2,1,0) 
## 
## Coefficients:
##          ar1      ar2
##       0.1197  -0.1678
## s.e.  0.0875   0.0873
## 
## sigma^2 = 1.986e+20:  log likelihood = -3147.08
## AIC=6300.17   AICc=6300.36   BIC=6308.7
## 
## Training set error measures:
##                      ME        RMSE        MAE       MPE     MAPE      MASE
## Training set 2186591024 13926542369 9737228005 0.8799691 4.298029 0.9644307
##                     ACF1
## Training set -0.04075674

lmtest::coeftest(model2.da) #seluruh parameter signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ar1  0.119743   0.087465  1.3690  0.17099  
## ar2 -0.167805   0.087315 -1.9218  0.05463 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(2,1,2)

model3.da=Arima(train.trans.ts, order=c(2,1,2),method="ML")
summary(model3.da) #AIC=1411.04

## Series: train.trans.ts 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2
##       -0.0580  -0.9622  0.1558  1.0000
## s.e.   0.0309   0.0262  0.0230  0.0563
## 
## sigma^2 = 1.858e+20:  log likelihood = -3144.09
## AIC=6298.17   AICc=6298.67   BIC=6312.39
## 
## Training set error measures:
##                      ME        RMSE        MAE       MPE     MAPE      MASE
## Training set 1909943338 13361486252 9142375183 0.7811361 4.043717 0.9055131
##                       ACF1
## Training set -0.0002119683

lmtest::coeftest(model3.da) #ar1 dan ma2 tidak signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.057966   0.030872  -1.8776   0.06044 .  
## ar2 -0.962185   0.026196 -36.7300 < 2.2e-16 ***
## ma1  0.155758   0.023021   6.7660 1.324e-11 ***
## ma2  0.999995   0.056346  17.7475 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,2)

model4.da=Arima(train.trans.ts, order=c(1,1,2),method="ML")
summary(model4.da) #AIC=11412.85 

## Series: train.trans.ts 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.4029  -0.3037  -0.1868
## s.e.  0.3343   0.3265   0.0817
## 
## sigma^2 = 2e+20:  log likelihood = -3147
## AIC=6302.01   AICc=6302.34   BIC=6313.38
## 
## Training set error measures:
##                      ME        RMSE        MAE       MPE     MAPE      MASE
## Training set 2439474309 13917820932 9835551287 0.9830953 4.337995 0.9741692
##                     ACF1
## Training set -0.02615259

lmtest::coeftest(model4.da) #ar1 dan ma1 tidak signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ar1  0.402933   0.334315  1.2052  0.22811  
## ma1 -0.303735   0.326456 -0.9304  0.35216  
## ma2 -0.186791   0.081685 -2.2867  0.02221 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,1)

model5.da=Arima(train.trans.ts, order=c(1,1,1),method="ML")
summary(model5.da) #AIC=11412.85 

## Series: train.trans.ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1     ma1
##       -0.3566  0.4921
## s.e.   0.4545  0.4224
## 
## sigma^2 = 2.023e+20:  log likelihood = -3148.25
## AIC=6302.49   AICc=6302.69   BIC=6311.03
## 
## Training set error measures:
##                      ME        RMSE        MAE       MPE     MAPE      MASE
## Training set 1894292511 14056840875 9730600741 0.7597473 4.289273 0.9637743
##                     ACF1
## Training set -0.03627121

lmtest::coeftest(model5.da) #ar1 dan ma1 tidak signifikan

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ar1 -0.35656    0.45448 -0.7845   0.4327
## ma1  0.49207    0.42239  1.1650   0.2440

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dan juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(2,1,2).

Analisis Sisaan

Model terbaik hasil identifikasi kemudian dicek asumsi sisaannya. Sisaan model ARIMA harus memenuhi asumsi normalitas, kebebasan sisaan, dan kehomogenan ragam. Diagnostik model dilakukan secara eksplorasi dan uji formal.

Eksplorasi Sisaan

#Eksplorasi 
sisaan.da <- model3.da$residuals 
par(mfrow=c(1,1)) 
qqnorm(sisaan.da) 
qqline(sisaan.da, col = "blue", lwd = 2) 

plot(c(1:length(sisaan.da)),sisaan.da) 

acf(sisaan.da) 

pacf(sisaan.da) 

par(mfrow = c(1,1))

Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan menyebar tidak normal ditandai dengan titik titik yang cenderung tidak mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki ragam yang homogen. Plot ACF dan PACF sisaan ARIMA(2,1,2) juga signifikan pada lag 11 yang menandakan tidak saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.

Uji Formal

#1) Sisaan Menyebar Normal 
ks.test(sisaan.da,"pnorm")  #tak tolak H0 > sisaan menyebar normal

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.da
## D = 0.60156, p-value < 2.2e-16
## alternative hypothesis: two-sided

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS). Hipotesis pada uji KS adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value sebesar 2.2e-16 yang kurang dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Hal ini tidak sesuai dengan hasil eksplorasi menggunakan plot kuantil-kuantil normal.

#2) Sisaan saling bebas/tidak ada autokorelasi 
Box.test(sisaan.da, type = "Ljung")  #tak tolak H0 > sisaan saling bebas

## 
##  Box-Ljung test
## 
## data:  sisaan.da
## X-squared = 5.887e-06, df = 1, p-value = 0.9981

Selanjutnya akan dilakukan uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box. Hipotesis yang digunakan adalah sebagai berikut.

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar 0.9981 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas.

#3) Sisaan homogen 
Box.test((sisaan.da)^2, type = "Ljung")  #tak tolak H0 > sisaan homogen

## 
##  Box-Ljung test
## 
## data:  (sisaan.da)^2
## X-squared = 2.3331, df = 1, p-value = 0.1266

Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.1266 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa ragam sisaan homogen.

#4) Nilai tengah sisaan sama dengan nol 
t.test(sisaan.da, mu = 0, conf.level = 0.95)  #tak tolak h0 > nilai tengah sisaan sama dengan 0

## 
##  One Sample t-test
## 
## data:  sisaan.da
## t = 1.6276, df = 127, p-value = 0.1061
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -412131316 4232017992
## sample estimates:
##  mean of x 
## 1909943338

Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol. Hipotesis yang diujikan sebagai berikut.

\(H_0\) : nilai tengah sisaan sama dengan 0

\(H_1\) : nilai tengah sisaan tidak sama dengan 0

Berdasarkan uji-ttersebut, didapat p-value sebesar 0.1061 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol.

Overfitting

Tahapan selanjutnya adalah overfitting dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(2,1,2) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(3,1,2) dan ARIMA(2,1,3).

Model ARIMA(3,1,2)

model.overfit1=Arima(train.trans.ts, order=c(3,1,2),method="ML")
summary(model.overfit1) #1433.78

## Series: train.trans.ts 
## ARIMA(3,1,2) 
## 
## Coefficients:
##           ar1      ar2     ar3     ma1     ma2
##       -0.0374  -0.9613  0.0237  0.1566  0.9965
## s.e.   0.0925   0.0267  0.0916  0.0237  0.0549
## 
## sigma^2 = 1.88e+20:  log likelihood = -3144.06
## AIC=6300.12   AICc=6300.82   BIC=6317.19
## 
## Training set error measures:
##                      ME        RMSE        MAE       MPE     MAPE      MASE
## Training set 1867658219 13384782785 9148055105 0.7641803 4.051503 0.9060757
##                     ACF1
## Training set -0.01900881

lmtest::coeftest(model.overfit1) #semua parameter signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.037368   0.092478  -0.4041    0.6862    
## ar2 -0.961287   0.026685 -36.0240 < 2.2e-16 ***
## ar3  0.023702   0.091634   0.2587    0.7959    
## ma1  0.156607   0.023679   6.6139 3.744e-11 ***
## ma2  0.996476   0.054898  18.1515 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model ARIMA(2,1,3)

model.overfit2=Arima(train.trans.ts, order=c(2,1,3),method="ML")
summary(model.overfit2) #1409.25

## Series: train.trans.ts 
## ARIMA(2,1,3) 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2     ma3
##       -0.0606  -0.9633  0.1845  1.0046  0.0299
## s.e.   0.0321   0.0263  0.1033  0.0573  0.1047
## 
## sigma^2 = 1.873e+20:  log likelihood = -3144.05
## AIC=6300.09   AICc=6300.79   BIC=6317.16
## 
## Training set error measures:
##                      ME        RMSE        MAE      MPE     MAPE      MASE
## Training set 1855973345 13361268920 9133766916 0.759841 4.045633 0.9046605
##                    ACF1
## Training set -0.0228522

lmtest::coeftest(model.overfit2) #ar1 tidak signifikan

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1 -0.060565   0.032147  -1.8840  0.05956 .  
## ar2 -0.963335   0.026316 -36.6060  < 2e-16 ***
## ma1  0.184494   0.103273   1.7865  0.07402 .  
## ma2  1.004626   0.057337  17.5215  < 2e-16 ***
## ma3  0.029945   0.104746   0.2859  0.77497    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model yang dipilih adalah model lama, yaitu ARIMA(2,1,2) karena parameter ar3 pada model ARIMA(3,1,2) dan ARIMA(2,1,3) tidak berbeda signifikan dengan 0. Selain itu, juga parameter ar1 dan ar2 pada model ARIMA(3,1,2) dan ARIMA(2,1,3) tidak jauh berbeda pada ARIMA(2,1,2) sehingga tetap dipilih model awal. Pada kedua model overfitting juga tidak semua parameternya signifikan sehingga tetap dipilih model awal yaitu ARIMA(2,1,2)

Peramalan

Peramalan dilakukan menggunakan fungsi forecast() . Contoh peramalan

ramalan<- forecast::forecast(model3.da, h =15) 
ramalan

##     Point Forecast        Lo 80        Hi 80        Lo 95        Hi 95
## 129   363561556675 345969503698 381153609652 336656833387 390466279963
## 130   354218176162 328109540321 380326812003 314288463375 394147888949
## 131   359441792231 326788537118 392095047344 309502948005 409380636457
## 132   368129065574 330863908148 405394223001 311136925423 425121205725
## 133   362599413028 321324609218 403874216839 299475047914 425723778143
## 134   354561176118 308975930232 400146422005 284844558726 424277793511
## 135   360347670555 310782710171 409912630940 284544605068 436150736043
## 136   367746526000 315044845852 420448206147 287146261609 448346790390
## 137   361749964326 306093234262 417406694390 276630341085 446869587566
## 138   354978489146 296042233449 413914744843 264843264271 445113714022
## 139   361140807825 299120281447 423161334202 266288597995 455993017654
## 140   367299018601 302751544897 431846492304 268582176335 466015860867
## 141   361012759340 293993040438 428032478242 258514944321 463510574359
## 142   355451806966 285679842213 425223771720 248744795151 462158818782
## 143   361822699458 289462899516 434182499401 251157935254 472487463662

data.ramalan <- ramalan$mean
Hasil1 <- (data.ramalan)^(1/3) 

AKURASI

perbandingan <- matrix(data=c(test.ts, Hasil1),
                     nrow = length(test.ts), ncol = 2)
colnames(perbandingan) <- c("Aktual","Hasil Forecast")
perbandingan

##         Aktual Hasil Forecast
##  [1,] 7178.583       7137.169
##  [2,] 7040.798       7075.497
##  [3,] 7026.783       7110.108
##  [4,] 6814.530       7166.933
##  [5,] 7017.771       7130.867
##  [6,] 7056.040       7077.780
##  [7,] 7045.527       7116.076
##  [8,] 7089.206       7164.450
##  [9,] 7082.181       7125.295
## [10,] 7053.150       7080.556
## [11,] 7019.639       7121.293
## [12,] 6715.118       7161.543
## [13,] 6812.193       7120.451
## [14,] 6800.673       7083.701
## [15,] 6835.808       7125.772

A <- accuracy(ts(Hasil1), head(data.test, n=length(data.test)))
print(A)

##                 ME     RMSE      MAE       MPE     MAPE
## Test set -147.2995 204.8412 152.8213 -2.151847 2.228769

MAPE ARIMA(2,1,2) sudah sangat baik yaitu sebesar 2,22. Akan tetapi, akan dibandingkan dengan model ARIMAX dengan peubah kovariat nya adalah kurs dolar.

MODEL ARIMAX

Import data variabel x

data_kurs <- import("https://raw.githubusercontent.com/divanm/mpdw/main/Data/kurs_dollar.csv")
kurs<- data_kurs$Terakhir
kurs <- (kurs)^3

Ubah data time series

kurs.ts <- ts(kurs)

Eksplorasi Data

Plot time series

plot(kurs.ts,xlab ="Periode", ylab = "Kurs Dollar", col="black", main = "Plot Data Histori Kurs Dollar")
points(kurs.ts)

Split Data

kurs.train <- kurs.ts[1:128]
traink.ts <- ts(kurs.train)
kurs.test <- kurs.ts[129:143]
testk.ts <- ts(kurs.test)

Plot time series variabel x

plot(traink.ts,xlab = "Periode", ylab = "Data Kurs Dollar", col="black", main = "Data Training Histori Kurs Dollar")
points(traink.ts)

Korelasi Data

Korelasi Antar Data Test

korelasitest <- c("IHSG vs Kurs Dollar" = cor(test.trans, kurs.test))
kort <- data.frame("Nilai Korelasi" = korelasitest)
kort

##                     Nilai.Korelasi
## IHSG vs Kurs Dollar     -0.3219294

Karena nilai korelasi antara nilai kurs USD dengan data testing IHSG sebesar -0.3219294 maka nilai kurs USD dengan data testing IHSG memiliki korelasi positif kuat. Hal ini mengindikasikan bahwa ketiga variabel tersebut memiliki hubungan yang dapat meningkatkan keakuratan peramalan dengan metode ARIMAX.

Pembentukan Model Regresi

reg1  <- lm(data.trans$data~kurs)
summary(reg1)

## 
## Call:
## lm(formula = data.trans$data ~ kurs)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.996e+11 -3.598e+10  1.187e+10  6.181e+10  1.384e+11 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 6.331e+10  7.504e+10   0.844    0.400  
## kurs        6.003e-02  2.412e-02   2.489    0.014 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.34e+10 on 141 degrees of freedom
## Multiple R-squared:  0.04209,    Adjusted R-squared:  0.0353 
## F-statistic: 6.196 on 1 and 141 DF,  p-value: 0.01397

Eksplorasi sisaan model regresi

#Eksplorasi 
sisaanx <- reg1$residuals 
par(mfrow=c(1,1)) 
qqnorm(sisaanx) 
qqline(sisaanx, col = "blue", lwd = 2) 

plot(c(1:length(sisaanx)),sisaanx) 

acf(sisaanx) 

pacf(sisaanx) 

par(mfrow = c(1,1))

Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan tidak menyebar normal ditandai dengan titik titik yang cenderung tidak mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung tidak sama menandakan bahwa sisaan memiliki ragam yang tidak homogen. Pada plot ACF juga menunjukan tails off slowly sehingga diidentifikasikan tidak stasioner dalam rataan.

Pengujian White Noise Regresi

H0 : ρ1 = ρ2 = …. = ρk = 0 untuk k = 1, 2, …, n (Sisaan saling bebas atau sisaan white noise)

H1 : minimal ada satu ρk ≠ 0 , untuk k = 1, 2, …, n (Sisaan tidak saling bebas atau sisaan tidak terjadi white noise)

Box.test(sisaanx, lag = 5)

## 
##  Box-Pierce test
## 
## data:  sisaanx
## X-squared = 564.52, df = 5, p-value < 2.2e-16

Box.test(sisaanx, lag = 10)

## 
##  Box-Pierce test
## 
## data:  sisaanx
## X-squared = 959.45, df = 10, p-value < 2.2e-16

Box.test(sisaanx, lag = 15)

## 
##  Box-Pierce test
## 
## data:  sisaanx
## X-squared = 1240.3, df = 15, p-value < 2.2e-16

Box.test(sisaanx, lag = 20)

## 
##  Box-Pierce test
## 
## data:  sisaanx
## X-squared = 1434.4, df = 20, p-value < 2.2e-16

Berdasarkan LJung-Box test, diperoleh p-value 2.2 x 10-16 < α = 0.05, maka tolak H0. Artinya, cukup bukti untuk menyatakan bahwa sisaan antara lag tidak saling bebas atau sisaan tidak terjadi white noise pada taraf nyata 5%.

Cek Stasioner Model Regresi

Plot ACF

acf(sisaanx, lag.max = 24, main = "Plot ACF Sisaan Regresi")

Berdasarkan Plot ACF Sisaan, nilai korelasi antar lag terlihat pada plot di atas menurun secara perlahan (tails off slowly). Hal tersebut mengindikasikan bahwa sisaan tidak stasioner. Perlu dilakukan uji formal untuk mengambil kesimpulan kestasioneran data

Plot PACF

pacf(sisaanx, lag.max = 24, main = "Plot PACF Sisaan Regresi")

Uji Formal

Secara formal, metode Augmented Dickey-Fuller (ADF) dapat memberikan hasil uji secara akurat untuk menentukan apakah sebuah data stasioner atau tidak. Namun, Uji ADF ini hanya mengukur tingkat stasioneritas berdasarkan nilai tengah saja. Dengan hipotesis yang diuji sebagai berikut : H0 : Nilai tengah sisaan tidak stasioner H1 : Nilai tengah sisaan stasioner α = 5% = 0.05

tseries::adf.test(sisaanx)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  sisaanx
## Dickey-Fuller = -1.9918, Lag order = 5, p-value = 0.5798
## alternative hypothesis: stationary

Berdasarkan hasil Augmented Dickey-Fuller Test (ADF Test) didapatkan p-value = 0.5517 > α, maka tak tolak H0. Artinya, tidak cukup bukti untuk mengatakan bahwa sisaan stasioner pada taraf nyata 5%. Sehingga, perlu dilakukan differencing sebelum melakukan penentuan model tentatif

Penanganan Ketidakstasioneran Data

Differencing 1

sisaan.dif <- diff(sisaanx, difference = 1)

Cek kestasioneran data

Uji Formal ADF

Pengujian menggunakan Augmented Dickey-Fuller Test H0: Nilai tengah sisaan tidak stasioner H1: Nilai tengah sisaan stasioner α= 5% = 0.05

tseries::adf.test(sisaan.dif)

## Warning in tseries::adf.test(sisaan.dif): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  sisaan.dif
## Dickey-Fuller = -5.8957, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan hasil Augmented Dickey-Fuller Test (ADF Test) didapatkan p-value = 0.01 < α, maka tolak H_0. Artinya, cukup bukti untuk mengatakan bahwa sisaan stasioner pada taraf nyata 5% setelah dilakukan differencing sebanyak 1 kali.

Identifikasi Model Sisaan Model Regresi

Plot ACF

acf(sisaan.dif, lag.max = 24, main = "Plot ACF Sisaan Setelah Differencing satu kali")

Berdasarkan plot ACF di atas,terlihat bahwa nilai korelasi antara data dengan lag seperti gambar di atas tidak turun secara perlahan, dimana pada plot ACF diperoleh cuts off pada lag ke-23 ## Plot PACF

pacf(sisaan.dif, lag.max = 24, main = "Plot PACF Sisaan Setelah Differencing satu kali")

Pada plot PACF diperoleh cut off di lag 21. Berdasarkan hasil eksplorasi di atas, model yang dapat dibentuk secara berurutan adalah ARIMAX(0,1,23) dan ARIMAX(23,1,0).

Plot EACF

eacf(sisaan.dif)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x o o o o o o o o o o  o  o  o 
## 4 x o x x o o o o o o o  o  o  o 
## 5 x o x o x o o o o o o  o  o  o 
## 6 o x x x o o o o o o o  o  o  o 
## 7 x x x x o o o o o o o  o  o  o

Dari Matriks EACF dapat diduga model yang cocok adalah model : 1. ARIMAX(0,1,1) 2. ARIMAX(1,1,1) 3. ARIMAX(0,1,2) 4. ARIMAX(2,1,2) 5. ARIMAX(1,1,2)

Identifikasi Model Tentatif

Berdasarkan plot ACF, PACF, dan matriks EACF, diperoleh 5 model tentatif beserta orde parameternya, sebagai berikut: 1. ARIMAX(0,1,1) 2. ARIMAX(1,1,1) 3. ARIMAX(0,1,2) 4. ARIMAX(2,1,2) 5. ARIMAX(1,1,2) 6. ARIMAX(2,1,1) 7. ARIMAX(3,1,1) 8. ARIMAX(3,1,2)

modelx1 <- Arima(train.trans, order = c(0,1,1), xreg = kurs.train, method = "ML")
modelx2 <- Arima(train.trans, order = c(1,1,1), xreg = kurs.train, method = "ML")
modelx3 <- Arima(train.trans, order = c(0,1,2), xreg = kurs.train, method = "ML")
modelx4 <- Arima(train.trans, order = c(2,1,2), xreg = kurs.train, method = "ML")
modelx5 <- Arima(train.trans, order = c(1,1,2), xreg = kurs.train, method = "ML")
modelx6 <- Arima(train.trans, order = c(2,1,1), xreg = kurs.train, method = "ML")
modelx7 <- Arima(train.trans, order = c(3,1,1), xreg = kurs.train, method = "ML")
modelx8 <- Arima(train.trans, order = c(3,1,2), xreg = kurs.train, method = "ML")

ARIMAX(0,1,1)

coeftest(modelx1)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)   
## ma1   0.093083   0.100138  0.9295 0.352606   
## xreg -0.044619   0.014401 -3.0984 0.001946 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(modelx1$aic, 3)

## [1] 6293.844

ARIMAX (1,1,1)

coeftest(modelx2)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)   
## ar1  -0.277892   0.475334 -0.5846 0.558799   
## ma1   0.374376   0.456421  0.8202 0.412078   
## xreg -0.043615   0.014401 -3.0287 0.002456 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(modelx2$aic, 3)

## [1] 6295.467

ARIMAX(0,1,2)

coeftest(modelx3)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)   
## ma1   0.059302   0.091617  0.6473 0.517449   
## ma2  -0.111329   0.090874 -1.2251 0.220542   
## xreg -0.041726   0.014153 -2.9481 0.003197 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(modelx3$aic, 3)

## [1] 6294.363

ARIMAX(2,1,2)

coeftest(modelx4)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  -0.066634   0.026213  -2.5421  0.011020 *  
## ar2  -0.973320   0.021676 -44.9040 < 2.2e-16 ***
## ma1   0.142927   0.022329   6.4011 1.543e-10 ***
## ma2   1.000000   0.045076  22.1850 < 2.2e-16 ***
## xreg -0.042431   0.013803  -3.0741  0.002111 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(modelx4$aic, 3)

## [1] 6291.003

ARIMAX(1,1,2)

coeftest(modelx5)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)   
## ar1   0.430940   0.381872  1.1285 0.259112   
## ma1  -0.363207   0.374908 -0.9688 0.332651   
## ma2  -0.153220   0.083622 -1.8323 0.066906 . 
## xreg -0.040790   0.013938 -2.9264 0.003429 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(modelx5$aic, 3)

## [1] 6295.48

ARIMAX(3,1,0)

coeftest(modelx6)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)   
## ar1   0.459247   0.434954  1.0559 0.291036   
## ar2  -0.172772   0.090095 -1.9177 0.055153 . 
## ma1  -0.382338   0.439440 -0.8701 0.384270   
## xreg -0.040041   0.013860 -2.8890 0.003865 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(modelx6$aic, 3)

## [1] 6295.069

ARIMAX(3,1,1)

coeftest(modelx7)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)   
## ar1  -0.205143   0.551459 -0.3720 0.709892   
## ar2  -0.105144   0.101139 -1.0396 0.298524   
## ar3  -0.133939   0.105545 -1.2690 0.204433   
## ma1   0.280516   0.551635  0.5085 0.611091   
## xreg -0.040963   0.013840 -2.9597 0.003079 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(modelx7$aic, 3)

## [1] 6296.475

ARIMAX(3,1,2)

coeftest(modelx8)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  -0.053927   0.092971  -0.5800  0.561887    
## ar2  -0.972845   0.021878 -44.4676 < 2.2e-16 ***
## ar3   0.013085   0.091849   0.1425  0.886712    
## ma1   0.142399   0.022886   6.2222 4.902e-10 ***
## ma2   0.999999   0.044814  22.3145 < 2.2e-16 ***
## xreg -0.042273   0.013823  -3.0582  0.002227 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(modelx8$aic, 3)

## [1] 6292.983

Hasil uji signifikansi hanya ARIMAX (2,1,2) karena terdapat tiga parameternya lebih kecil dari taraf nyata 5%. Selain itu ARIMAX (2,1,2) memiliki AIC yang tergolong kecil dari pada model lainnya yaitu 6291.003. Selanjutnya model ARIMAX (2,1,2) akan di-overfitting.

Diagnostik Model ARIMAX(2,1,2)

Eksplorasi Sisaan

#Eksplorasi 
sisaan.max <- modelx4$residuals 
par(mfrow=c(1,1)) 
qqnorm(sisaan.max) 
qqline(sisaan.max, col = "blue", lwd = 2) 

plot(c(1:length(sisaan.max)),sisaan.max) 

acf(sisaan.max) 

pacf(sisaan.max) 

par(mfrow = c(1,1))

Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan tidak menyebar normal ditandai dengan titik titik yang cenderung tidak mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki ragam yang homogen. Plot ACF dan PACF sisaan ARIMA(2,1,2) juga tidak signifikan pada 20 lag awal yang menandakan sisaan saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.

Uji Formal

1) Sisaan Menyebar Normal

ks.test(sisaan.max,"pnorm")  #tak tolak H0 > sisaan menyebar normal

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.max
## D = 0.58594, p-value < 2.2e-16
## alternative hypothesis: two-sided

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS). Hipotesis pada uji KS adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value sebesar 2.2e-16 yang kurang dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Hal ini sesuai dengan hasil eksplorasi menggunakan plot kuantil-kuantil normal.

2) Sisaan saling bebas/tidak ada autokorelasi

Box.test(sisaan.max, type = "Ljung")  #tak tolak H0 > sisaan saling bebas

## 
##  Box-Ljung test
## 
## data:  sisaan.max
## X-squared = 0.0098505, df = 1, p-value = 0.9209

Selanjutnya akan dilakukan uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box. Hipotesis yang digunakan adalah sebagai berikut.

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar 0.9209 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas.

3) Sisaan homogen

Box.test((sisaan.max)^2, type = "Ljung")  #tak tolak H0 > sisaan homogen

## 
##  Box-Ljung test
## 
## data:  (sisaan.max)^2
## X-squared = 2.4357, df = 1, p-value = 0.1186

Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.1186 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa ragam sisaan homogen.

4) Nilai tengah sisaan sama dengan nol

t.test(sisaan.max, mu = 0, conf.level = 0.95)  #tak tolak h0 > nilai tengah sisaan sama dengan 0

## 
##  One Sample t-test
## 
## data:  sisaan.max
## t = 1.553, df = 127, p-value = 0.1229
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -482944359 4005962385
## sample estimates:
##  mean of x 
## 1761509013

Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol. Hipotesis yang diujikan sebagai berikut.

\(H_0\) : nilai tengah sisaan sama dengan 0

\(H_1\) : nilai tengah sisaan tidak sama dengan 0

Berdasarkan uji-ttersebut, didapat p-value sebesar 0.1229 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol.

Pengujian White Noise

H0: ρ1= ρ2= …. = ρk= 0 untuk k = 1, 2, …, n (Sisaan saling bebas atau sisaan white noise)

H1: minimal ada satu ρk ≠ 0 , untuk k = 1, 2, …, n (Sisaan tidak saling bebas atau sisaan tidak terjadi white noise)

Box.test(sisaan.max, lag = 5)

## 
##  Box-Pierce test
## 
## data:  sisaan.max
## X-squared = 2.1906, df = 5, p-value = 0.8222

Box.test(sisaan.max, lag = 5)

## 
##  Box-Pierce test
## 
## data:  sisaan.max
## X-squared = 2.1906, df = 5, p-value = 0.8222

Box.test(sisaan.max, lag = 10)

## 
##  Box-Pierce test
## 
## data:  sisaan.max
## X-squared = 4.5286, df = 10, p-value = 0.9204

Box.test(sisaan.max, lag = 15)

## 
##  Box-Pierce test
## 
## data:  sisaan.max
## X-squared = 10.032, df = 15, p-value = 0.8177

Berdasarkan LJung-Box test, untuk setiap lag diperoleh p-value > α = 0.05, maka terima H0. Artinya, tidak cukup bukti untuk menyatakan bahwa sisaan antara lag tidak saling bebas atau sisaan tidak terjadi white noise pada taraf nyata 5%.

Overfitting

ARIMAX(3,1,2)

model1x.over <- Arima(train.trans, order = c(3,1,2), xreg = kurs.train, method = "ML")
coeftest(model1x.over)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  -0.053927   0.092971  -0.5800  0.561887    
## ar2  -0.972845   0.021878 -44.4676 < 2.2e-16 ***
## ar3   0.013085   0.091849   0.1425  0.886712    
## ma1   0.142399   0.022886   6.2222 4.902e-10 ***
## ma2   0.999999   0.044814  22.3145 < 2.2e-16 ***
## xreg -0.042273   0.013823  -3.0582  0.002227 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(model1x.over$aic, 3)

## [1] 6292.983

ARIMAX(2,1,3)

model2x.over <- Arima(train.trans, order = c(2,1,3), xreg = kurs.train, method = "ML")
coeftest(model2x.over)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  -0.067477   0.026724  -2.5250  0.011571 *  
## ar2  -0.973822   0.021790 -44.6907 < 2.2e-16 ***
## ma1   0.157309   0.098377   1.5990  0.109812    
## ma2   1.002047   0.046939  21.3480 < 2.2e-16 ***
## ma3   0.014996   0.099568   0.1506  0.880280    
## xreg -0.042251   0.013832  -3.0545  0.002255 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(model2x.over$aic, 3)

## [1] 6292.98

Karena parameter ar1 dan ar2 pada model ARIMAX(3,1,2) tidak berubah secara signifikan dari parameter ar1 dan parameter ar2 pada model ARIMAX(2,1,2), Pada parameter ma1 dan ma2 pada model ARIMAX(3,1,2) tidak berubah secara signifikan dari parameter ma1 dan parameter ma2 pada model ARIMAX(2,1,2). Selain itu diperoleh bahwa model ARIMAX(2,1,3) dan ARIMAX(3,1,2) tidak signifikan pada beberapa parameter dan nilai AIC lebih besar daripada ARIMAX(2,1,2).

Maka model ARIMAX(2,1,2) tidak terjadi overfitting atau merupakan model yang baik untuk digunakan. (Cryer dan Chan. 2008)

Prediksi ARIMAX

#ramalanx <- forecast::forecast((modelx4), xreg = cbind(kurs.test))
#ram<- ramalanx$mean
ramalan <- predict(modelx4,  n.ahead = 15, newxreg = data.frame(kurs.test))

hasil_arimax<- ramalan$pred
Hasil2 <- (hasil_arimax)^(1/3)
ts.plot(data.ts,xlab = "Periode", ylab = "Harga Saham IHSG", col="black",lwd=2,main="Forecasting ARIMAX(2,1,2)",gpars = list(col.main="black",col.axis="black",col.sub="black"))
lines(Hasil2, col = "blue",lwd=2)
lines(test.ts, col = "red", lwd =2)
legend("bottomright", 100,20,legend = c("Data Training", "Data Testing", "Data Forecast ARIMAX(2,1,2)"), 
       lwd=2, col=c("black","red","blue"), cex=0.8)
box(col="black",lwd=2)

Akurasi ARIMAX

perbandingan <- data.frame(Aktual=c(test.ts),
                           Predik = Hasil2)
perbandingan

##      Aktual   Predik
## 1  7178.583 7120.776
## 2  7040.798 7040.214
## 3  7026.783 7070.358
## 4  6814.530 7075.787
## 5  7017.771 6996.366
## 6  7056.040 6975.944
## 7  7045.527 6975.269
## 8  7089.206 7060.545
## 9  7082.181 6979.809
## 10 7053.150 6952.320
## 11 7019.639 7044.771
## 12 6715.118 7039.026
## 13 6812.193 6994.682
## 14 6800.673 6972.269
## 15 6835.808 7019.839

A <- accuracy(ts(Hasil2), head(data.test, n=length(data.test)))
print(A)

##                 ME     RMSE      MAE       MPE     MAPE
## Test set -48.66502 143.5123 110.2667 -0.733301 1.604003

par(mfrow=c(1,2))
ts.plot(data.ts,xlab = "Periode", ylab = "Harga Saham IHSG", col="black",lwd=2,main="Forecasting ARIMA(2,1,2)",ylim = c(4000, max(Hasil1)),gpars = list(col.main="black",col.axis="black",col.sub="black"))
lines(Hasil1, col = "blue",lwd=2)
lines(test.ts, col = "red", lwd =2)
legend("bottomright", 20,20,legend = c("Data Training", "Data Testing", "Data Forecast ARIMA(2,1,0)"), 
       lwd=2, col=c("black","red","blue"), cex=0.5)
box(col="black",lwd=2)

ts.plot(data.ts,xlab = "Periode", ylab = "Harga Saham IHSG", col="black",lwd=2,main="Forecasting ARIMAX(2,1,2)",ylim = c(4000, max(Hasil1)),gpars = list(col.main="black",col.axis="black",col.sub="black"))
lines(Hasil2, col = "blue",lwd=2)
lines(test.ts, col = "red", lwd =2)
legend("bottomright", 20,20,legend = c("Data Training", "Data Testing", "Data Forecast ARIMAX(2,1,2)"), 
       lwd=2, col=c("black","red","blue"), cex=0.5)
box(col="black",lwd=2)

Perbandingan ARIMA dan ARIMAX

perbandingan_arima_arimax <- data.frame(Aktual=c(test.ts),
                           Predik_ARIMAX = Hasil2,
                           Predik_ARIMA = Hasil1)
perbandingan

##      Aktual   Predik
## 1  7178.583 7120.776
## 2  7040.798 7040.214
## 3  7026.783 7070.358
## 4  6814.530 7075.787
## 5  7017.771 6996.366
## 6  7056.040 6975.944
## 7  7045.527 6975.269
## 8  7089.206 7060.545
## 9  7082.181 6979.809
## 10 7053.150 6952.320
## 11 7019.639 7044.771
## 12 6715.118 7039.026
## 13 6812.193 6994.682
## 14 6800.673 6972.269
## 15 6835.808 7019.839

akurasi1<- data.frame(accuracy(ts(Hasil1), head(data.test, n=length(data.test))))
akurasi2<- data.frame(accuracy(ts(Hasil2), head(data.test, n=length(data.test))))
akurasi <- rbind(akurasi1,akurasi2)
rownames(akurasi)<- c("ARIMA","ARIMAX")
print(akurasi)

##                ME     RMSE      MAE       MPE     MAPE
## ARIMA  -147.29945 204.8412 152.8213 -2.151847 2.228769
## ARIMAX  -48.66502 143.5123 110.2667 -0.733301 1.604003

Berdasarkan nilai akurasi kedua model di atas, terlihat MAPE ARIMAX(2,1,2) menghasilkan MAPE yang lebih kecil dibandingan ARIMA(2,1,2). Oleh karena itu, model ARIMAX(2,1,2) dianggap lebih baik dalam meramalkan IHSG dengan peubah kovariat nya kurs dolar.