1.1.1 Import Data

install_load('rio')
raw.data <- import("https://raw.githubusercontent.com/Zen-Rofiqy/STA1341-MPDW/main/Data/MAANG%20Stock%20Prices.csv")

1.1.2 Data Checking

Cek Tipe data.

str(raw.data)

## 'data.frame':    8812 obs. of  8 variables:
##  $ Name     : chr  "AMZN" "AMZN" "AMZN" "AMZN" ...
##  $ Date     : chr  "7/30/18" "7/31/18" "8/1/18" "8/2/18" ...
##  $ Open     : chr  "91.366501" "89.324501" "89.199997" "89.438499" ...
##  $ High     : chr  "91.474998" "90.091499" "89.921997" "91.828003" ...
##  $ Low      : chr  "88.301003" "86.966003" "88.801003" "89.300003" ...
##  $ Close    : chr  "88.960999" "88.872002" "89.858498" "91.716499" ...
##  $ Adj Close: chr  "88.960999" "88.872002" "89.858498" "91.716499" ...
##  $ Volume   : chr  "131246000" "114774000" "83062000" "87094000" ...

Semua data Karakter, harus diubah.

Cek Data kosong.

sum(is.na(raw.data))

## [1] 0

Tidak ada data kosong.

1.1.3 Penyesuaian Tipe Data

Semua tipe data masih berupa character. Harus diubah menjadi tipe data yang sesuai.

install_load('dplyr')
data <- raw.data %>%  
  mutate(
    Date = as.Date(raw.data[, 2], format = "%m/%d/%y"), #Mengubah menjadi Date 
    across(3:ncol(raw.data), as.numeric)                #Mengubah menjadi Numerik
  )
str(data)

## 'data.frame':    8812 obs. of  8 variables:
##  $ Name     : chr  "AMZN" "AMZN" "AMZN" "AMZN" ...
##  $ Date     : Date, format: "2018-07-30" "2018-07-31" ...
##  $ Open     : num  91.4 89.3 89.2 89.4 91.9 ...
##  $ High     : num  91.5 90.1 89.9 91.8 92.1 ...
##  $ Low      : num  88.3 87 88.8 89.3 91.1 ...
##  $ Close    : num  89 88.9 89.9 91.7 91.2 ...
##  $ Adj Close: num  89 88.9 89.9 91.7 91.2 ...
##  $ Volume   : num  1.31e+08 1.15e+08 8.31e+07 8.71e+07 6.92e+07 ...

1.1.4 Rechecking Data 

Cek kembali data kosong.

cat('Banyaknya Data Kosong', sum(is.na(data)))

## Banyaknya Data Kosong 42

Melihat baris, kolom mana data yang kosong.

# Mencari indeks baris dan kolom yang mengandung NA
na.idx <- which(is.na(data), arr.ind = TRUE)

# Menampilkan data raw dengan baris dan kolom yang mengandung NA
install_load('DT')

## Warning: package 'DT' was built under R version 4.2.3

datatable(raw.data[                        # Subsetting
                unique(na.idx[, 1]),   # Vektor indeks baris yang mengandung NA
                unique(na.idx[, 2])  ] # Vektor indeks kolom yang mengandung NA
          )  

Ternyata pada baris tersebut ada data karakter text yang merupakan label dari tiap kolomnya. Sehingga ketika diubah ke numerik akan menjadi NA. Maka saya akan menghapus baris tersebut.

1.1.5 Data Cleaned

data <- data %>%
  filter(!row_number() %in% unique(na.idx[, 1]))
datatable(data, filter = 'top', 
          options = list(pageLength = 5))

1.2.1 Time Series MAANG

Melihat keseluruhan Time Series data saham.

