Analisis Semiparametrik Regresi Linear–Fourier untuk Pemodelan dan Prediksi Data Meteorologi Harian

Zahra Aura Hisani, Amelia Kusuma Wardani, Zulfa Septilia, Fadilah Andriana Putri Ayunda, M. Ilham Nurul Azka

8 Januari 2026

1 STUDI KASUS

Studi kasus ini menerapkan regresi semiparametrik, yang mengombinasikan regresi linear parametrik dan Deret Fourier nonparametrik, untuk memodelkan dan meramalkan data atmosfer harian. Variabel respon yang digunakan adalah kelembapan spresifik udara pada ketinggian 2 meter(g/kg) (QV2M), sedangkan variabel prediktor meliputi curah hujan harian (mm/hari) (PRECTOTCORR), suhu bola basah pada 2 meter (°C) (T2MWET), tekanan udara permukaan (PS), dan kecepatan angin pada ketinggian 10 meter (m/s) (WS10M).

Regresi linear digunakan untuk memodelkan hubungan linier antara QV2M dan variabel atmosfer, sementara Deret Fourier berfungsi sebagai komponen nonparametrik untuk menangkap pola musiman tahunan pada data deret waktu. Evaluasi dilakukan menggunakan pendekatan rolling forecasting (Xt−1) pada lima observasi terakhir, yang menunjukkan bahwa model mampu menghasilkan prediksi yang stabil dengan tingkat kesalahan yang relatif kecil.

2 TINJAUAN PUSTAKA

2.1 Regresi Linear (Ordinary Least Squares / OLS)

Regresi linear berganda digunakan untuk memodelkan hubungan antara variabel respon \(Y\) dan beberapa variabel penjelas \(X_1, X_2, \ldots, X_p\). Pendugaan parameter regresi dilakukan menggunakan metode Ordinary Least Squares (OLS), yaitu dengan meminimalkan jumlah kuadrat galat antara nilai aktual dan nilai prediksi (Montgomery et al., 2021).

Model regresi linear berganda secara umum dinyatakan sebagai:

\[ Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots + \beta_p X_{pi} + \varepsilon_i \]

Pendugaan parameter dilakukan dengan meminimalkan fungsi objektif:

\[ \min_{\boldsymbol{\beta}} \sum_{i=1}^{n} \left( Y_i - \hat{Y}_i \right)^2 \]

Solusi OLS dalam bentuk matriks adalah:

\[ \hat{\boldsymbol{\beta}} = \left( \mathbf{X}^\top \mathbf{X} \right)^{-1} \mathbf{X}^\top \mathbf{Y} \]

Formulasi tersebut merupakan dasar dari regresi linear klasik dan banyak digunakan dalam analisis meteorologi untuk memodelkan hubungan antarvariabel atmosfer (Kim et al., 2020).

2.2 Deret Fourier untuk Pemodelan Musiman

Deret Fourier digunakan untuk merepresentasikan pola periodik dalam data deret waktu melalui kombinasi fungsi sinus dan kosinus. Pendekatan ini sangat efektif dalam memodelkan pola musiman yang berulang, seperti siklus tahunan pada data meteorologi harian (Bloomfield, 2004).

Komponen Fourier untuk periode \(P\) dan harmonik ke-\(k\) didefinisikan sebagai:

\[ \sin\left(\frac{2\pi k t}{P}\right), \quad \cos\left(\frac{2\pi k t}{P}\right) \]

Deret Fourier memungkinkan pemodelan pola musiman secara kontinu tanpa memerlukan variabel dummy musiman, sehingga lebih efisien dan stabil untuk data deret waktu yang panjang (Shumway & Stoffer, 2017).

2.3 Regresi Linear dengan Komponen Fourier (Harmonic Regression)

Regresi linear dengan komponen Fourier, atau harmonic regression, merupakan pengembangan dari regresi linear dengan menambahkan fungsi sinus dan kosinus untuk menangkap pola musiman. Model ini banyak digunakan dalam analisis iklim dan hidrologi (Bloomfield, 2004).

Bentuk umum model harmonic regression dinyatakan sebagai:

\[ Y_t = \beta_0 + \sum_{j=1}^{p} \beta_j X_{jt} + \sum_{k=1}^{K} \left[ \alpha_k \sin\left(\frac{2\pi k t}{P}\right) + \gamma_k \cos\left(\frac{2\pi k t}{P}\right) \right] + \varepsilon_t \]

Model ini tetap bersifat linear terhadap parameter, sehingga dapat diestimasi menggunakan metode Ordinary Least Squares (OLS). Namun, dibandingkan regresi linear biasa, model ini lebih fleksibel karena mampu menangkap variasi musiman secara eksplisit (Shumway & Stoffer, 2017).

2.4 Evaluasi Kinerja Model Prediksi

Evaluasi kinerja model prediksi dilakukan menggunakan ukuran kesalahan yang umum digunakan dalam analisis peramalan, yaitu Mean Squared Error (MSE), Root Mean Squared Error (RMSE), dan Mean Absolute Percentage Error (MAPE) (Hyndman & Koehler, 2006).

Mean Squared Error (MSE)

Mean Squared Error mengukur rata-rata kuadrat selisih antara nilai aktual dan nilai prediksi, dan dirumuskan sebagai:

\[ \mathrm{MSE} = \frac{1}{n} \sum_{t=1}^{n} \left( Y_t - \hat{Y}_t \right)^2 \]

Root Mean Squared Error (RMSE)

Root Mean Squared Error merupakan akar kuadrat dari MSE, sehingga memiliki satuan yang sama dengan variabel respon:

\[ \mathrm{RMSE} = \sqrt{\mathrm{MSE}} \]

Mean Absolute Percentage Error (MAPE)

Mean Absolute Percentage Error mengukur rata-rata galat absolut dalam bentuk persentase, sehingga lebih mudah diinterpretasikan secara praktis:

\[ \mathrm{MAPE} = \frac{100\%}{n} \sum_{t=1}^{n} \left| \frac{Y_t - \hat{Y}_t}{Y_t} \right| \]

MSE dan RMSE bersifat sensitif terhadap galat berukuran besar, sedangkan MAPE memberikan interpretasi kesalahan dalam bentuk persentase sehingga lebih mudah dipahami dalam konteks evaluasi model peramalan (Hyndman & Koehler, 2006).

2.5 Uji Diagnostik Model

Model regresi linear dan regresi linear dengan komponen musiman berbasis deret Fourier yang diestimasi menggunakan metode Ordinary Least Squares (OLS) memerlukan pengujian asumsi dasar agar hasil estimasi dan inferensi statistik bersifat valid (Montgomery et al., 2021). Oleh karena itu, dilakukan uji diagnostik yang meliputi uji normalitas residual, uji heteroskedastisitas, uji multikolinearitas, dan uji autokorelasi residual.

Uji Normalitas Residual (Shapiro–Wilk)

Uji Shapiro–Wilk digunakan untuk menguji apakah residual model berdistribusi normal (Shapiro & Wilk, 1965). Statistik uji dirumuskan sebagai:

\[ W = \frac{\left( \sum_{i=1}^{n} a_i x_{(i)} \right)^2}{\sum_{i=1}^{n} (x_i - \bar{x})^2} \]

Asumsi normalitas residual penting terutama untuk mendukung validitas inferensi parameter regresi.

Uji Heteroskedastisitas (Breusch–Pagan)

Uji Breusch–Pagan digunakan untuk mendeteksi apakah varians residual bersifat konstan (Breusch & Pagan, 1979). Statistik uji dinyatakan sebagai:

