R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

install.packages("rlang", type = "binary")

Error in install.packages : Updating loaded packages

install.packages(c("readr", "dplyr", "lubridate", "forecast", "tseries", "ggplot2"))

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing packages into ‘C:/Users/alice/AppData/Local/R/win-library/4.4’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/readr_2.2.0.zip'
Content type 'application/zip' length 1194020 bytes (1.1 MB)
downloaded 1.1 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/dplyr_1.2.1.zip'
Content type 'application/zip' length 1618178 bytes (1.5 MB)
downloaded 1.5 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/lubridate_1.9.5.zip'
Content type 'application/zip' length 991003 bytes (967 KB)
downloaded 967 KB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/forecast_9.0.2.zip'
Content type 'application/zip' length 1965798 bytes (1.9 MB)
downloaded 1.9 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/tseries_0.10-61.zip'
Content type 'application/zip' length 402806 bytes (393 KB)
downloaded 393 KB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/ggplot2_4.0.3.zip'
Content type 'application/zip' length 8464184 bytes (8.1 MB)
downloaded 8.1 MB

package ‘readr’ successfully unpacked and MD5 sums checked
package ‘dplyr’ successfully unpacked and MD5 sums checked
package ‘lubridate’ successfully unpacked and MD5 sums checked
package ‘forecast’ successfully unpacked and MD5 sums checked
package ‘tseries’ successfully unpacked and MD5 sums checked
package ‘ggplot2’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\alice\AppData\Local\Temp\RtmpYJ1XGg\downloaded_packages

library(readr)

Warning: package ‘readr’ was built under R version 4.4.3

library(dplyr)

Warning: package ‘dplyr’ was built under R version 4.4.3
Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(lubridate)

Warning: package ‘lubridate’ was built under R version 4.4.3
Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

library(forecast)

Warning: package ‘forecast’ was built under R version 4.4.3

library(tseries)

Warning: package ‘tseries’ was built under R version 4.4.3Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

    ‘tseries’ version: 0.10-61

    ‘tseries’ is a package for time series analysis and computational finance.

    See ‘library(help="tseries")’ for details.

library(ggplot2)

Warning: package ‘ggplot2’ was built under R version 4.4.3

data <- read.csv("C:/Users/alice/Downloads/Employment on manufacturing sector.csv")

str(data)

'data.frame':   1848 obs. of  4 variables:
 $ series    : chr  "abs" "abs" "abs" "abs" ...
 $ sector    : chr  "p0" "p1" "p2" "p3" ...
 $ date      : chr  "1/01/2015" "1/01/2015" "1/01/2015" "1/01/2015" ...
 $ employment: num  13947 1857 78 2382 391 ...

data_p0 <- data %>%
  filter(series == "abs", sector == "p0") %>%
  mutate(date = dmy(date)) %>%
  arrange(date)

View(data_p0)
nrow(data_p0)

[1] 44

employment_ts <- ts(
  data_p0$employment,
  start = c(2015, 1),
  frequency = 4
)

employment_ts

      Qtr1  Qtr2  Qtr3  Qtr4
2015 13947 13944 13950 14128
2016 14096 14162 14202 14261
2017 14356 14415 14484 14581
2018 14684 14767 14857 14933
2019 15010 15078 15162 15255
2020 15243 14884 15096 15162
2021 15237 15207 15275 15441
2022 15575 15701 15831 15942
2023 16062 16146 16251 16347
2024 16401 16591 16715 16793
2025 16703 16849 16967 17100

frequency(employment_ts)

[1] 4

start(employment_ts)

[1] 2015    1

end(employment_ts)

[1] 2025    4

plot(data_p0$date,
     data_p0$employment,
     type = "o",
     col = "blue",
     pch = 16,
     cex = 1.2,
     lwd = 2.5,
     xaxt = "n",
     cex.main = 1.5,
     cex.lab = 1.3,
     cex.axis = 1.1,
     main = "Quarterly Employment in Malaysia's Manufacturing Sector (2015 Q1–2025 Q4)",
     xlab = "Year",
     ylab = "Employment")

axis.Date(1,
          at = seq(min(data_p0$date),
                   max(data_p0$date),
                   by = "year"),
          format = "%Y",
          cex.axis = 1.1)

table51 <- data.frame(
  Parameter = c("Alpha (α)",
                "Beta (β)",
                "Initial Level (l0)",
                "Initial Trend (b0)"),
  Estimate = c(0.9137,
               0.0001,
               13845.7654,
               73.7071)
)

table51

# Pastikan modelnya didefinisikan dulu di dalam chunk ini
Holt_est <- holt(employment_ts, h = 1)

acc_holt <- accuracy(Holt_est)

MSE_holt <- acc_holt[,"RMSE"]^2
MAPE_holt <- acc_holt[,"MAPE"]