install_load('viridis','ggrepel')
#Plot
chart <-
ggplot(data, aes(x=Date, y=`Adj Close`, color=Name, alpha=Name)) + #Data
  geom_line(aes(color=Name), linewidth=1.5) + #Timeseries
  #Color
  scale_color_manual(values = c(NVDA="green4", NFLX="red4", MSFT="steelblue3",
                                META="royalblue4", AAPL="lightskyblue4",
                                GOOG="yellow3", AMZN="orange2") ) +
  scale_alpha_manual(values = c("NVDA" = .25, "NFLX" = .25, "MSFT" = .25, 
                                "META" = .25, "AAPL" = .25, "GOOG" = .25, 
                                "AMZN" = 1)) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Time Series MAANG",
       subtitle = "Seperti apa sih pola deret waktu saham MAANG?\n") +
  # Label / legend
  geom_text_repel(
    data=data[data$Date == max(data$Date),], #Posisi di ujung data
    aes(color = Name, label = Name), #Warna garis & label saham
    size = 8, #Ukuran text
    nudge_x = 80, #Posisi Text (kanan 50)
    hjust = 0, #Ujung
    segment.size = 1,               #Ukuran garis
    segment.alpha = .75,             #transparasi garis
    segment.linetype = "dotted",    #Time garis
    box.padding = .4, #Biar label saham nggak dempetan
    segment.curvature = -0.1, #biar garis mulus
    segment.ncp = 8, 
    segment.angle = 60 
  ) +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date("2018-07-30"), as.Date("2023-12-28"))
    #Tampilin lebih dari 20023-07-28 agar label saham bisa masuk
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
    annotate( #Buat nandain batas data
    "text", x = as.Date("2023-7-28"), y = 50, 
    label = "28 Juli", size=6
  ) +
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date("2023-07-28")), 
             linetype = "dotted", color = "red")
chart

#Export Chart
ggsave("01_Time Series MAANG.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 23)

Jika dilihat dari tahun 2019-2022, semua saham cenderung memiliki pola trend naik. Lalu dari 2021-2023 polanya cenderung trend turun. Untuk tugas praktikum kali ini, saya hanya akan menggunakan rentang tahun 2022-2023 dengan tren cenderung turun. Agar pengerjaannya tidak terlalu sulit, karena masih tahap awal pembelajaran.

Ada yang menarik perhatian saya. Kenapa dulu Pendiri Amazon, Jeff Bezos yang pernah menjadi orang terkaya di dunia pada tahun 2017 lalu, harga saham sekarang tidak setinggi yang saya kira?. Oleh karena itu saya memutuskan untuk menggunakan data saham Amazon untuk praktikum kali ini.

amzn <- data %>%
  select(1, 2, 7) %>%  # Memilih kolom 1, 2, dan 7
  filter(Name == "AMZN", Date >= as.Date("2022-01-01"))  # Filter data saham Amazon tahun 2022 ke atas

rownames(amzn) <- NULL
str(amzn)

## 'data.frame':    394 obs. of  3 variables:
##  $ Name     : chr  "AMZN" "AMZN" "AMZN" "AMZN" ...
##  $ Date     : Date, format: "2022-01-03" "2022-01-04" ...
##  $ Adj Close: num  170 168 164 163 163 ...

datatable(amzn)

Mengubah Ajd Close Menjadi Time series.

amzn.ts <- ts(amzn[,3])

Ringkasan Data Ajd CLose.

summary(amzn.ts)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   81.82  102.01  115.77  119.81  134.88  170.40

1.2.2 Time Series Amazon

min_value <- min(amzn$`Adj Close`)
min_date <- amzn$Date[which.min(amzn$`Adj Close`)]
percentage <- (which.min(amzn$`Adj Close`) / nrow(amzn)) * 100

chart <-
ggplot(amzn, aes(x=Date, y=`Adj Close`)) + 
  geom_line(aes(color=Name), linewidth=2) +
  scale_color_manual(values = c("orange2")) +
  labs(x = "\nPeriode (Tahun)", y='Saham Harga penutup',
       title = "Time Series Saham Amazon",
       subtitle = "Seperti apa sih pola deret waktu saham Amazon?\n") +
  theme(legend.position = "none") +
  theme.ts1 + 
  geom_vline(xintercept = as.numeric(min_date), 
             linetype = "dotted", color = "cyan4", linewidth = 1.5) +
  geom_text(aes(x = min_date-1*40, y = max(`Adj Close`)/1.59, label = 
                  paste0("Titik Potong\n","(",round(percentage, 2), "%)",
                         "   ",min_date)), 
            vjust = -1.5, hjust = 0, size = 7, color = "cyan4") 
chart

#Export Chart
ggsave("02_Time Series Amazon.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 20)

Berdasarkan chart data deret waktu diatas, terlihat bahwa data cenderung memiliki trend naik dan turun dengan titik potong di tanggal 20222-12-28. Namun titik potongnya masih jauh dari angka \(80\%\). Sehingga tidak akan dijadikan acuan untuk pemisahan data training dan testing.

Pembagian Data Training Dan Test.

#membagi 80% data latih (training) dan 20% data uji (testing)
training <- amzn[1: round(nrow(amzn) *80/100),]
testing <- amzn[round(nrow(amzn) *80/100): nrow(amzn),]
train.ts <- ts(training[,3])
test.ts <- ts(testing[,3])

chart <-
ggplot() + 
  geom_line(data = training, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "Data Latih")) +
  geom_line(data = testing, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "Data Uji")) +
  labs(x = "\nPeriode (Tahun)", y='Saham Harga penutup',
       title = "Time Series Saham Amazon",
       subtitle = "Pembagian Data Training dan Test\n") +
  theme(legend.position = "none") +
  scale_colour_manual(name="Keterangan:", 
                      breaks = c("Data Latih", "Data Uji"),
                      values = c("orange", "cyan4")) + theme.ts1
