# Import data
setwd("D:/IPB UNIVERSITY/PERKULIAHAN/Semester 6/Analisis Model Prediktif/Tugas Akhir/Data AMP")
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
faktor <- read_excel("kelompok_8_data_variabel.xlsx")
View(faktor)
# Sesuaikan format desimal
faktor$gdp_agr <- as.numeric(gsub(",", ".", faktor$gdp_agr))
faktor$oil_price <- as.numeric(gsub(",", ".", faktor$oil_price))
faktor$coal_price <- as.numeric(gsub(",", ".", faktor$coal_price))
faktor$cpo_price <- as.numeric(gsub(",", ".", faktor$cpo_price))
str(faktor)
## tibble [180 × 17] (S3: tbl_df/tbl/data.frame)
## $ date : POSIXct[1:180], format: "2010-01-01" "2010-02-01" ...
## $ usd_idr : num [1:180] 9350 9337 9090 9012 9175 ...
## $ inflasion_mom: num [1:180] 0.84 0.3 -0.14 0.15 0.29 0.97 1.57 0.76 0.44 0.06 ...
## $ bi_rate : num [1:180] 6.5 6.5 6.5 6.5 6.5 6.5 6.5 6.5 6.5 6.5 ...
## $ fed_rate : num [1:180] 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ trade_balance: num [1:180] 2.11 1.67 1.8 0.8 2.64 0.57 -0.14 1.55 2.53 2.28 ...
## $ cad_devisa : num [1:180] 69.6 69.7 71.8 78.6 74.6 ...
## $ ekspor : num [1:180] 11.6 11.2 12.8 12 12.6 ...
## $ impor : num [1:180] 9.49 9.5 11 11.2 9.98 11.8 12.6 12.2 9.65 12.1 ...
## $ m1 : num [1:180] 490 494 495 514 545 ...
## $ m2 : num [1:180] 2066 2112 2116 2143 2231 ...
## $ ihsg : num [1:180] 2611 2549 2777 2971 2797 ...
## $ utang_luar : num [1:180] 181 181 181 183 183 183 196 196 196 202 ...
## $ gdp_agr : num [1:180] 5.69 5.69 5.69 6.17 6.17 6.17 5.8 5.8 5.8 6.9 ...
## $ oil_price : num [1:180] 72.9 79.7 83.8 86.2 74 ...
## $ cpo_price : num [1:180] 2445 2595 2556 2558 2436 ...
## $ coal_price : num [1:180] 77.4 87.8 86.6 86.6 92.1 ...
# Standarisasi data
faktor_baru <- faktor[, -1]
X <- as.matrix(faktor_baru[, -which(names(faktor_baru) == "usd_idr")])
y <- faktor_baru$usd_idr
X_standar <- scale(X)
# Pemodelan OLS
model_ols <- lm(y ~ ., data = as.data.frame(X_standar))
summary(model_ols)
##
## Call:
## lm(formula = y ~ ., data = as.data.frame(X_standar))
##
## Residuals:
## Min 1Q Median 3Q Max
## -946.85 -302.33 -32.85 273.64 1026.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12867.531 30.718 418.891 < 2e-16 ***
## inflasion_mom -61.452 33.782 -1.819 0.070724 .
## bi_rate 260.980 52.204 4.999 1.47e-06 ***
## fed_rate -290.895 98.866 -2.942 0.003730 **
## trade_balance 3.492 164.975 0.021 0.983138
## cad_devisa -567.537 92.114 -6.161 5.39e-09 ***
## ekspor -99.052 373.280 -0.265 0.791069
## impor 28.343 276.427 0.103 0.918460
## m1 -3412.363 720.874 -4.734 4.74e-06 ***
## m2 6082.443 854.483 7.118 3.25e-11 ***
## ihsg 139.915 132.584 1.055 0.292845
## utang_luar 211.436 260.540 0.812 0.418237
## gdp_agr 143.492 53.935 2.660 0.008579 **
## oil_price -249.217 64.408 -3.869 0.000157 ***
## cpo_price -206.510 84.959 -2.431 0.016147 *
## coal_price 125.795 61.337 2.051 0.041868 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 412.1 on 164 degrees of freedom
## Multiple R-squared: 0.9711, Adjusted R-squared: 0.9685
## F-statistic: 367.7 on 15 and 164 DF, p-value: < 2.2e-16
# Uji asumsi multikolinearitas model OLS
#install.packages("car")
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.3
vif(model_ols)
## inflasion_mom bi_rate fed_rate trade_balance cad_devisa
## 1.202731 2.872122 10.301107 28.683280 8.942238
## ekspor impor m1 m2 ihsg
## 146.845644 80.529178 547.659340 769.481452 18.525726
## utang_luar gdp_agr oil_price cpo_price coal_price
## 71.538579 3.065781 4.371946 7.606891 3.964893
# Packages regresi LASSO dan ridge
#install.packages("glmnet")
library(glmnet)
## Warning: package 'glmnet' was built under R version 4.3.3
## Loading required package: Matrix
## Loaded glmnet 4.1-8
# Function rolling window cross-validation
rolling_cv_detailed <- function(X, y, window_size = 36, alpha = 1, lambda_seq = NULL) {
n <- nrow(X)
rmse_vec <- c()
rmse_sd_vec <- c()
for (lambda in lambda_seq) {
errors <- c()
for (i in (window_size + 1):n) {
train_X <- X[(i - window_size):(i - 1), ]
train_y <- y[(i - window_size):(i - 1)]
test_X <- X[i, , drop = FALSE]
test_y <- y[i]
model <- glmnet(train_X, train_y, alpha = alpha, lambda = lambda)
pred <- predict(model, newx = test_X)
err <- (test_y - pred)^2
errors <- c(errors, err)
}
rmse_vec <- c(rmse_vec, mean(sqrt(errors)))
rmse_sd_vec <- c(rmse_sd_vec, sd(sqrt(errors)))
}
best_lambda <- lambda_seq[which.min(rmse_vec)]
list(
lambda_seq = lambda_seq,
rmse = rmse_vec,
rmse_sd = rmse_sd_vec,
lambda.min = best_lambda
)
}
# Pemodelan regresi LASSO
set.seed(123)
lambda_seq <- 10^seq(4, -4, length = 100)
lasso_result <- rolling_cv_detailed(X_standar, y, window_size = 36, alpha = 1, lambda_seq = lambda_seq)
best_lambda_lasso <- lambda_seq[which.min(lasso_result$rmse)]
lasso_model <- glmnet(X_standar, y, alpha = 1, lambda = best_lambda_lasso)
# Buat Plot
lambda_vals <- lasso_result$lambda_seq
log_lambda <- log(lambda_vals)
rmse_vals <- lasso_result$rmse
rmse_sd <- lasso_result$rmse_sd
cap_width <- 0.05
plot(log_lambda, rmse_vals, type = "p", pch = 20, col = "blue",
xlab = expression(Log(lambda)), ylab = "Root Mean-Squared Error",
ylim = c(min(rmse_vals - rmse_sd), max(rmse_vals + rmse_sd)), cex.axis = 0.9)
segments(log_lambda, rmse_vals - rmse_sd,
log_lambda, rmse_vals + rmse_sd, col = "grey")
segments(log_lambda - cap_width, rmse_vals + rmse_sd,
log_lambda + cap_width, rmse_vals + rmse_sd, col = "grey")
segments(log_lambda - cap_width, rmse_vals - rmse_sd,
log_lambda + cap_width, rmse_vals - rmse_sd, col = "grey")
abline(v = log(lasso_result$lambda.min), lty = 2)
tick_pos <- seq(floor(min(log_lambda)), ceiling(max(log_lambda)), by = 1)
axis(side = 1, at = tick_pos, labels = round(tick_pos, 1), cex.axis = 0.9)
model <- glmnet(X_standar, y, alpha = 1, lambda = lambda_vals, standardize = FALSE)
nzero <- model$df
axis(3, at = log_lambda, labels = nzero, las = 1, cex.axis = 0.7, hadj = 0.5)
# Pemodelan regresi ridge
set.seed(123)
lambda_seq <- 10^seq(4, -4, length = 100)
ridge_result <- rolling_cv_detailed(X_standar, y, window_size = 36, alpha = 0, lambda_seq = lambda_seq)
best_lambda_ridge <- lambda_seq[which.min(ridge_result$rmse)]
ridge_model <- glmnet(X_standar, y, alpha = 0, lambda = best_lambda_ridge)
# Buat Plot
lambda_vals <- ridge_result$lambda_seq
log_lambda <- log(lambda_vals)
rmse_vals <- ridge_result$rmse
rmse_sd <- ridge_result$rmse_sd
cap_width <- 0.05
plot(log_lambda, rmse_vals, type = "p", pch = 20, col = "blue",
xlab = expression(Log(lambda)), ylab = "Root Mean-Squared Error",
ylim = c(min(rmse_vals - rmse_sd), max(rmse_vals + rmse_sd)), cex.axis = 0.9)
segments(log_lambda, rmse_vals - rmse_sd,
log_lambda, rmse_vals + rmse_sd, col = "grey")
segments(log_lambda - cap_width, rmse_vals + rmse_sd,
log_lambda + cap_width, rmse_vals + rmse_sd, col = "grey")
segments(log_lambda - cap_width, rmse_vals - rmse_sd,
log_lambda + cap_width, rmse_vals - rmse_sd, col = "grey")
abline(v = log(ridge_result$lambda.min), lty = 2)
tick_pos <- seq(floor(min(log_lambda)), ceiling(max(log_lambda)), by = 1)
axis(side = 1, at = tick_pos, labels = round(tick_pos, 1), cex.axis = 0.9)
model <- glmnet(X_standar, y, alpha = 0, lambda = lambda_vals, standardize = FALSE)
nzero <- model$df
axis(3, at = log_lambda, labels = nzero, las = 1, cex.axis = 0.7, hadj = 0.5)
# Ekstrak koefisien
coef_ridge <- coef(ridge_model)
coef_lasso <- coef(lasso_model)
# Koefisien Lasso yang tidak nol (variabel terpilih)
selected_lasso <- coef_lasso[coef_lasso[, 1] != 0, ]
selected_lasso
## (Intercept) inflasion_mom bi_rate fed_rate cad_devisa
## 12867.53056 -78.26657 316.00950 -346.08018 -611.11340
## ekspor impor m2 ihsg utang_luar
## -93.57950 -11.07751 1906.82704 319.89165 1015.17464
## gdp_agr oil_price cpo_price coal_price
## 136.79641 -355.46692 -215.26206 72.73320
# Koefisien ridge
coef_ridge
## 16 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 12867.530556
## inflasion_mom -86.437259
## bi_rate 329.092233
## fed_rate -332.518845
## trade_balance 9.853282
## cad_devisa -585.252069
## ekspor -117.881332
## impor -11.090091
## m1 438.539377
## m2 1248.650366
## ihsg 432.469229
## utang_luar 1078.096102
## gdp_agr 113.956219
## oil_price -393.905853
## cpo_price -194.082827
## coal_price 89.552570
library(Metrics)
## Warning: package 'Metrics' was built under R version 4.3.3
# Prediksi
pred_ridge <- predict(ridge_model, s = best_lambda_ridge, newx = X_standar)
pred_lasso <- predict(lasso_model, s = best_lambda_lasso, newx = X_standar)
# Metrik evaluasi
rmse_ridge <- rmse(y, pred_ridge)
rmse_lasso <- rmse(y, pred_lasso)
mape_ridge_persen <- mape(y, pred_ridge)*100
mape_lasso_persen <- mape(y, pred_lasso)*100
r2_ridge <- 1 - sum((y - pred_ridge)^2) / sum((y - mean(y))^2)
r2_lasso <- 1 - sum((y - pred_lasso)^2) / sum((y - mean(y))^2)
# Tampilkan hasil
data.frame(
Model = c("Ridge", "Lasso"),
RMSE = c(rmse_ridge, rmse_lasso),
MAPE = c(mape_ridge_persen, mape_lasso_persen),
R2 = c(r2_ridge, r2_lasso)
)