table52 <- data.frame(
  Measure = c("MSE", "MAPE"),
  Value = c(MSE_holt, MAPE_holt)
)

table52

forecast_holt <- forecast(Holt_est, h = 1)

forecast_holt

autoplot(forecast_holt) +
  ggtitle("Holt's Linear Trend One-Step-Ahead Forecast") +
  xlab("Year") +
  ylab("Employment")

table53 <- data.frame(
  Forecast_Period = "2026 Q1",
  Forecast_Value = as.numeric(forecast_holt$mean)
)

table53

model_arima120 <- Arima(employment_ts, order = c(1,2,0))

summary(model_arima120)

Series: employment_ts 
ARIMA(1,2,0) 

Coefficients:
          ar1
      -0.4972
s.e.   0.1305

sigma^2 = 12903:  log likelihood = -257.99
AIC=519.97   AICc=520.28   BIC=523.45

Training set error measures:
                   ME    RMSE      MAE        MPE      MAPE      MASE       ACF1
Training set 4.146971 109.652 62.03111 0.02674741 0.4047804 0.2006181 -0.1713952

acc_arima120 <- accuracy(model_arima120)

MSE_arima120 <- acc_arima120[,"RMSE"]^2
MAPE_arima120 <- acc_arima120[,"MAPE"]

model_arima021 <- Arima(employment_ts, order = c(0,2,1))

summary(model_arima021)

Series: employment_ts 
ARIMA(0,2,1) 

Coefficients:
          ma1
      -0.9606
s.e.   0.1002

sigma^2 = 8503:  log likelihood = -250.34
AIC=504.68   AICc=504.99   BIC=508.16

Training set error measures:
                   ME     RMSE      MAE        MPE      MAPE      MASE        ACF1
Training set 14.88872 89.01417 54.69445 0.09416823 0.3561064 0.1768902 -0.09355206

acc_arima021 <- accuracy(model_arima021)

MSE_arima021 <- acc_arima021[,"RMSE"]^2
MAPE_arima021 <- acc_arima021[,"MAPE"]

table54 <- data.frame(
  Model = c("ARIMA(1,2,0)",
            "ARIMA(0,2,1)"),
  AIC = c(AIC(model_arima120),
          AIC(model_arima021)),
  BIC = c(BIC(model_arima120),
          BIC(model_arima021)),
  MSE = c(MSE_arima120,
          MSE_arima021),
  MAPE = c(MAPE_arima120,
           MAPE_arima021)
)

table54

forecast_arima021 <- forecast(model_arima021, h = 1)

forecast_arima021

autoplot(forecast_arima021) +
  ggtitle("ARIMA(0,2,1) One-Step-Ahead Forecast") +
  xlab("Year") +
  ylab("Employment")

table55 <- data.frame(
  Forecast_Period = "2026 Q1",
  Forecast_Value = as.numeric(forecast_arima021$mean)
)

table55

table56 <- data.frame(
  Model = c("Holt's Linear Trend",
            "ARIMA(1,2,0)",
            "ARIMA(0,2,1)"),
  MSE = c(MSE_holt,
          MSE_arima120,
          MSE_arima021),
  MAPE = c(MAPE_holt,
           MAPE_arima120,
           MAPE_arima021)
)

table56

model

Call:
lm(formula = electric ~ t, data = data)

Coefficients:
(Intercept)            t  
      97.84         0.47  

table57 <- data.frame(
  Best_Model = "Holt's Linear Trend",
  Forecast_Period = "2026 Q1",
  Forecast_Value = as.numeric(forecast_holt$mean)
)

table57

fit_part <- window(employment_ts, end = c(2022, 3))
holdout_part <- window(employment_ts, start = c(2022, 4))

length(fit_part)

[1] 31

length(holdout_part)

[1] 13

# Plot fitted and hold-out data
plot(employment_ts,
     main = "Fitted and Hold-out Data",
     xlab = "Year",
     ylab = "Employment",
     col = "black",
     lwd = 2)

lines(fit_part, col = "blue", lwd = 2)

lines(holdout_part, col = "red", lwd = 2)

legend("topleft",
       legend = c("Full Series", "Fitted Part", "Hold-out Part"),
       col = c("black", "blue", "red"),
       lwd = 2)

adf.test(fit_part)

    Augmented Dickey-Fuller Test

data:  fit_part
Dickey-Fuller = -1.7685, Lag order = 3, p-value = 0.6619
alternative hypothesis: stationary

fit_diff1 <- diff(fit_part, differences = 1)
adf.test(fit_diff1)

    Augmented Dickey-Fuller Test

data:  fit_diff1
Dickey-Fuller = -1.9581, Lag order = 3, p-value = 0.5888
alternative hypothesis: stationary

fit_diff2 <- diff(fit_part, differences = 2)
adf.test(fit_diff2)

Warning: p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  fit_diff2
Dickey-Fuller = -4.622, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

par(mfrow = c(1,2))

acf(fit_diff2,
    main = "ACF of Second-Differenced Fitted Series")

pacf(fit_diff2,
     main = "PACF of Second-Differenced Fitted Series")

par(mfrow = c(1,1))

# -----------------------------------------------------
# 5.2.1 Model Identification
# -----------------------------------------------------

# 5.2.1.1 Data Splitting (70% Fitted, 30% Hold-out)
fit_part <- window(employment_ts, end = c(2022, 3))
holdout_part <- window(employment_ts, start = c(2022, 4))

# Plot Pembagian Data
plot(employment_ts, main = "Fitted and Hold-out Data Separation", 
     xlab = "Year", ylab = "Employment", col = "black", lwd = 2)
lines(fit_part, col = "blue", lwd = 2)

lines(holdout_part, col = "red", lwd = 2)
legend("topleft", legend = c("Full Series", "Fitted Part (70%)", "Hold-out Part (30%)"), 
       col = c("black", "blue", "red"), lwd = 2)

# 5.2.1.2 Stationary Test (ADF Test pada data Fit)
print("ADF Test: Original Fit Data")

[1] "ADF Test: Original Fit Data"

adf.test(fit_part)

    Augmented Dickey-Fuller Test

data:  fit_part
Dickey-Fuller = -1.7685, Lag order = 3, p-value = 0.6619
alternative hypothesis: stationary

# Melakukan Regular Differencing ke-1 (d=1) karena data memiliki trend
fit_diff1 <- diff(fit_part, differences = 1)
print("ADF Test: First-Differenced Fit Data")

[1] "ADF Test: First-Differenced Fit Data"

adf.test(fit_diff1)

    Augmented Dickey-Fuller Test

data:  fit_diff1
Dickey-Fuller = -1.9581, Lag order = 3, p-value = 0.5888
alternative hypothesis: stationary

# Melakukan Seasonal Differencing ke-1 (D=1, lag=4) untuk membersihkan efek kuartalan
fit_diff1_diff4 <- diff(fit_diff1, lag = 4)
print("ADF Test: Regular + Seasonal Differenced Data")

[1] "ADF Test: Regular + Seasonal Differenced Data"

adf.test(fit_diff1_diff4)

    Augmented Dickey-Fuller Test

data:  fit_diff1_diff4
Dickey-Fuller = -2.1325, Lag order = 2, p-value = 0.5219
alternative hypothesis: stationary

# 5.2.1.3 ACF and PACF Plots (Untuk mengidentifikasi model kandidat)
par(mfrow = c(1,2))
acf(fit_diff1_diff4, main = "ACF of Differenced Fit Series", lag.max = 12)
pacf(fit_diff1_diff4, main = "PACF of Differenced Fit Series", lag.max = 12)

par(mfrow = c(1,1))

# 5.2.2.1 Estimate Three ARIMA Models (menggunakan fit_part)
arima_cand1 <- Arima(fit_part, order = c(0,1,1), seasonal = c(0,1,0), include.drift = TRUE)

Warning: No drift term fitted as the order of difference is 2 or more.

arima_cand2 <- Arima(fit_part, order = c(1,1,0), seasonal = c(0,1,0), include.drift = TRUE)

Warning: No drift term fitted as the order of difference is 2 or more.

arima_cand3 <- Arima(fit_part, order = c(1,1,1), seasonal = c(0,1,0), include.drift = TRUE)

Warning: No drift term fitted as the order of difference is 2 or more.

print("Summary of Candidate Model Parameters 1")

[1] "Summary of Candidate Model Parameters 1"

summary(arima_cand1)

Series: fit_part 
ARIMA(0,1,1)(0,1,0)[4] 

Coefficients:
          ma1
      -0.0976
s.e.   0.1817

sigma^2 = 16920:  log likelihood = -162.93
AIC=329.86   AICc=330.38   BIC=332.37

Training set error measures:
                   ME     RMSE      MAE        MPE      MAPE     MASE        ACF1
Training set 13.65468 116.8131 69.95149 0.08780466 0.4664527 0.258831 -0.02582113

print("Summary of Candidate Model Parameters 2")

[1] "Summary of Candidate Model Parameters 2"

summary(arima_cand2)

Series: fit_part 
ARIMA(1,1,0)(0,1,0)[4] 

Coefficients:
          ar1
      -0.1097
s.e.   0.1930

sigma^2 = 16896:  log likelihood = -162.91
AIC=329.82   AICc=330.34   BIC=332.34

Training set error measures:
                  ME     RMSE      MAE        MPE     MAPE      MASE       ACF1
Training set 13.7262 116.7283 69.56205 0.08826414 0.463901 0.2573901 -0.0135894

print("Summary of Candidate Model Parameters 3")

[1] "Summary of Candidate Model Parameters 3"

summary(arima_cand3)

Series: fit_part 
ARIMA(1,1,1)(0,1,0)[4] 

Coefficients:
          ar1     ma1
      -0.2242  0.1142
s.e.   0.7993  0.7996

sigma^2 = 17587:  log likelihood = -162.9
AIC=331.81   AICc=332.9   BIC=335.58

Training set error measures:
                   ME     RMSE      MAE        MPE      MAPE     MASE        ACF1
Training set 13.66149 116.6866 69.28987 0.08784001 0.4621329 0.256383 -0.01444378

lb_cand1 <- Box.test(residuals(arima_cand1), lag = 8, type = "Ljung-Box")
lb_cand2 <- Box.test(residuals(arima_cand2), lag = 8, type = "Ljung-Box")
lb_cand3 <- Box.test(residuals(arima_cand3), lag = 8, type = "Ljung-Box")

table_comparison_in <- data.frame(
  Model = c("ARIMA(0,1,1)(0,1,0)[4] w/ drift", 
            "ARIMA(1,1,0)(0,1,0)[4] w/ drift", 
            "ARIMA(1,1,1)(0,1,0)[4] w/ drift"),
  AIC  = c(AIC(arima_cand1), AIC(arima_cand2), AIC(arima_cand3)),
  BIC  = c(BIC(arima_cand1), BIC(arima_cand2), BIC(arima_cand3)),
  Ljung_Box_p_value = c(lb_cand1$p.value, lb_cand2$p.value, lb_cand3$p.value)
)

print("Table 5.2.2.2: In-Sample Model Comparison")

[1] "Table 5.2.2.2: In-Sample Model Comparison"

table_comparison_in

h_holdout <- length(holdout_part)

f_cand1 <- forecast(arima_cand1, h = h_holdout)
f_cand2 <- forecast(arima_cand2, h = h_holdout)
f_cand3 <- forecast(arima_cand3, h = h_holdout)

# Menghitung akurasi Out-of-Sample (pilih baris ke-2 [Test set])
acc_cand1 <- accuracy(f_cand1, holdout_part)
acc_cand2 <- accuracy(f_cand2, holdout_part)
acc_cand3 <- accuracy(f_cand3, holdout_part)

table_accuracy_out <- data.frame(
  Model = c("ARIMA(0,1,1)(0,1,0)[4] w/ drift", 
            "ARIMA(1,1,0)(0,1,0)[4] w/ drift", 
            "ARIMA(1,1,1)(0,1,0)[4] w/ drift"),
  RMSE_Test = c(acc_cand1[2, "RMSE"], acc_cand2[2, "RMSE"], acc_cand3[2, "RMSE"]),
  MAE_Test  = c(acc_cand1[2, "MAE"],  acc_cand2[2, "MAE"],  acc_cand3[2, "MAE"]),
  MAPE_Test = c(acc_cand1[2, "MAPE"], acc_cand2[2, "MAPE"], acc_cand3[2, "MAPE"])
)

print("--- Table 5.2.2.3: Out-of-Sample Accuracy (Hold-out Evaluation) ---")

[1] "--- Table 5.2.2.3: Out-of-Sample Accuracy (Hold-out Evaluation) ---"

table_accuracy_out

# 5.2.2.4 Comparison with the Best Univariate Model
# Membandingkan model ARIMA terbaik hasil split data dengan model Holt di Chapter 5(a)
table5224 <- data.frame(
  Model = c("Best ARIMA Model (Split Data)", "Holt's Linear Trend (Full Data)"),
  MSE  = c(acc_cand1[1, "RMSE"]^2, MSE_holt),  # Menggunakan nilai in-sample (Training) agar adil
  MAPE = c(acc_cand1[1, "MAPE"], MAPE_holt)
)

print("Table 5.2.2.4: Comparison with the Best Univariate Model")

[1] "Table 5.2.2.4: Comparison with the Best Univariate Model"

table5224

# 5.2.2.5 Best Forecasting Model & Final Plot
# Berdasarkan laporan, ARIMA(0,1,1)(0,1,0)[4] dengan drift terpilih sebagai model terbaik.
plot(f_cand1, main = "Best ARIMA Model Forecast vs Actual Hold-out Data", 
     xlab = "Year", ylab = "Employment", col = "blue")
lines(holdout_part, col = "red", lwd = 2)

legend("topleft", legend = c("Forecasted", "Actual Hold-out Data"), 
       col = c("blue", "red"), lwd = 2)