chart

#Export Chart
ggsave("02_TSA_train-test.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 20)

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

library(ggplot2)
library(tsibble)
library(tseries)

5.2.1 ARIMA(0,1,1)

model1.da=Arima(train.diff, order=c(0,1,1),method="ML")
summary(model1.da)

## Series: train.diff 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -1.0000
## s.e.   0.0131
## 
## sigma^2 = 13.87:  log likelihood = -858.02
## AIC=1720.04   AICc=1720.08   BIC=1727.53
## 
## Training set error measures:
##                     ME     RMSE      MAE     MPE     MAPE      MASE        ACF1
## Training set 0.2263963 3.711957 2.701246 102.999 130.2316 0.7194769 -0.01544094

lmtest::coeftest(model1.da)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.999994   0.013122 -76.205 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar \(1720.04\) dan seluruh parameter signifikan.

Sesuai dengan function

5.2.2 ARIMA(1,1,1)

model2.da=Arima(train.diff, order=c(1,1,1),method="ML")
summary(model2.da)

## Series: train.diff 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.0099  -1.0000
## s.e.   0.0566   0.0133
## 
## sigma^2 = 13.91:  log likelihood = -858
## AIC=1722.01   AICc=1722.08   BIC=1733.25
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.2282761 3.711646 2.704446 102.3478 130.0943 0.7203291
##                      ACF1
## Training set -0.004925212

lmtest::coeftest(model2.da)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value Pr(>|z|)    
## ar1 -0.0099178  0.0565844  -0.1753   0.8609    
## ma1 -0.9999999  0.0133448 -74.9355   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar \(1722.01\) dan parameter ar1 tidak signifikan.

Sesuai dengan function

5.2.3 ARIMA(2,1,2)

model3.da=Arima(train.diff, order=c(2,1,2),method="ML")
summary(model3.da)

## Series: train.diff 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1     ar2      ma1      ma2
##       -0.3613  0.0603  -0.6457  -0.3542
## s.e.   1.1394  0.0650   1.1429   1.1427
## 
## sigma^2 = 13.94:  log likelihood = -857.28
## AIC=1724.55   AICc=1724.75   BIC=1743.28
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.2203966 3.703627 2.703416 95.80655 127.7273 0.7200548
##                      ACF1
## Training set -0.005979234

lmtest::coeftest(model3.da)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.361269   1.139416 -0.3171   0.7512
## ar2  0.060255   0.064954  0.9277   0.3536
## ma1 -0.645719   1.142861 -0.5650   0.5721
## ma2 -0.354238   1.142704 -0.3100   0.7566

Didapatkan nilai AIC sebesar \(1724.55\) dan semua parameter tidak signifikan.

Sesuai dengan function

5.2.4 ARIMA(3,1,3)

model4.da=Arima(train.diff, order=c(3,1,3),method="ML")
summary(model4.da)

## Series: train.diff 
## ARIMA(3,1,3) 
## 
## Coefficients:
##           ar1     ar2     ar3      ma1      ma2     ma3
##       -0.1041  0.3262  0.0017  -0.9038  -0.3547  0.2585
## s.e.   0.9780  0.4728  0.0887   0.9764   0.7967  0.4738
## 
## sigma^2 = 14.02:  log likelihood = -857.17
## AIC=1728.33   AICc=1728.7   BIC=1754.56
## 
## Training set error measures:
##                     ME    RMSE      MAE      MPE     MAPE      MASE
## Training set 0.2154983 3.70268 2.705672 96.69477 128.6789 0.7206558
##                      ACF1
## Training set -0.005229142

lmtest::coeftest(model4.da)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)
## ar1 -0.1040802  0.9780405 -0.1064   0.9153
## ar2  0.3261986  0.4727589  0.6900   0.4902
## ar3  0.0017341  0.0886766  0.0196   0.9844
## ma1 -0.9038021  0.9763744 -0.9257   0.3546
## ma2 -0.3546642  0.7967347 -0.4451   0.6562
## ma3  0.2584697  0.4737918  0.5455   0.5854