\[ LM = nR^2 \]

Model dikatakan memenuhi asumsi homoskedastisitas apabila tidak terdapat bukti heteroskedastisitas yang signifikan.

Uji Multikolinearitas (Variance Inflation Factor / VIF)

Multikolinearitas diuji menggunakan Variance Inflation Factor (VIF), yang dirumuskan sebagai (Kutner et al., 2005):

\[ \text{VIF}_j = \frac{1}{1 - R_j^2} \]

Nilai VIF yang mendekati satu menunjukkan tidak adanya multikolinearitas yang signifikan antarvariabel bebas.

Uji Autokorelasi Residual (Durbin–Watson)

Uji Durbin–Watson digunakan untuk mendeteksi autokorelasi residual, khususnya pada data deret waktu (Durbin & Watson, 1951). Statistik uji dirumuskan sebagai:

\[ DW = \frac{\sum_{t=2}^{n} (e_t - e_{t-1})^2}{\sum_{t=1}^{n} e_t^2} \]

Nilai statistik Durbin–Watson yang mendekati dua mengindikasikan tidak adanya autokorelasi residual.

3 METODE PENELITIAN

3.1 Library

# Library
library(readxl)
library(readr)
library(dplyr)
library(ggplot2)
library(car)
library(ResourceSelection)
library(caret)
library(knitr)
library(lmtest)
library(np)
library(corrplot)
library(sandwich)

3.2 Input Data

Data yang digunakan merupakan data meteorologi harian yang bersumber dari National Aeronautics and Space Administration (NASA).

# Input data
data <- read_excel(
  "C:/Users/User/Downloads/SEMIPAR/POWER_Point_Daily_20100101_20251130_000d00N_102d00E_LST (2).xlsx"
)

3.3 Pra-pemrosesan Data

Pada tahap pra-pemrosesan data, dilakukan penyiapan data sebelum analisis lebih lanjut. Proses ini meliputi pembentukan variabel tanggal berdasarkan informasi hari dalam tahun (DOY), konversi variabel respon dan variabel penjelas ke dalam bentuk numerik, pengurutan data berdasarkan waktu, serta pemilihan variabel yang relevan untuk analisis. Tahap ini bertujuan untuk memastikan data berada dalam format yang sesuai dan siap digunakan pada tahap eksplorasi dan pemodelan selanjutnya.

data <- data %>%
mutate(
Date = as.Date(DOY - 1, origin = paste0(YEAR, "-01-01")),
Y = as.numeric(QV2M),
X1 = as.numeric(PRECTOTCORR),
X2 = as.numeric(T2MWET),
X3 = as.numeric(PS),
X4 = as.numeric(WS10M)
) %>%
arrange(Date)

dat <- data %>%
select(Date, Y, X1, X2, X3, X4)

head(dat)

## # A tibble: 6 × 6
##   Date           Y    X1    X2    X3    X4
##   <date>     <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2010-01-01  18.0  9.31  24.4  100.  0.53
## 2 2010-01-02  18.4 22.7   24.4  100.  0.67
## 3 2010-01-03  19.3 22.6   25.2  100.  0.53
## 4 2010-01-04  19.1  0.66  25.2  100.  0.58
## 5 2010-01-05  17.6  4.92  24.0  100.  0.7 
## 6 2010-01-06  18.7 14.6   24.8  100.  0.78

str(dat)

## tibble [5,813 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Date: Date[1:5813], format: "2010-01-01" "2010-01-02" ...
##  $ Y   : num [1:5813] 18 18.4 19.3 19.1 17.6 ...
##  $ X1  : num [1:5813] 9.31 22.68 22.6 0.66 4.92 ...
##  $ X2  : num [1:5813] 24.4 24.4 25.2 25.2 24.1 ...
##  $ X3  : num [1:5813] 100 100 100 100 100 ...
##  $ X4  : num [1:5813] 0.53 0.67 0.53 0.58 0.7 0.78 0.66 0.76 0.99 0.69 ...

summary(dat)

##       Date                  Y               X1                X2       
##  Min.   :2010-01-01   Min.   :15.29   Min.   :  0.000   Min.   :21.62  
##  1st Qu.:2013-12-24   1st Qu.:18.53   1st Qu.:  0.910   1st Qu.:24.80  
##  Median :2017-12-16   Median :19.03   Median :  3.780   Median :25.19  
##  Mean   :2017-12-16   Mean   :19.06   Mean   :  7.533   Mean   :25.21  
##  3rd Qu.:2021-12-08   3rd Qu.:19.57   3rd Qu.: 10.370   3rd Qu.:25.63  
##  Max.   :2025-11-30   Max.   :21.96   Max.   :113.770   Max.   :27.27  
##        X3               X4        
##  Min.   : 99.93   Min.   :0.3100  
##  1st Qu.:100.28   1st Qu.:0.6600  
##  Median :100.37   Median :0.7400  
##  Mean   :100.37   Mean   :0.7457  
##  3rd Qu.:100.45   3rd Qu.:0.8300  
##  Max.   :100.83   Max.   :1.3000

3.4 Eksplorasi Data: Scatter Plot dan Korelasi

Analisis ini bertujuan untuk melihat hubungan antara QV2M dan PRECTOTCORR menggunakan scatter plot dan koefisien korelasi Pearson. Scatter plot digunakan untuk mengamati pola hubungan secara visual, sedangkan korelasi Pearson memberikan ukuran numerik kekuatan dan arah hubungan linier.

# Korelasi Y vs X1 

r1 <- cor(data$Y, data$X1, use = "complete.obs", method = "pearson")

ggplot(data, aes(x = X1, y = Y)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(
x = "QV2M (Kelembapan Spesifik)",
y = "PRECTOTCORR (Curah Hujan)",
title = paste0("Korelasi PRECTOTCORR vs QV2M (r = ", round(r1, 3), ")")
) +
theme_minimal()

Interpretasi:

Terdapat hubungan linear positif antara QV2M dan PRECTOTCORR dengan nilai korelasi Pearson sebesar r = 0,208, yang menunjukkan hubungan lemah. Sebaran data yang cukup luas mengindikasikan bahwa hubungan antarvariabel tidak sepenuhnya linear dan dipengaruhi oleh faktor lain.

# Korelasi Y vs X2

r2 <- cor(data$Y, data$X2, use = "complete.obs", method = "pearson")

ggplot(data, aes(x = X2, y = Y)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(
x = "T2MWET (Suhu Bola Basah)",
y = "PRECTOTCORR (Curah Hujan)",
title = paste0("Korelasi PRECTOTCORR vs T2MWET (r = ", round(r2, 3), ")")
) +
theme_minimal()

Interpretasi:

Terdapat hubungan linear positif yang sangat kuat antara QV2M dan T2MWET dengan nilai korelasi Pearson sebesar r = 0,917, yang menunjukkan bahwa suhu bola basah memiliki pengaruh dominan terhadap variasi kelembapan spesifik.

# Korelasi Y vs X3

r3 <- cor(data$Y, data$X3, use = "complete.obs", method = "pearson")

ggplot(data, aes(x = X3, y = Y)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(
x = "PS (Tekanan Permukaan)",
y = "PRECTOTCORR (Curah Hujan)",
title = paste0("Korelasi PRECTOTCORR vs PS (r = ", round(r3, 3), ")")
) +
theme_minimal()

Interpretasi:

Terdapat hubungan linear negatif antara QV2M dan tekanan permukaan (PS) dengan nilai korelasi Pearson sebesar r = −0,300, yang menunjukkan hubungan lemah hingga sedang.

# Korelasi Y vs X4

r4 <- cor(data$Y, data$X4, use = "complete.obs", method = "pearson")

ggplot(data, aes(x = X4, y = Y)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(
x = "WS10M (Kecepatan Angin 10 m)",
y = "PRECTOTCORR (Curah Hujan)",
title = paste0("Korelasi PRECTOTCORR vs WS10M (r = ", round(r4, 3), ")")
) +
theme_minimal()

Interpretasi:

Terdapat hubungan linear negatif yang lemah antara QV2M dan WS10M dengan nilai korelasi Pearson sebesar r = −0,162.

3.4.1 Matriks Korelasi

Pada tahap ini dilakukan analisis matriks korelasi Pearson untuk melihat kekuatan dan arah hubungan linier antara variabel respon QV2M (Y) dan variabel prediktor X1–X4. Matriks korelasi digunakan sebagai ringkasan numerik dari hubungan antarvariabel yang sebelumnya telah dieksplorasi melalui scatter plot.

cor_matrix <- cor(
  dat %>% select(Y, X1, X2, X3, X4),
  use = "complete.obs",
  method = "pearson"
)

round(cor_matrix, 3)

##         Y     X1     X2     X3     X4
## Y   1.000  0.208  0.917 -0.300 -0.162
## X1  0.208  1.000  0.042 -0.043 -0.188
## X2  0.917  0.042  1.000 -0.276 -0.118
## X3 -0.300 -0.043 -0.276  1.000  0.101
## X4 -0.162 -0.188 -0.118  0.101  1.000

corrplot(cor_matrix, method = "number", type = "upper")

Interpretasi:

Matriks korelasi Pearson menunjukkan bahwa T2MWET (X2) memiliki hubungan linier positif sangat kuat dengan QV2M (Y) dengan nilai korelasi sebesar 0,92, sehingga merupakan prediktor utama dalam model. PRECTOTCORR (X1) memiliki korelasi positif lemah terhadap QV2M sebesar 0,21.

Sementara itu, PS (X3) dan WS10M (X4) menunjukkan korelasi negatif lemah dengan QV2M, masing-masing sebesar −0,30 dan −0,16. Korelasi antarvariabel prediktor relatif kecil, sehingga tidak mengindikasikan adanya multikolinearitas yang kuat.

3.5 Split Data (70:30) untuk Model Parametrik

Pada tahap ini, data dibagi menjadi data latih (training) dan data uji (testing) dengan proporsi 70% : 30% secara acak. Pembagian ini bertujuan untuk membangun model parametrik menggunakan data latih dan mengevaluasi kinerja model menggunakan data uji yang tidak terlibat dalam proses pemodelan. Penetapan seed dilakukan untuk memastikan hasil pembagian data dapat direproduksi.

set.seed(123)
n <- nrow(dat)
n_train <- round(0.7 * n)
idx_train <- sample(1:n, size = n_train)

train_par <- dat[idx_train, ]
test_par <- dat[-idx_train, ]

3.6 Estimasi Parametrik (Regresi Linear)

Pada tahap ini dilakukan estimasi model regresi linear untuk memodelkan hubungan antara variabel respon QV2M (Y) dan variabel prediktor X1–X4. Estimasi parameter dilakukan menggunakan data latih hasil pembagian sebelumnya, dengan tujuan memperoleh hubungan linier antarvariabel serta mengevaluasi signifikansi masing-masing prediktor.

mod_par <- lm(
  Y ~ X1 + X2 + X3 + X4,
  data = train_par
)
summary(mod_par)

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = train_par)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.42764 -0.18352  0.01095  0.20535  0.72480 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.2879086  3.7572525   4.867 1.17e-06 ***
## X1           0.0134636  0.0004605  29.238  < 2e-16 ***
## X2           1.1287344  0.0074454 151.601  < 2e-16 ***
## X3          -0.2760044  0.0369281  -7.474 9.47e-14 ***
## X4          -0.1008388  0.0370444  -2.722  0.00651 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2883 on 4064 degrees of freedom
## Multiple R-squared:   0.87,  Adjusted R-squared:  0.8699 
## F-statistic:  6798 on 4 and 4064 DF,  p-value: < 2.2e-16

Interpretasi:

Hasil regresi linear menunjukkan bahwa seluruh variabel prediktor X1, X2, X3, dan X4 berpengaruh signifikan terhadap QV2M pada taraf signifikansi 5%. Variabel X1 dan X2 berpengaruh positif, sedangkan X3 dan X4 berpengaruh negatif terhadap QV2M. Nilai R² sebesar 0,87 menunjukkan bahwa sekitar 87% variasi QV2M dapat dijelaskan oleh model, sehingga model memiliki kecocokan yang sangat baik.

3.7 Uji Signifikansi Parameter (ANOVA)

Pada tahap ini dilakukan uji signifikansi parameter secara simultan menggunakan Analisis Varians (ANOVA). Uji ini bertujuan untuk mengetahui apakah seluruh variabel prediktor dalam model regresi linear secara bersama-sama berpengaruh signifikan terhadap variabel respon QV2M.

mod_par <- lm(
  Y ~ X1 + X2 + X3 + X4,
  data = train_par
)
summary(mod_par)

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = train_par)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.42764 -0.18352  0.01095  0.20535  0.72480 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.2879086  3.7572525   4.867 1.17e-06 ***
## X1           0.0134636  0.0004605  29.238  < 2e-16 ***
## X2           1.1287344  0.0074454 151.601  < 2e-16 ***
## X3          -0.2760044  0.0369281  -7.474 9.47e-14 ***
## X4          -0.1008388  0.0370444  -2.722  0.00651 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2883 on 4064 degrees of freedom
## Multiple R-squared:   0.87,  Adjusted R-squared:  0.8699 
## F-statistic:  6798 on 4 and 4064 DF,  p-value: < 2.2e-16

Interpretasi:

Hasil uji ANOVA menunjukkan bahwa model regresi linear signifikan secara simultan, ditunjukkan oleh nilai F-statistic yang besar dengan p-value < 2.2×10⁻¹⁶. Hal ini mengindikasikan bahwa variabel prediktor X1, X2, X3, dan X4 secara bersama-sama berpengaruh signifikan terhadap variabel respon QV2M. Dengan demikian, model regresi yang dibangun layak digunakan untuk analisis dan pemodelan lebih lanjut.

3.8 Robust Standard Error

Digunakan untuk menghitung dan menguji koefisien regresi dengan robust standard error agar hasil uji statistik tetap valid meskipun terjadi heteroskedastisitas pada data. Dengan metode ini, nilai standar error, t-statistik, dan p-value menjadi lebih reliabel dibandingkan standar error OLS biasa.

coeftest(mod_par, vcov = vcovHC(mod_par, type = "HC1"))

## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) 18.28790859  3.69839331   4.9448 7.929e-07 ***
## X1           0.01346355  0.00064191  20.9742 < 2.2e-16 ***
## X2           1.12873444  0.00780658 144.5876 < 2.2e-16 ***
## X3          -0.27600442  0.03640676  -7.5811 4.219e-14 ***
## X4          -0.10083877  0.03808266  -2.6479  0.008131 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.9 Uji Kebaikan Model (Diagnostik dan Evaluasi Akurasi)

Pada tahap ini dilakukan pengujian kebaikan model regresi linear parametrik. Pengujian meliputi diagnostik asumsi klasik (normalitas, heteroskedastisitas, multikolinearitas, dan autokorelasi) serta evaluasi akurasi prediksi menggunakan MSE, RMSE, dan MAPE.

3.9.1 Diagnostik Residual

Digunakan untuk memeriksa pola residual, normalitas residual, serta kemungkinan adanya pencilan dan pelanggaran asumsi regresi.

par(mfrow = c(2, 2))
plot(mod_par)

par(mfrow = c(1, 1))

Interpretasi:

Berdasarkan plot diagnostik, residual menyebar di sekitar nol tanpa pola yang kuat, meskipun terdapat indikasi nonlinieritas ringan. Plot Q–Q menunjukkan penyimpangan normalitas pada bagian ekor, namun masih dapat ditoleransi karena ukuran sampel yang besar. Plot Scale–Location memperlihatkan varians residual yang relatif homogen, sehingga asumsi homoskedastisitas cukup terpenuhi. Sementara itu, plot Residuals vs Leverage menunjukkan beberapa titik dengan leverage tinggi, namun tidak ada pengamatan yang sangat berpengaruh secara ekstrem.

3.9.2 Uji Normalitas Residual (Shapiro–Wilk)

Uji Shapiro–Wilk digunakan untuk menguji apakah residual model berdistribusi normal.

shapiro_res <- shapiro.test(residuals(mod_par))
shapiro_res

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_par)
## W = 0.99097, p-value = 2.002e-15

Interpretasi:

Hasil uji Shapiro–Wilk menghasilkan nilai p-value < 0,05, sehingga residual tidak berdistribusi normal secara statistik. Namun, mengingat ukuran sampel yang besar, penyimpangan normalitas ini masih dapat ditoleransi dan tidak secara signifikan memengaruhi keandalan estimasi model regresi.

3.9.3 Uji Heteroskedastisitas (Breusch–Pagan)

Uji Breusch–Pagan digunakan untuk mengetahui apakah varians residual bersifat konstan (homoskedastis).

bp_res <- bptest(mod_par)
bp_res

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_par
## BP = 35.015, df = 4, p-value = 4.613e-07

Interpretasi:

Hasil uji Breusch–Pagan menunjukkan nilai p-value < 0,05, sehingga terdapat indikasi heteroskedastisitas pada model regresi. Hal ini menunjukkan bahwa varians residual tidak konstan, sehingga asumsi homoskedastisitas tidak sepenuhnya terpenuhi.

3.9.4 Uji Multikolinearitas (VIF)

Uji multikolinearitas dilakukan menggunakan nilai Variance Inflation Factor (VIF).

vif(mod_par)

##       X1       X2       X3       X4 
## 1.038322 1.086426 1.083333 1.056167

Interpretasi:

Nilai VIF untuk seluruh variabel prediktor (X1–X4) berada sangat dekat dengan 1, yaitu antara 1,04–1,09. Hal ini menunjukkan tidak terdapat multikolinearitas yang signifikan antarvariabel bebas. Dengan demikian, hubungan antarvariabel prediktor relatif independen dan estimasi parameter regresi dapat dianggap stabil dan reliabel.

3.9.5 Uji Autokorelasi (Durbin–Watson)

Uji Durbin–Watson digunakan untuk mendeteksi adanya autokorelasi pada residual model.

dwtest(mod_par)

## 
##  Durbin-Watson test
## 
## data:  mod_par
## DW = 1.9818, p-value = 0.2802
## alternative hypothesis: true autocorrelation is greater than 0

Interpretasi:

Hasil uji Durbin–Watson menunjukkan nilai p-value > 0,05, sehingga tidak terdapat indikasi autokorelasi pada residual model. Dengan demikian, asumsi independensi residual dapat dianggap terpenuhi.

3.9.6 Evaluasi Akurasi Model (MSE, RMSE, dan MAPE)

Evaluasi dilakukan untuk mengukur kemampuan prediksi model pada data latih dan data uji.

mape_safe <- function(actual, pred, eps = 1e-8) {
  ok <- !is.na(actual) & !is.na(pred) & abs(actual) > eps
  if (sum(ok) == 0) return(NA_real_)
  mean(abs((actual[ok] - pred[ok]) / actual[ok])) * 100
}

# Prediksi
pred_train_par <- predict(mod_par, newdata = train_par)
pred_test_par  <- predict(mod_par, newdata = test_par)

# MSE
MSE_train_par <- mean((train_par$Y - pred_train_par)^2, na.rm = TRUE)
MSE_test_par  <- mean((test_par$Y - pred_test_par)^2, na.rm = TRUE)

# RMSE
RMSE_train_par <- sqrt(MSE_train_par)
RMSE_test_par  <- sqrt(MSE_test_par)

# MAPE
MAPE_train_par <- mape_safe(train_par$Y, pred_train_par)
MAPE_test_par  <- mape_safe(test_par$Y, pred_test_par)

MSE_train_par; MSE_test_par

## [1] 0.0830059

## [1] 0.08101609

RMSE_train_par; RMSE_test_par

## [1] 0.2881075

## [1] 0.2846333

MAPE_train_par; MAPE_test_par

## [1] 1.212907

## [1] 1.184458

Interpretasi:

Berdasarkan hasil evaluasi akurasi model, diperoleh nilai MSE pada data latih sebesar 0,083 dan pada data uji sebesar 0,081. Nilai yang relatif kecil dan hampir sama ini menunjukkan bahwa kesalahan prediksi model tergolong rendah serta tidak terdapat indikasi overfitting.

Nilai RMSE pada data latih (0,288) dan data uji (0,285) juga menunjukkan konsistensi yang baik, sehingga model memiliki kestabilan dalam memprediksi data baru.

Selanjutnya, nilai MAPE pada data latih sebesar 1,21% dan pada data uji sebesar 1,18% tergolong sangat kecil. Hal ini mengindikasikan bahwa secara persentase, kesalahan prediksi model sangat rendah.

Secara keseluruhan, model regresi parametrik memiliki tingkat akurasi yang baik dan kemampuan generalisasi yang memadai terhadap data uji.

3.10 Estimasi Deret Fourier & Evaluasi (Testing = 5 data terakhir)

Pada tahap ini dilakukan pemodelan menggunakan deret Fourier untuk menangkap pola musiman pada data deret waktu. Komponen Fourier digunakan sebagai pendekatan nonparametrik, sedangkan variabel prediktor cuaca dimodelkan secara parametrik melalui regresi linear. Evaluasi dilakukan dengan menggunakan 5 observasi terakhir sebagai data uji, sehingga mencerminkan skenario peramalan aktual.

3.10.1 Pembentukan Data Deret Waktu dan Komponen Fourier

Data terlebih dahulu diurutkan berdasarkan waktu untuk menjaga konsistensi struktur deret waktu. Selanjutnya, dibentuk indeks waktu dan komponen Fourier berupa fungsi sinus dan kosinus dengan periode tahunan (365 hari) dan jumlah harmonik tertentu.

# Pastikan urut waktu
dat_ts <- dat %>% arrange(Date)

# Periode musiman (harian -> tahunan)
P <- 365

# Jumlah harmonik
K <- 3

# Indeks waktu
dat_ts <- dat_ts %>% mutate(t = row_number())

# Fitur Fourier
for (k in 1:K) {
  dat_ts[[paste0("sin", k)]] <- sin(2 * pi * k * dat_ts$t / P)
  dat_ts[[paste0("cos", k)]] <- cos(2 * pi * k * dat_ts$t / P)
}

3.10.2 Pembagian Data Latih dan Data Uji

Data dibagi dengan menyisakan 5 observasi terakhir sebagai data uji, sementara data sebelumnya digunakan sebagai data latih. Pendekatan ini umum digunakan dalam analisis deret waktu untuk mengevaluasi kemampuan peramalan satu langkah ke depan.

# Split data: sisakan 5 observasi terakhir
n2 <- nrow(dat_ts)
n_test <- 5
n_train2 <- n2 - n_test

train_fourier <- dat_ts[1:n_train2, ]
test_fourier  <- dat_ts[(n_train2 + 1):n2, ]

3.10.3 Spesifikasi Model Fourier

Model dibentuk dengan mengombinasikan variabel prediktor parametrik (X1–X4) dan komponen Fourier nonparametrik.

# Formula Fourier otomatis
fourier_terms <- paste(
  c(paste0("sin", 1:K), paste0("cos", 1:K)),
  collapse = " + "
)

form_fourier <- as.formula(
  paste("Y ~ X1 + X2 + X3 + X4 +", fourier_terms)
)

3.10.4 Peramalan dengan Rolling / Expanding Window (Xt−1)

Untuk mensimulasikan kondisi peramalan aktual, digunakan pendekatan rolling (expanding window). Pada setiap langkah, model diestimasi ulang menggunakan data latih yang terus bertambah, kemudian digunakan untuk memprediksi satu periode ke depan. Nilai variabel prediktor yang digunakan saat prediksi adalah nilai pada waktu sebelumnya (Xt−1).

rolling_fourier <- function(data_ts, form_fourier, n_test) {
  
  n_total <- nrow(data_ts)
  hasil <- list()
  
  for (i in 1:n_test) {
    
    # Data training bertambah (expanding window)
    train_i <- data_ts[1:(n_total - n_test + i - 1), ]
    test_i  <- data_ts[n_total - n_test + i, ]
    
    # Estimasi ulang model
    mod_i <- lm(form_fourier, data = train_i)
    
    # Gunakan Xt-1
    prev_row <- data_ts[n_total - n_test + i - 1, ]
    test_i$X1 <- prev_row$X1
    test_i$X2 <- prev_row$X2
    test_i$X3 <- prev_row$X3
    test_i$X4 <- prev_row$X4
    
    # Prediksi
    yhat <- predict(mod_i, newdata = test_i)
    
    # Simpan hasil
    hasil[[i]] <- data.frame(
      Date = test_i$Date,
      Y_Aktual = test_i$Y,
      Y_Prediksi = yhat
    )
    
    # Tampilkan model
    cat("====================================\n")
    cat("Model ke-", i, "(Tanggal:", as.character(test_i$Date), ")\n")
    print_model(mod_i)
  }
  
  return(do.call(rbind, hasil))
}

mod_fourier <- lm(form_fourier, data = train_fourier)

summary(mod_fourier)

## 
## Call:
## lm(formula = form_fourier, data = train_fourier)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35124 -0.17275  0.01251  0.18610  0.80076 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.5513481  2.9706421   2.205   0.0275 *  
## X1           0.0135778  0.0003653  37.166  < 2e-16 ***
## X2           1.0984612  0.0064385 170.609  < 2e-16 ***
## X3          -0.1518751  0.0292257  -5.197 2.10e-07 ***
## X4          -0.0495190  0.0298329  -1.660   0.0970 .  
## sin1         0.0786077  0.0051130  15.374  < 2e-16 ***
## sin2        -0.0845308  0.0054477 -15.517  < 2e-16 ***
## sin3        -0.0123528  0.0049966  -2.472   0.0135 *  
## cos1        -0.0428596  0.0052156  -8.218 2.54e-16 ***
## cos2         0.0831836  0.0051223  16.239  < 2e-16 ***
## cos3         0.0466504  0.0049692   9.388  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2672 on 5797 degrees of freedom
## Multiple R-squared:  0.8889, Adjusted R-squared:  0.8887 
## F-statistic:  4639 on 10 and 5797 DF,  p-value: < 2.2e-16

anova(mod_fourier)

## Analysis of Variance Table
## 
## Response: Y
##             Df  Sum Sq Mean Sq    F value    Pr(>F)    
## X1           1  160.69  160.69  2250.9251 < 2.2e-16 ***
## X2           1 3077.80 3077.80 43114.7252 < 2.2e-16 ***
## X3           1    7.01    7.01    98.1513 < 2.2e-16 ***
## X4           1    1.66    1.66    23.2698 1.444e-06 ***
## sin1         1   15.49   15.49   216.9433 < 2.2e-16 ***
## sin2         1   17.61   17.61   246.6440 < 2.2e-16 ***
## sin3         1    0.40    0.40     5.5828   0.01817 *  
## cos1         1    5.78    5.78    80.9578 < 2.2e-16 ***
## cos2         1   18.94   18.94   265.2530 < 2.2e-16 ***
## cos3         1    6.29    6.29    88.1321 < 2.2e-16 ***
## Residuals 5797  413.83    0.07                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretasi:

Model regresi dengan penambahan komponen Deret Fourier menunjukkan kinerja yang sangat baik, dengan nilai Adjusted R-squared sebesar 0,8887, yang berarti sekitar 88,9% variasi Y dapat dijelaskan oleh model.

Variabel X1 dan X2 berpengaruh positif dan signifikan terhadap Y, sedangkan X3 berpengaruh negatif dan signifikan. Variabel X4 memiliki pengaruh yang lebih lemah. Komponen Fourier (sin dan cos) hingga harmonik ke-3 sebagian besar signifikan, yang menandakan adanya pola musiman tahunan yang kuat dalam data.

Hasil uji ANOVA menunjukkan bahwa model secara keseluruhan signifikan (p-value < 0,05). Dengan demikian, model Deret Fourier yang bersifat semiparametrik ini mampu menangkap pola linier dan musiman secara simultan, serta layak digunakan untuk peramalan.

3.10.5 Evaluasi n (5) data testing (Xt-1) - LOOPING DINAMIS

Evaluasi dilakukan dengan menghitung nilai MSE, RMSE, dan MAPE berdasarkan hasil prediksi satu langkah ke depan menggunakan pendekatan Xt−1.

evaluate_lag1 <- function(model, data_ts, n_test) {
  
  cat("Model yang digunakan untuk evaluasi testing:\n")
  print_model(model)
  
  n_total <- nrow(data_ts)
  n_train <- n_total - n_test
  
  pred <- numeric(n_test)
  
  for (i in 1:n_test) {
    prev_row <- data_ts[n_train + i - 1, ]
    curr_row <- data_ts[n_train + i, ]
    
    newdata <- curr_row
    newdata$X1 <- prev_row$X1
    newdata$X2 <- prev_row$X2
    newdata$X3 <- prev_row$X3
    newdata$X4 <- prev_row$X4
    
    pred[i] <- predict(model, newdata = newdata)
  }
  
  hasil <- data.frame(
    Date = data_ts$Date[(n_train + 1):n_total],
    Y_Aktual = data_ts$Y[(n_train + 1):n_total],
    Y_Prediksi = pred
  )
  
  hasil$Error <- hasil$Y_Aktual - hasil$Y_Prediksi
  
  list(
    hasil = hasil,
    MSE  = mean(hasil$Error^2),
    RMSE = sqrt(mean(hasil$Error^2)),
    MAPE = mape_safe(hasil$Y_Aktual, hasil$Y_Prediksi)
  )
}

3.10.6 Rolling Fourier untuk 5 Data Testing

Pada bagian ini dilakukan evaluasi kemampuan prediksi model regresi linear dengan komponen Fourier menggunakan pendekatan rolling forecasting pada 5 data uji terakhir. Metode rolling dilakukan dengan cara memperbarui jendela data latih secara bertahap dan menghasilkan prediksi satu langkah ke depan (one-step ahead forecast) pada setiap periode pengujian. Pendekatan ini bertujuan untuk menilai kestabilan dan konsistensi kinerja model dalam memprediksi data yang belum diamati, serta mencerminkan kondisi peramalan yang lebih realistis pada data deret waktu.

# Rolling
rolling_result <- rolling_fourier(
  data_ts = dat_ts,
  form_fourier = form_fourier,
  n_test = 5
)

## ====================================
## Model ke- 1 (Tanggal: 2025-11-26 )
## 
## Call:
## lm(formula = form_fourier, data = train_i)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35124 -0.17275  0.01251  0.18610  0.80076 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.5513481  2.9706421   2.205   0.0275 *  
## X1           0.0135778  0.0003653  37.166  < 2e-16 ***
## X2           1.0984612  0.0064385 170.609  < 2e-16 ***
## X3          -0.1518751  0.0292257  -5.197 2.10e-07 ***
## X4          -0.0495190  0.0298329  -1.660   0.0970 .  
## sin1         0.0786077  0.0051130  15.374  < 2e-16 ***
## sin2        -0.0845308  0.0054477 -15.517  < 2e-16 ***
## sin3        -0.0123528  0.0049966  -2.472   0.0135 *  
## cos1        -0.0428596  0.0052156  -8.218 2.54e-16 ***
## cos2         0.0831836  0.0051223  16.239  < 2e-16 ***
## cos3         0.0466504  0.0049692   9.388  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2672 on 5797 degrees of freedom
## Multiple R-squared:  0.8889, Adjusted R-squared:  0.8887 
## F-statistic:  4639 on 10 and 5797 DF,  p-value: < 2.2e-16
## 
## ====================================
## Model ke- 2 (Tanggal: 2025-11-27 )
## 
## Call:
## lm(formula = form_fourier, data = train_i)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35113 -0.17294  0.01283  0.18608  0.80064 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.4932297  2.9708195   2.186   0.0289 *  
## X1           0.0135801  0.0003654  37.167  < 2e-16 ***
## X2           1.0985686  0.0064390 170.612  < 2e-16 ***
## X3          -0.1513084  0.0292275  -5.177 2.33e-07 ***
## X4          -0.0516140  0.0298082  -1.732   0.0834 .  
## sin1         0.0786397  0.0051137  15.378  < 2e-16 ***
## sin2        -0.0842990  0.0054465 -15.478  < 2e-16 ***
## sin3        -0.0121803  0.0049961  -2.438   0.0148 *  
## cos1        -0.0429687  0.0052159  -8.238  < 2e-16 ***
## cos2         0.0831717  0.0051230  16.235  < 2e-16 ***
## cos3         0.0466464  0.0049699   9.386  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2672 on 5798 degrees of freedom
## Multiple R-squared:  0.8889, Adjusted R-squared:  0.8887 
## F-statistic:  4639 on 10 and 5798 DF,  p-value: < 2.2e-16
## 
## ====================================
## Model ke- 3 (Tanggal: 2025-11-28 )
## 
## Call:
## lm(formula = form_fourier, data = train_i)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35121 -0.17290  0.01276  0.18664  0.80067 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.4990007  2.9706934   2.188   0.0287 *  
## X1           0.0135792  0.0003654  37.167  < 2e-16 ***
## X2           1.0985042  0.0064381 170.624  < 2e-16 ***
## X3          -0.1513559  0.0292262  -5.179 2.31e-07 ***
## X4          -0.0507277  0.0297829  -1.703   0.0886 .  
## sin1         0.0786283  0.0051135  15.377  < 2e-16 ***
## sin2        -0.0844138  0.0054440 -15.506  < 2e-16 ***
## sin3        -0.0122615  0.0049947  -2.455   0.0141 *  
## cos1        -0.0429264  0.0052154  -8.231 2.28e-16 ***
## cos2         0.0831796  0.0051228  16.237  < 2e-16 ***
## cos3         0.0466538  0.0049697   9.388  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2672 on 5799 degrees of freedom
## Multiple R-squared:  0.8889, Adjusted R-squared:  0.8887 
## F-statistic:  4640 on 10 and 5799 DF,  p-value: < 2.2e-16
## 
## ====================================
## Model ke- 4 (Tanggal: 2025-11-29 )
## 
## Call:
## lm(formula = form_fourier, data = train_i)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35142 -0.17286  0.01272  0.18682  0.80064 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.4786792  2.9708159   2.181   0.0292 *  
## X1           0.0135802  0.0003654  37.167  < 2e-16 ***
## X2           1.0983620  0.0064375 170.618  < 2e-16 ***
## X3          -0.1511225  0.0292273  -5.171 2.41e-07 ***
## X4          -0.0500058  0.0297793  -1.679   0.0932 .  
## sin1         0.0786072  0.0051137  15.372  < 2e-16 ***
## sin2        -0.0845911  0.0054426 -15.542  < 2e-16 ***
## sin3        -0.0123928  0.0049939  -2.482   0.0131 *  
## cos1        -0.0428643  0.0052154  -8.219 2.51e-16 ***
## cos2         0.0832107  0.0051230  16.243  < 2e-16 ***
## cos3         0.0466745  0.0049699   9.391  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2672 on 5800 degrees of freedom
## Multiple R-squared:  0.8889, Adjusted R-squared:  0.8887 
## F-statistic:  4640 on 10 and 5800 DF,  p-value: < 2.2e-16
## 
## ====================================
## Model ke- 5 (Tanggal: 2025-11-30 )
## 
## Call:
## lm(formula = form_fourier, data = train_i)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35113 -0.17305  0.01262  0.18687  0.80063 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.5239442  2.9704812   2.196   0.0281 *  
## X1           0.0135780  0.0003654  37.162  < 2e-16 ***
## X2           1.0983830  0.0064375 170.622  < 2e-16 ***
## X3          -0.1515828  0.0292238  -5.187 2.21e-07 ***
## X4          -0.0495045  0.0297752  -1.663   0.0964 .  
## sin1         0.0786519  0.0051135  15.381  < 2e-16 ***
## sin2        -0.0845150  0.0054421 -15.530  < 2e-16 ***
## sin3        -0.0123016  0.0049931  -2.464   0.0138 *  
## cos1        -0.0429327  0.0052150  -8.233 2.24e-16 ***
## cos2         0.0831454  0.0051226  16.231  < 2e-16 ***
## cos3         0.0466580  0.0049699   9.388  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2672 on 5801 degrees of freedom
## Multiple R-squared:  0.8889, Adjusted R-squared:  0.8887 
## F-statistic:  4640 on 10 and 5801 DF,  p-value: < 2.2e-16

rolling_result

##          Date Y_Aktual Y_Prediksi
## 1  2025-11-26    18.06   18.83710
## 11 2025-11-27    18.43   18.49098
## 12 2025-11-28    18.40   18.23664
## 13 2025-11-29    18.47   18.06175
## 14 2025-11-30    18.29   18.74304

Interpretasi:

Evaluasi menggunakan rolling forecasting pada lima data testing terakhir menunjukkan bahwa prediksi model Deret Fourier mendekati nilai aktual dengan kesalahan yang relatif kecil. Nilai MSE, RMSE, dan MAPE yang rendah menandakan bahwa model mampu menangkap pola musiman data dengan baik dan layak digunakan untuk peramalan jangka pendek.

# Evaluasi Prediksi Dinamis (Lag-1)
eval_5 <- evaluate_lag1(
  model   = mod_fourier,
  data_ts = dat_ts,
  n_test  = 5
)

## Model yang digunakan untuk evaluasi testing:
## 
## Call:
## lm(formula = form_fourier, data = train_fourier)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35124 -0.17275  0.01251  0.18610  0.80076 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.5513481  2.9706421   2.205   0.0275 *  
## X1           0.0135778  0.0003653  37.166  < 2e-16 ***
## X2           1.0984612  0.0064385 170.609  < 2e-16 ***
## X3          -0.1518751  0.0292257  -5.197 2.10e-07 ***
## X4          -0.0495190  0.0298329  -1.660   0.0970 .  
## sin1         0.0786077  0.0051130  15.374  < 2e-16 ***
## sin2        -0.0845308  0.0054477 -15.517  < 2e-16 ***
## sin3        -0.0123528  0.0049966  -2.472   0.0135 *  
## cos1        -0.0428596  0.0052156  -8.218 2.54e-16 ***
## cos2         0.0831836  0.0051223  16.239  < 2e-16 ***
## cos3         0.0466504  0.0049692   9.388  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2672 on 5797 degrees of freedom
## Multiple R-squared:  0.8889, Adjusted R-squared:  0.8887 
## F-statistic:  4639 on 10 and 5797 DF,  p-value: < 2.2e-16

eval_5$hasil

##         Date Y_Aktual Y_Prediksi       Error
## 1 2025-11-26    18.06   18.83710 -0.77709998
## 2 2025-11-27    18.43   18.49245 -0.06245282
## 3 2025-11-28    18.40   18.23740  0.16259954
## 4 2025-11-29    18.47   18.06155  0.40845087
## 5 2025-11-30    18.29   18.74317 -0.45317015

Interpretasi:

Hasil evaluasi rolling testing dengan skema prediksi dinamis Lag-1 menunjukkan bahwa model Deret Fourier mampu menghasilkan prediksi yang mendekati nilai aktual pada 5 data testing terakhir. Nilai error relatif kecil dan tidak menunjukkan pola kesalahan yang konsisten. Nilai MSE, RMSE, dan MAPE yang rendah mengindikasikan bahwa model memiliki akurasi yang baik dan stabil dalam memprediksi data secara rolling.

# Nilai Akurasi Rolling Forecast
eval_5$MSE

## [1] 0.2012837

eval_5$RMSE

## [1] 0.4486466

eval_5$MAPE

## [1] 2.042912

Interpretasi:

Hasil evaluasi rolling forecast menunjukkan nilai MSE sebesar 0,2013, RMSE 0,4486, dan MAPE 2,04%. Nilai kesalahan yang relatif kecil ini menandakan bahwa model Deret Fourier memiliki akurasi yang baik dalam memprediksi 5 data testing terakhir dan mampu mengikuti pola data secara stabil.

# Evaluasi Akurasi Model pada Data Latih dan Data Uji
# Prediksi
pred_train_fourier <- predict(mod_fourier, newdata = train_fourier)
pred_test_fourier  <- predict(mod_fourier, newdata = test_fourier)

# MSE
MSE_train_fourier <- mean((train_fourier$Y - pred_train_fourier)^2, na.rm = TRUE)
MSE_test_fourier  <- mean((test_fourier$Y  - pred_test_fourier)^2,  na.rm = TRUE)

# RMSE
RMSE_train_fourier <- sqrt(MSE_train_fourier)
RMSE_test_fourier  <- sqrt(MSE_test_fourier)

# MAPE
MAPE_train_fourier <- mape_safe(train_fourier$Y, pred_train_fourier)
MAPE_test_fourier  <- mape_safe(test_fourier$Y,  pred_test_fourier)

MSE_train_fourier; MSE_test_fourier

## [1] 0.07125101

## [1] 0.1029486

RMSE_train_fourier; RMSE_test_fourier

## [1] 0.2669289

## [1] 0.3208561

MAPE_train_fourier; MAPE_test_fourier

## [1] 1.117683

## [1] 1.702283

Interpretasi:

Hasil evaluasi menunjukkan bahwa model Deret Fourier memiliki kinerja prediksi yang baik dan stabil. Pada data latih diperoleh nilai MSE = 0,071, RMSE = 0,267, dan MAPE = 1,12%, sedangkan pada data uji diperoleh MSE = 0,103, RMSE = 0,321, dan MAPE = 1,70%.

Nilai kesalahan pada data uji yang sedikit lebih besar dibandingkan data latih masih dalam batas wajar, sehingga tidak menunjukkan overfitting. Nilai MAPE yang rendah (di bawah 2%) mengindikasikan bahwa model mampu memprediksi dengan tingkat akurasi yang tinggi baik pada data latih maupun data uji.

3.11 Tabel Ringkasan Evaluasi Model

Untuk membandingkan kinerja model Parametrik (Linear Regression) dan model Deret Fourier, dilakukan peringkasan nilai akurasi pada data latih dan data uji. Ukuran evaluasi yang digunakan meliputi Mean Squared Error (MSE), Root Mean Squared Error (RMSE), dan Mean Absolute Percentage Error (MAPE). Tabel ini digunakan untuk menilai model dengan kesalahan prediksi paling kecil dan kinerja paling stabil.

eval_table <- data.frame(
  Model = c(
    "Parametrik (LM)",
    paste0("Fourier (LM + Fourier, K = ", K, ")")
  ),
  MSE_train  = c(MSE_train_par,  MSE_train_fourier),
  RMSE_train = c(RMSE_train_par, RMSE_train_fourier),
  MAPE_train = c(MAPE_train_par, MAPE_train_fourier),
  MSE_test   = c(MSE_test_par,   MSE_test_fourier),
  RMSE_test  = c(RMSE_test_par,  RMSE_test_fourier),
  MAPE_test  = c(MAPE_test_par,  MAPE_test_fourier)
)

eval_table

##                           Model  MSE_train RMSE_train MAPE_train   MSE_test
## 1               Parametrik (LM) 0.08300590  0.2881075   1.212907 0.08101609
## 2 Fourier (LM + Fourier, K = 3) 0.07125101  0.2669289   1.117683 0.10294861
##   RMSE_test MAPE_test
## 1 0.2846333  1.184458
## 2 0.3208561  1.702283

Interpretasi:

Berdasarkan Tabel Ringkasan Evaluasi Model, model Deret Fourier menunjukkan kinerja yang lebih baik pada data latih, ditunjukkan oleh nilai MSE, RMSE, dan MAPE yang lebih kecil dibandingkan model parametrik. Hal ini menandakan bahwa penambahan komponen Fourier mampu menangkap pola musiman data dengan lebih baik.

Pada data uji, model parametrik memiliki nilai kesalahan yang sedikit lebih kecil dibandingkan model Fourier. Hal ini menunjukkan bahwa meskipun model Fourier lebih akurat pada data latih, model parametrik cenderung lebih stabil dalam memprediksi data uji. Secara keseluruhan, model Fourier unggul dalam pemodelan pola, sedangkan model parametrik menunjukkan kemampuan generalisasi yang lebih baik.

3.12 Plot Aktual vs Prediksi dan Prediksi Satu Periode ke Depan

Pada bagian ini dilakukan visualisasi perbandingan antara nilai aktual dan nilai prediksi hasil model Deret Fourier menggunakan skema rolling forecasting (Xt−1) pada 5 data testing terakhir. Plot ini bertujuan untuk melihat secara visual kedekatan pola prediksi terhadap data aktual.

Selain itu, dilakukan pula prediksi satu periode ke depan menggunakan model Fourier yang telah diestimasi, dengan memanfaatkan informasi observasi terakhir (Xt−1). Hasil prediksi ini memberikan gambaran kemampuan model dalam melakukan peramalan jangka pendek.

# Plot Aktual vs Prediksi (5 Data Testing, Xt−1)
# Ambil prediksi lag-1 dari evaluasi
pred_test_lag1 <- eval_5$hasil$Y_Prediksi

plot(
  test_fourier$Date,
  test_fourier$Y,
  type = "l",
  col = "black",
  lwd = 2,
  xlab = "Date",
  ylab = "Y",
  main = "Aktual vs Prediksi (Fourier K = 3, Xt-1)"
)

lines(
  test_fourier$Date,
  pred_test_lag1,
  col = "red",
  lwd = 2,
  lty = 2
)

legend(
  "topleft",
  legend = c("Aktual", "Prediksi"),
  col = c("black", "red"),
  lty = c(1, 2),
  lwd = 2,
  bty = "n"
)

Interpretasi:

Grafik menunjukkan bahwa prediksi model Deret Fourier (K = 3) dengan pendekatan rolling forecasting (Xt-1) mampu mengikuti pola pergerakan data aktual pada 5 data uji terakhir. Garis prediksi secara umum berada dekat dengan nilai aktual, meskipun terdapat deviasi pada beberapa periode, khususnya saat terjadi perubahan nilai yang relatif tajam. Hal ini menunjukkan bahwa model cukup baik dalam menangkap pola musiman dan dinamika jangka pendek, namun masih memiliki keterbatasan dalam memprediksi fluktuasi ekstrem pada data uji.

# Prediksi Satu Periode ke Depan (Xt−1)
predict_next_lag1 <- function(model, data_ts, K, P) {
  
  last_row <- data_ts[nrow(data_ts), ]
  next_t <- last_row$t + 1
  next_date <- last_row$Date + 1
  
  newdata <- last_row
  newdata$t <- next_t
  newdata$Date <- next_date
  newdata$Y <- NA_real_
  
  for (k in 1:K) {
    newdata[[paste0("sin", k)]] <- sin(2 * pi * k * next_t / P)
    newdata[[paste0("cos", k)]] <- cos(2 * pi * k * next_t / P)
  }
  
  yhat <- predict(model, newdata = newdata)
  
  cat("Model yang digunakan:\n")
  print_model(model)
  
  data.frame(
    Date = next_date,
    Y_Prediksi = yhat
  )
}

pred_1day <- predict_next_lag1(mod_fourier, dat_ts, K, P)

## Model yang digunakan:
## 
## Call:
## lm(formula = form_fourier, data = train_fourier)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35124 -0.17275  0.01251  0.18610  0.80076 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.5513481  2.9706421   2.205   0.0275 *  
## X1           0.0135778  0.0003653  37.166  < 2e-16 ***
## X2           1.0984612  0.0064385 170.609  < 2e-16 ***
## X3          -0.1518751  0.0292257  -5.197 2.10e-07 ***
## X4          -0.0495190  0.0298329  -1.660   0.0970 .  
## sin1         0.0786077  0.0051130  15.374  < 2e-16 ***
## sin2        -0.0845308  0.0054477 -15.517  < 2e-16 ***
## sin3        -0.0123528  0.0049966  -2.472   0.0135 *  
## cos1        -0.0428596  0.0052156  -8.218 2.54e-16 ***
## cos2         0.0831836  0.0051223  16.239  < 2e-16 ***
## cos3         0.0466504  0.0049692   9.388  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2672 on 5797 degrees of freedom
## Multiple R-squared:  0.8889, Adjusted R-squared:  0.8887 
## F-statistic:  4639 on 10 and 5797 DF,  p-value: < 2.2e-16

pred_1day

##         Date Y_Prediksi
## 1 2025-12-01   18.61437

Interpretasi:

Berdasarkan hasil prediksi satu periode ke depan menggunakan model Deret Fourier (K = 3, Xt-1), diperoleh nilai prediksi sebesar 18,61 untuk tanggal 1 Desember 2025. Prediksi ini dihasilkan dari model yang telah teruji baik pada data latih dan data uji, sehingga mencerminkan pola musiman dan dinamika data historis. Dengan demikian, model Fourier mampu memberikan estimasi jangka pendek yang stabil dan dapat digunakan untuk peramalan satu periode ke depan.

4 KESIMPULAN

Berdasarkan hasil penelitian, dapat disimpulkan bahwa model regresi linear mampu menjelaskan hubungan linier antara variabel atmosfer, namun kurang optimal dalam menangkap pola musiman yang terdapat pada data deret waktu.

Penambahan komponen Deret Fourier sebagai bagian nonparametrik terbukti meningkatkan kemampuan model dalam merepresentasikan pola musiman tahunan. Model regresi semiparametrik (linear–Fourier) menunjukkan kinerja prediksi yang baik dan stabil, baik pada data latih maupun data uji, dengan tingkat kesalahan yang rendah berdasarkan nilai MSE, RMSE, dan MAPE.

Pendekatan rolling forecasting memperlihatkan bahwa model mampu melakukan peramalan satu langkah ke depan secara konsisten dan realistis. Oleh karena itu, regresi semiparametrik berbasis Deret Fourier nonparametrik merupakan metode yang efektif untuk pemodelan dan peramalan data atmosfer harian yang memiliki pola musiman kuat.

5 DAFTAR PUSTAKA

Bloomfield, P. (2004). Fourier analysis of time series: An introduction (2nd ed.).
John Wiley & Sons.
Breusch, T. S., & Pagan, A. R. (1979).
A simple test for heteroskedasticity and random coefficient variation.
Econometrica, 47(5), 1287–1294.
https://doi.org/10.2307/1911963
Durbin, J., & Watson, G. S. (1951).
Testing for serial correlation in least squares regression, I.
Biometrika, 38(1–2), 159–178.
https://doi.org/10.1093/biomet/38.1-2.159
Hyndman, R. J., & Koehler, A. B. (2006).
Another look at measures of forecast accuracy.
International Journal of Forecasting, 22(4), 679–688.
https://doi.org/10.1016/j.ijforecast.2006.03.001
Kim, S. W., Singh, V. P., & Park, J. (2020).
Development of a multiple linear regression model for drought index estimation.
Water, 12(12), 3393.
https://doi.org/10.3390/w12123393
Kutner, M. H., Nachtsheim, C. J., Neter, J., & Li, W. (2005).
Applied linear statistical models (5th ed.).
McGraw-Hill/Irwin.
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2021).
Introduction to linear regression analysis (6th ed.).
John Wiley & Sons.
Shapiro, S. S., & Wilk, M. B. (1965).
An analysis of variance test for normality (complete samples).
Biometrika, 52(3–4), 591–611.
https://doi.org/10.1093/biomet/52.3-4.591
Shumway, R. H., & Stoffer, D. S. (2017).
Time series analysis and its applications: With R examples (4th ed.).
Springer.
https://doi.org/10.1007/978-3-319-52452-8