Didapatkan nilai AIC sebesar \(1728.33\) dan semua parameter tidak signifikan.

Sesuai dengan function

5.2.5 ARIMA(4,1,5)

model5.da=Arima(train.diff, order=c(4,1,5),method="ML")
summary(model5.da)

## Series: train.diff 
## ARIMA(4,1,5) 
## 
## Coefficients:
##           ar1     ar2      ar3      ar4      ma1      ma2     ma3     ma4
##       -0.2750  0.4927  -0.2505  -0.8488  -0.6975  -0.7685  0.7276  0.7147
## s.e.   0.0519  0.0554   0.0493   0.0482   0.0348   0.0276  0.0424  0.0297
##           ma5
##       -0.9762
## s.e.   0.0317
## 
## sigma^2 = 13.22:  log likelihood = -849.02
## AIC=1718.04   AICc=1718.77   BIC=1755.5
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set 0.207726 3.577948 2.622015 107.5843 161.1924 0.6983738 -0.01911452

lmtest::coeftest(model5.da)

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.275029   0.051891  -5.3001 1.157e-07 ***
## ar2  0.492659   0.055361   8.8991 < 2.2e-16 ***
## ar3 -0.250506   0.049259  -5.0855 3.667e-07 ***
## ar4 -0.848818   0.048201 -17.6100 < 2.2e-16 ***
## ma1 -0.697516   0.034758 -20.0675 < 2.2e-16 ***
## ma2 -0.768509   0.027565 -27.8799 < 2.2e-16 ***
## ma3  0.727565   0.042374  17.1699 < 2.2e-16 ***
## ma4  0.714674   0.029735  24.0347 < 2.2e-16 ***
## ma5 -0.976213   0.031736 -30.7602 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

best.model.da <- model5.da

6.2.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.27186, 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.

6.2.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 = 0.11582, df = 1, p-value = 0.7336

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.7336\) yang lebih besar dari taraf nyata \(5\%\) sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini sesuai dengan hasil eksplorasi menggunakan plot ACF dan PACF.

6.2.3 Sisaan homogen

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

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

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.01499\) yang lebih kecil dari taraf nyata \(5\%\) sehingga tolak \(H_0\) dan menandakan bahwa ragam sisaan tidak 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.0289, df = 313, p-value = 0.3043
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1895198  0.6049718
## sample estimates:
## mean of x 
##  0.207726

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.3043\) yang lebih besar dari taraf nyata \(5\%\) sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol.