Econometrics

Finals Exam


Kontak : \(\downarrow\)
Email
Instagram https://www.instagram.com/arifin.alicia/
RPubs https://rpubs.com/aliciaarifin/
Nama Alicia Arifin
NIM 20214920001
Prodi Statistika, 2021

Data For Sales Analysis

The dataset below represents the monthly sales data of a company, including various factors that might influence product sales. The data is presented in a tabular format with each row representing a month and each column representing a variable

library(tibble)
sales_data <- tibble::tibble(
  Month = seq.Date(from = as.Date("2019-01-01"), to = as.Date("2024-05-01"), by = "1 month"),
Advertising_Expense = c(50, 75, 60, 65, 80, 85, 90, 95, 100, 105,110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170, 175, 180, 185, 190, 195, 200, 205, 210, 215, 220,225, 230, 235, 240, 245, 250, 255, 260, 265, 270, 275, 280, 285, 290, 295, 300, 305, 310, 315, 320, 325, 330, 335,340, 345, 350, 355, 360, 365, 370, 375, 380),
Product_Quality = c(20, 25, 22, 23, 28, 30, 32, 35, 38, 40, 42, 45, 48, 50, 55, 60, 65, 68, 70, 72, 75, 78, 80, 82, 85, 88,90, 92, 95, 98, 100, 102, 105, 108, 110, 112, 115, 118, 120, 122, 125, 128, 130, 132, 135, 138, 140, 142, 145, 148,150, 152, 155, 158, 160, 162, 165, 168, 170, 172, 175, 178, 180, 182, 185),
Product_Price = c(95, 112, 107,100, 110, 105, 108, 115, 120, 125,130, 135, 140, 145, 150, 155, 160, 165, 170, 175, 180, 185, 190, 195, 200, 205, 210, 215, 220, 225, 230, 235, 240,245, 250, 255, 260, 265, 270, 275, 280, 285, 290, 295, 300, 305, 310, 315, 320, 325, 330, 335, 340, 345, 350, 355,360, 365, 370, 375, 380, 385, 390, 395, 400),
Sales_Promotion = c(15, 20, 18, 19, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64,66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110, 112, 114, 116, 118, 120, 122, 124, 126, 128, 130, 132, 134, 136, 138, 140,142),
Online_Marketing = c(76, 77, 78, 79, 80, 85, 82, 84, 88, 90, 92,94, 96, 98, 100, 102, 104, 106, 108, 110, 112, 114, 116, 118, 120, 122, 124, 126, 128, 130, 132, 134, 136, 138, 140,142, 144, 146, 148, 150, 152, 154, 156, 158, 160, 162, 164, 166, 168, 170, 172, 174, 176, 178, 180, 182, 184, 186,188, 190, 192, 194, 196, 198, 200),
Offline_Marketing = c(84,85, 87, 89, 89, 90, 95, 92, 94, 98, 100,102, 104, 106, 108, 110, 112, 114, 116, 118, 120, 122, 124, 126, 128, 130, 132, 134, 136, 138, 140, 142, 144, 146,148, 150, 152, 154, 156, 158, 160, 162, 164, 166, 168, 170, 172, 174, 176, 178, 180, 182, 184, 186, 188, 190, 192,194, 196, 198, 200, 202, 204, 206, 208),
Product_Sales = NA)

num = c(0:58)
num_1 = 250 +num*5
num = c(200,220,210,215,230,240)
tn <- c(num,num_1)
sales_data$Product_Sales = tn

library(data.table)
data.table(sales_data)

This dataset includes the following variables:

  • Month: The month of observation.

  • Advertising_Expense: The amount spent on advertising (in dollars).

  • Product_Quality: The quality score of the product (on a scale from 1 to 100).

  • Product_Price: The price of the product (in dollars).

  • Sales_Promotion: The expenditure on sales promotion (in dollars).

  • Online_Marketing: The expenditure on online marketing (in dollars).

  • Offline_Marketing: The expenditure on offline marketing (in dollars).

  • Product_Sales: The number of products sold

Regression Analysis

Pada analisis regresi kita akan cek, semua variabel untuk memprediksi Product Sales.

lm1 <- lm(Product_Sales~., data=sales_data)
summary(lm1)
## 
## Call:
## lm(formula = Product_Sales ~ ., data = sales_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.91093 -0.07192  0.00035  0.06233  2.77371 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -430.21314  345.90017  -1.244  0.21868    
## Month                  0.04153    0.01971   2.107  0.03956 *  
## Advertising_Expense    0.51579    0.09814   5.256 2.29e-06 ***
## Product_Quality        0.12570    0.04469   2.812  0.00673 ** 
## Product_Price         -0.73378    0.05896 -12.445  < 2e-16 ***
## Sales_Promotion        3.88609    0.62879   6.180 7.28e-08 ***
## Online_Marketing      -0.92713    0.17745  -5.225 2.57e-06 ***
## Offline_Marketing     -0.70134    0.16460  -4.261 7.72e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7266 on 57 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 1.616e+05 on 7 and 57 DF,  p-value: < 2.2e-16

Dari hasil regresi di atas, didapatkan bahwa r square dari model senilai 0.9999. Semua variabel memiliki dampak yang siknifikan terhadap product sales. perlu disayangkan bahwa intercept dari model tidak siknifikan. Lalu, mari kita coba jika kita tidak melibatkan variabel waktu (month).

reg_data <- sales_data[,-1]
lm2 <- lm(Product_Sales~., data= reg_data) 
summary(lm2)
## 
## Call:
## lm(formula = Product_Sales ~ ., data = reg_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6031 -0.0697 -0.0051  0.0594  3.3179 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         297.58588   17.69853  16.814  < 2e-16 ***
## Advertising_Expense   0.46830    0.09830   4.764 1.31e-05 ***
## Product_Quality       0.12252    0.04597   2.665 0.009955 ** 
## Product_Price        -0.82782    0.03965 -20.880  < 2e-16 ***
## Sales_Promotion       4.68575    0.51594   9.082 9.71e-13 ***
## Online_Marketing     -0.85343    0.17905  -4.767 1.30e-05 ***
## Offline_Marketing    -0.58501    0.15959  -3.666 0.000537 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7478 on 58 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 1.78e+05 on 6 and 58 DF,  p-value: < 2.2e-16

Dari hasil regresi dengan tidak menggunakan variabel waktu, semua variabel independen memiliki pengaruh siknifikan terhadap product sales. R-square yang dihasilkan juga hampir sempurna. Mari gunakan model ini untuk uji asumsi di bawah.

Assumption

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

Residual vs Fitted

Plot ini digunakan untuk mendeteksi apakah data memiliki sifat linearitas dan kehomogenan. Kehomogenan suatu data tidak akan berbentuk pola-pola tertentu. semakin tersebar dengan baik, maka semakin homogen. Homogen adalah varians pada model konstan atau tidak. Regresi linear membutuhkan varians yang konstan agar mengurangi error. Linearitas bisa didekeksi dengan melihat persebaran titik dengan residunya. Data akan semakin linear jika titik bersebar mendekati 0.
Jika dilihat dari Plot, dikatakan bahwa di titik-titik awal (titik 1-20an) yang bersebar lalu berbentuk di laris lurus diantara titik 0 residuals. Data bisa disebut linear. kehomogenitasnya bisa diragukan karena cukup banyak data yang beda sabaran di range fitted values 200-300.

Q-Q Residuals

Plot ini digunakan untuk mendeteksi data berdistribusi normal dari data.

shapiro.test(rstandard(lm2))
## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(lm2)
## W = 0.64578, p-value = 2.993e-11

dari hasil tersebut pvalue < \(\alpha\), dikatakan bahwa data tidak berdistribusi dengan normal

hist(rstandard(lm2))

dari grafik histogram juga dilihat bahwa tidak berdistribusi dnegan normal. Data skewed ke kanan.

Scale-Location

Breusch-Pagan Test Hypothesis:

  • Null Hypothesis (H0): The residuals are homoscedastic (i.e. evenly spread)

  • Alternative Hypothesis (HA): The residuals are heteroscedastic (i.e. not evenly spread)

library(lmtest)
bptest(lm2)
## 
##  studentized Breusch-Pagan test
## 
## data:  lm2
## BP = 40.3, df = 6, p-value = 3.977e-07

dari hasil plot scale location, kita bisa melihat apakah data memiliki varians yang tetap atau tidak. persebaran bagusnya berada di antara satu garis horizontal yang lurus. bisa dilihat bahwa garis merahnya seperti melandai terlebih dahulu, lalu konstan. Dari hasil itu kita bisa ketahui bahwa data tidak homogen, didukung dengan Breusch-Pagan Test (p-vale <\(\alpha\)).
Jika dilihat dari plot, dan dilihat banyak yang anggapannya variansnya konstan, kita akan lanjut saja untuk membuat model nya.

Standarized Residuals

Plot ini digunakan untuk melihat apakah adanya outlier. Outlier dideteksi ada 3. Maka kita akan cari dimana outlier itu berada.

outliers <- which(abs(resid(lm2)) > 2 * sd(resid(lm2)))
outliers
## 3 4 8 
## 3 4 8

terdapat 3 outlier/pencilan yaitu data ke 3,4 dan 8. nanti akan dihilangkan pencilan.

Independence Test

library(lmtest)
dwtest(lm2)
## 
##  Durbin-Watson test
## 
## data:  lm2
## DW = 2.2033, p-value = 0.6809
## alternative hypothesis: true autocorrelation is greater than 0

Dari hasil Durbin-Watson Statistic, didapatkan bahwa p-value > \(\alpha\) 5%. H0 di terima, tidak ada autokorelasi antar residual.

Multicollinearity

library(car)
m_data <- sales_data[,c(2:8)]
cor(m_data)
##                     Advertising_Expense Product_Quality Product_Price
## Advertising_Expense           1.0000000       0.9990361     0.9986571
## Product_Quality               0.9990361       1.0000000     0.9980254
## Product_Price                 0.9986571       0.9980254     1.0000000
## Sales_Promotion               0.9997737       0.9990737     0.9993783
## Online_Marketing              0.9993106       0.9988327     0.9994418
## Offline_Marketing             0.9991459       0.9986337     0.9994426
## Product_Sales                 0.9996214       0.9987346     0.9970366
##                     Sales_Promotion Online_Marketing Offline_Marketing
## Advertising_Expense       0.9997737        0.9993106         0.9991459
## Product_Quality           0.9990737        0.9988327         0.9986337
## Product_Price             0.9993783        0.9994418         0.9994426
## Sales_Promotion           1.0000000        0.9998151         0.9997317
## Online_Marketing          0.9998151        1.0000000         0.9997606
## Offline_Marketing         0.9997317        0.9997606         1.0000000
## Product_Sales             0.9990594        0.9983740         0.9981576
##                     Product_Sales
## Advertising_Expense     0.9996214
## Product_Quality         0.9987346
## Product_Price           0.9970366
## Sales_Promotion         0.9990594
## Online_Marketing        0.9983740
## Offline_Marketing       0.9981576
## Product_Sales           1.0000000
vif(lm2)
## Advertising_Expense     Product_Quality       Product_Price     Sales_Promotion 
##           9994.0793            592.2988           1549.4610          43328.2612 
##    Online_Marketing   Offline_Marketing 
##           5182.5024           4096.6433

Dari hasil VIF didapatkan bahwa semua variabel independen mempunyai nilai VIF lebih dari 10 dan nilai korelasinya hampir 100%. Unutk menyelesaikan masalah multikolinearitas, kita akan membuat regresi ridge. Lalu bagaimana ketika ingin membuat regresi linear biasa? Kita akan menggunakan variabel yang memiliki nilai korelasi yang terbesar dengan Product_Sales yaitu variabel Advertising_Expense sebesar 99.96%.

Regresi Linear

Regresi Linear akan menggunakan Advertising_Expense untuk model satu variabel.

mo_data = m_data[-c(3,4,8),]
lm3 <- lm(Product_Sales~Advertising_Expense, data= mo_data)
summary(lm3)
## 
## Call:
## lm(formula = Product_Sales ~ Advertising_Expense, data = mo_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.5046  -0.5370   0.3925   1.3220   2.3125 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.566e+02  8.239e-01   190.1   <2e-16 ***
## Advertising_Expense 1.012e+00  3.373e-03   300.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.41 on 60 degrees of freedom
## Multiple R-squared:  0.9993, Adjusted R-squared:  0.9993 
## F-statistic: 9.005e+04 on 1 and 60 DF,  p-value: < 2.2e-16

Model regresi linear yang telah dibuat adalah
\(Product Sales = 1.566e+02 + 1.012e+00 Advertising Expense\)
Dari hasil summary model, didapatkan p-value model < \(\alpha\) yang berarti Advertising_Expense memiliki pengaruh siknifikan terhadap Product_Sales. Nilai dari R-squarednya adlah 0.9993, yang berarti variabel Advertising_Expense dapat menjelaskan atau sudah masuk ke model bebanyak 99.93%.
coba kita diagnostic plotnya

plot(sales_data$Advertising_Expense,sales_data$Product_Sales)

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

Asumsi linearitas terpenuhi, Dari resiudal vs fitted dan scale-location juga bisa dibilang hampir homogen (masih tidak di satu garis lurus). Dari QQ-Residuals bisa dibilang linear dan tidak terlalu banyak skewednya. dan memiliki outlier 3 titik pada plot Residuals vs leverage.

Regresi Ridge

Pada data, didapatkan bahwa semua variabel memiliki multikolinearitas. Untuk itu, kita akan hilangkan multikolinearitasnya.

# Load libraries, get data & set seed for reproducibility ---------------------
set.seed(1993)    # seef for reproducibility
library(glmnet)  # for ridge regression
library(dplyr)   # for data cleaning
library(psych)   # for function tr() to compute trace of a matrix

y <- mo_data %>% select(Product_Sales) %>% scale(center = TRUE, scale = FALSE) %>% as.matrix()
X <- mo_data %>% select(-c(Product_Sales)) %>% as.matrix()


# Perform 10-fold cross-validation to select lambda ---------------------------
lambdas_to_try <- 10^seq(-3, 5, length.out = 100)
# Setting alpha = 0 implements ridge regression
ridge_cv <- cv.glmnet(X, y, alpha = 0, lambda = lambdas_to_try,
                      standardize = TRUE, nfolds = 10)
# Plot cross-validation results
plot(ridge_cv)

# Best cross-validated lambda
lambda_cv <- ridge_cv$lambda.min
# Fit final model, get its sum of squared residuals and multiple R-squared
model_cv <- glmnet(X, y, alpha = 0, lambda = lambda_cv)
y_hat_cv <- predict(model_cv, X)
ssr_cv <- t(y - y_hat_cv) %*% (y - y_hat_cv)
rsq_ridge_cv <- cor(y, y_hat_cv)^2



# Use information criteria to select lambda -----------------------------------
X_scaled <- scale(X)
aic <- c()
bic <- c()



for (lambda in seq(lambdas_to_try)) {
  # Run model
  model <- glmnet(X, y, alpha = 0, lambda = lambdas_to_try[lambda], standardize = TRUE)
  # Extract coefficients and residuals (remove first row for the intercept)
  betas <- as.vector((as.matrix(coef(model))[-1, ]))
  resid <- y - (X_scaled %*% betas)
  # Compute hat-matrix and degrees of freedom
  ld <- lambdas_to_try[lambda] * diag(ncol(X_scaled))
  H <- X_scaled %*% solve(t(X_scaled) %*% X_scaled + ld) %*% t(X_scaled)
  df <- tr(H)
  # Compute information criteria
  aic[lambda] <- nrow(X_scaled) * log(t(resid) %*% resid) + 2 * df
  bic[lambda] <- nrow(X_scaled) * log(t(resid) %*% resid) + 2 * df * log(nrow(X_scaled))
}
# Optimal lambdas according to both criteria
lambda_aic <- lambdas_to_try[which.min(aic)]
lambda_bic <- lambdas_to_try[which.min(bic)]
lambda_aic==lambda_bic # melihat apakah lamda aic dan bic yang terkecil sama atau tidak
## [1] FALSE

lamdanya berbeda sehingga kita akan membuat 2 model dengan lamda aic dan lambda bic.

ridge_aic <- lmridge::lmridge(Product_Sales~., data = mo_data, lambda = lambda_aic)
summary(ridge_aic)
## 
## Call:
## lmridge.default(formula = Product_Sales ~ ., data = mo_data, 
##     lambda = lambda_aic)
## 
## 
## Coefficients: for Ridge parameter K= 0 
##                       Estimate Estimate (Sc) StdErr (Sc) t-value (Sc) Pr(>|t|)
## Intercept             339.2721    48839.7114  13010.6916       3.7538   0.0004
## Advertising_Expense     0.3517      251.2129     38.3348       6.5531   <2e-16
## Product_Quality         0.0701       26.0468      7.7709       3.3518   0.0015
## Product_Price          -0.9345     -657.3995     13.4417     -48.9075   <2e-16
## Sales_Promotion         5.8533     1665.5297     88.8824      18.7386   <2e-16
## Online_Marketing       -1.1650     -330.7348     27.5838     -11.9902   <2e-16
## Offline_Marketing      -0.8183     -232.1433     26.9862      -8.6023   <2e-16
##                        
## Intercept           ***
## Advertising_Expense ***
## Product_Quality     ** 
## Product_Price       ***
## Sales_Promotion     ***
## Online_Marketing    ***
## Offline_Marketing   ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Ridge Summary
##           R2       adj-R2     DF ridge            F          AIC          BIC 
##      1.00000      1.00000      5.99996 817367.25422   -133.06349    135.58156 
## Ridge minimum MSE= 11099.81 at K= 0 
## P-value for F-test ( 5.99996 , 56.00007 ) = 1.784767e-136 
## -------------------------------------------------------------------
rsq_ridge_aic <- cor(mo_data$Product_Sales, predict(ridge_aic))^2 ; rsq_ridge_aic 
## [1] 0.9999886
ridge_bic <- lmridge::lmridge(Product_Sales~., data = mo_data, lambda = lambda_bic)
summary(ridge_bic)
## 
## Call:
## lmridge.default(formula = Product_Sales ~ ., data = mo_data, 
##     lambda = lambda_bic)
## 
## 
## Coefficients: for Ridge parameter K= 0 
##                       Estimate Estimate (Sc) StdErr (Sc) t-value (Sc) Pr(>|t|)
## Intercept             339.2721    48839.7114  13010.6916       3.7538   0.0004
## Advertising_Expense     0.3517      251.2129     38.3348       6.5531   <2e-16
## Product_Quality         0.0701       26.0468      7.7709       3.3518   0.0015
## Product_Price          -0.9345     -657.3995     13.4417     -48.9075   <2e-16
## Sales_Promotion         5.8533     1665.5297     88.8824      18.7386   <2e-16
## Online_Marketing       -1.1650     -330.7348     27.5838     -11.9902   <2e-16
## Offline_Marketing      -0.8183     -232.1433     26.9862      -8.6023   <2e-16
##                        
## Intercept           ***
## Advertising_Expense ***
## Product_Quality     ** 
## Product_Price       ***
## Sales_Promotion     ***
## Online_Marketing    ***
## Offline_Marketing   ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Ridge Summary
##           R2       adj-R2     DF ridge            F          AIC          BIC 
##      1.00000      1.00000      5.99996 817367.25422   -133.06349    135.58156 
## Ridge minimum MSE= 11099.81 at K= 0 
## P-value for F-test ( 5.99996 , 56.00007 ) = 1.784767e-136 
## -------------------------------------------------------------------
rsq_ridge_bic <- cor(mo_data$Product_Sales, predict(ridge_bic))^2 ; rsq_ridge_bic
## [1] 0.9999886

Jika dilihat, regresi ridge aic dan bic memiliki hasil nilai summary yang sama hingga koefisien dan rsquarenya. maka dari itu kita akn pilih salah satu saja karena sama modelnya.


Implementasi Regresi RIdge :
Pada regresi ridge, nilai p-valuenya kurang dari \(\alpha\), yang berarti model memiliki setidaknya satu variabel independen yang siknifikan. Semua variabel independen siknifikan terhadap variabel dependen. Nilai R-squarenya 0.9999886 atau bisa dibulatkan 1. Model hampir sempurna. 


coef(ridge_aic)
##           Intercept Advertising_Expense     Product_Quality       Product_Price 
##           339.27206             0.35166             0.07008            -0.93449 
##     Sales_Promotion    Online_Marketing   Offline_Marketing 
##             5.85329            -1.16500            -0.81833

Persamaan regresinya adalah
Product Sales = 339.27206 + 0.35166 Advertising_Expense + 0.07008 Product_Quality - 0.93449 Product_Price + 5.85329 Sales_Promotion - 1.16500 Online_Marketing - 0.81833 Offline_Marketing

Regression Model Selection

r_error <- data.frame(
  actual = mo_data$Product_Sales
)

r_error$lm <- abs(predict(lm3) - mo_data$Product_Sales) / mo_data$Product_Sales *100
r_error$ridge=  abs(predict(ridge_aic) - mo_data$Product_Sales) / mo_data$Product_Sales *100


data.table(r_error)
mean(r_error$lm)
## [1] 0.484986
mean(r_error$ridge)
## [1] 0.05324759

Berikut hasil nilai error dalam bentuk Percentage Error (%).

Karena Regresi linear

library(highcharter)
hchart(sales_data[-c(3,4,8),], type = "line", 
         hcaes(x= Advertising_Expense, y= Product_Sales), name="Sales Product by Month")%>%
  hc_add_series(sales_data[-c(3,4,8),], type="line", hcaes(x=Advertising_Expense, y= predict(lm3)), name="Linear Regression")%>% 
  hc_add_series(sales_data[-c(3,4,8),], type="line", hcaes(x=Advertising_Expense, y= predict(ridge_aic)),name="Ridge Regression")%>% 
  hc_yAxis(opposite = TRUE,
         labels = list(format = "{value}"))%>%
  hc_title(    text = "Product Sales 2019-2024",
    margin = 20,
    align = "left",
    style = list(color = "#240750", useHTML = TRUE))

Jika dilihat dari model, dapat ditentukan bahwa kedua regresi cocok dengan memiliki nilai MAPE kurang dari 1%.

Time Series Analysis

Plot Time Series

sales_data %>%
  hchart(., type = "line", 
         hcaes(x= Month, y= Product_Sales), name="Sales Product by Month")%>%
  hc_yAxis(opposite = TRUE,
         labels = list(format = "{value}"))%>%
  hc_title(    text = "Product Sales 2019-2024",
    margin = 20,
    align = "left",
    style = list(color = "#240750", useHTML = TRUE))

Possible Time-Series Modeling

Model yang akan coba kita bandingkan adalah Holt’s Linear, dan metode BOx-Jenkins.

Stasionerity

Kita perlu menguji stasioner karena metode Box-Jenkins membutuhkan data yang stasioner, stasioner terhadap rata-rata maupun terhadap varians. Dari hasil plot time-series yang telah ada, dapat dilihat bahwa plot menunjukkan adanya sifat linear. Dari bulan ke bulan, penjualan produk semakin naik (trend naik). Jika dilihat dari Musiman atau siklusnya, tidak ada Siklus dan tidak ada musiman. Untuk mengetahui lebih lanjut mari kita uji asumsi terlebih dahulu.
Pada subbab ini, kita akan menguji kestasioneritas data. Uji stasioner dibagi menjadi 2 yaitu stasioner terhadap rata-rata dan stasioner terhadap varians. Cara untuk mengstasionerkan data yang tidak stasioner terhadap rata-rata adalah transformasi. Cara untuk mengstasionerkan data yang tidak stasioner terhadap varians adalah differencing.


Pada grafik ACF kita bisa menentukan apakah terdapat seasonal (musiman) atau stasioner terhadap varians. Ciri-ciri data seasonal adalah ketika grafik ACF plotnya akan membentuk suatu pola yang berulang dan jeda antar waktunya sama dengan pola sebelumnya. Data seasonal dan data siklus memiliki perbedaan sedikit, Data siklus juga akan membentuk pola yang berulang, tetapi jeda waktu yang diperlukan berbeda. Grafik ACF juga bisa menentukan data telah stasioner terhadap varians atau tidak dilihat apakah nilai-nilai korelasi lag-nya ada di antara critical value atau tidak. 95% diantara critical value = stasioner terhadap varians (white noise).
Box-Jenkins method mempunyai notasi ARIMA(\(p,d,q\)). Banyaknya jumlah data didifferencing akan menentukan berapa nilai \(d\) pada I(\(d\)). Jumlah nilai yang lebih dari batas critical value digunakan untuk mencari berapa AR(p) dan MA(q). Cara menentukan p dan q sesuai dengan tabel di bawah ini :

Tabel 1. AR dan MA di plot ACF dan PACF Cara mudah untuk menentukannya adalah dengan menghitung berapa lonjakan yang di luar batasan critical value.
Mari kita langsung coba cek plot ACF-nya :

library(tseries) 
ts_test =sales_data$Product_Sales
print(acf(ts_test,lag=65,pl=TRUE)) 

## 
## Autocorrelations of series 'ts_test', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.949  0.904  0.853  0.803  0.756  0.710  0.666  0.622  0.579  0.536 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.493  0.451  0.410  0.369  0.328  0.288  0.249  0.211  0.173  0.136  0.100 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.065  0.030 -0.003 -0.035 -0.067 -0.097 -0.126 -0.154 -0.180 -0.206 -0.230 
##     33     34     35     36     37     38     39     40     41     42     43 
## -0.253 -0.274 -0.294 -0.313 -0.330 -0.345 -0.359 -0.371 -0.382 -0.391 -0.398 
##     44     45     46     47     48     49     50     51     52     53     54 
## -0.403 -0.406 -0.408 -0.407 -0.405 -0.401 -0.394 -0.386 -0.375 -0.362 -0.347 
##     55     56     57     58     59     60     61     62     63     64 
## -0.330 -0.311 -0.289 -0.264 -0.238 -0.207 -0.172 -0.132 -0.090 -0.048

Jika dilihat dari plot ACF, dapat lihat bahwa tidak ada pola berulang yang terbentuk, maka data tidak memiliki musiman. korelasi Lag-nya banyak yang melebihi batas critical value, maka data tidak stasioner terhadap varians. Maka dari itu, kita perlu lalukan differencing dan plot grafiknya ulang.
Differencing di R bisa menggunakan coding diff().

acf(diff(ts_test),lag=65,pl=TRUE)

Pada grafik ACF, terdapat 2 lonjakan pada lag 0 dan 1. sehingga AR yang cocok adalah 0-2. Grafik di atas membentuk model white noise (ketika 95% dari lonjakan berada di dalam critical value).

pacf(diff(ts_test),lag=65,pl=TRUE)

Pada grafik PACF, terdapat 2 lonjakan pada lag 1 dan 2. Sehingga MA yang cocok adalah antara 0-2. Grafik juga membentuk model white noise.

Model Building

Pembentukan model

Holt’s Linear Method

Holt’s Linear adalah metode lanjutan dari single exponential smoothing yang bisa digunakan untuk memprediksi data dengan trend. Metode ini menggunakan 2 parameter, yaitu \(\alpha\)(parameter pemulusan) dan \(\beta\) (parameter trend).

library(forecast)
HL_model <- holt(sales_data$Product_Sales, h=2)
summary(HL_model)
## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = sales_data$Product_Sales, h = 2) 
## 
##   Smoothing parameters:
##     alpha = 0.573 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 200.6437 
##     b = 5.2411 
## 
##   sigma:  2.9234
## 
##      AIC     AICc      BIC 
## 416.6651 417.6821 427.5371 
## 
## Error measures:
##                       ME     RMSE     MAE        MPE      MAPE      MASE
## Training set -0.03212918 2.832047 1.29566 0.02839918 0.5085973 0.2303395
##                     ACF1
## Training set -0.02768528
## 
## Forecasts:
##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 66       545.4205 541.6740 549.1670 539.6907 551.1503
## 67       550.6614 546.3432 554.9797 544.0573 557.2656

Dari hasil didapatkan \(\alpha\) = 0,573 dan \(\beta\) = 0,0001. Dapat dilihat bahwa nilai \(\beta\) sangat kecil yang berarti slope dari model berganti sangatlah tidak mungkin.

Dari hasil prediksi Holt’s Linear, untuk bulan berikutnya \(F_{t+1}\) adalah 545.42 dan forecast 2 bulan setelah itu \(F_{t+2}\) adalah 550.66.

Box Jenkins Method

Dari hasil grafik ACF dan PACF, didapatkan hasil model yang mungkin :

  • AR : 0-2

  • I : 1

  • MA : 0-2


# menentukan model arima terbaik dengan auto.arima()
arima_model <-auto.arima(sales_data$Product_Sales, max.p = 2, d=1, max.q = 2)
summary(arima_model)
## Series: sales_data$Product_Sales 
## ARIMA(2,1,2) with drift 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2   drift
##       -1.3317  -0.9433  1.8308  0.8406  5.2913
## s.e.   0.1659   0.0541  0.2004  0.2023  0.2992
## 
## sigma^2 = 4.932:  log likelihood = -144.16
## AIC=300.32   AICc=301.79   BIC=313.27
## 
## Training set error measures:
##                       ME     RMSE     MAE        MPE      MAPE      MASE
## Training set -0.01087948 2.115745 1.02334 0.02892078 0.3669924 0.1819271
##                    ACF1
## Training set -0.3067911

Model Selection

Pada subbab ini, kita akan menentukan model dengan menggunakan MAPE atau Mean Absolute Square Error. MAPE merupakan rata-rata dari nilai absolute dari persentase error yang dihasilkan. Persentase error didapatkan dari nilai forecast dikurangi nilai sebenarnya dan dibagi dengan nilai sebenarnya.

error_actual <- data.frame(
  time <- sales_data$Month,
  actual = sales_data$Product_Sales
)
error_actual$HoltsLinear <- abs(HL_model$residuals/ error_actual$actual)*100
error_actual$ARIMA <- abs(arima_model$residuals/error_actual$actual)*100
data.table(error_actual) # nilai percentage error time series
data.table(r_error) # nilai percentage error regression

untuk menentukan MAPE-nya akan menggunakan split data 80-20. 80 % initial data dan 20% tes data.

p20 <- nrow(error_actual)- nrow(error_actual)*0.2

MAPE <- data.frame(
  "Method" <- c("Holt's Linear", "ARIMA(2,1,2)"),
  "MAPE" <- c(mean(error_actual$HoltsLinear[p20:nrow(error_actual)]), mean(error_actual$ARIMA[p20:nrow(error_actual)]))
)
colnames(MAPE) <-c("Method","MAPE(%)")

Menghitung MAPE regression

MAPE_r <- data.frame(
  Method = c("Simple Lm", "Ridge Regrssion"),
  MAPE = c(mean(r_error$lm) , mean(r_error$ridge))
)
colnames(MAPE_r) <-c("Method","MAPE(%)")
MAPE.R <- data.table(rbind(MAPE,MAPE_r))

MAPE.R$"Rsq(%)" =c(
cor(error_actual$actual, fitted.values(HL_model))^2,
cor(error_actual$actual, fitted.values(arima_model))^2,
cor(error_actual$actual[-c(3,4,8)], fitted.values(lm3))^2,
cor(error_actual$actual[-c(3,4,8)], fitted.values(ridge_aic))^2 )


data.table(MAPE.R)

Jika dilihat dari hasil MAPE, ARIMA(2,1,2) dan regresi ridge memiliki nilai mape yang perbedaannya 0.001%. Karena terdapat multikolinearitas, semua variabel digunakan dan dimasukkan ke dalam regresi ridge. Regresi ridge memiliki r-square sebesar 0.9999886. Nilai R-square dari semua model lebih dari 0.99, yang berarti tidak sampai 1% variansi yang hilang dalam model. Nilai Rsq

Jika dilihat dari errornya, errornya sangat kecil untuk kedua model. bisa dikatakan model ARIMA(2,1,2) dan Regression Ridge merupakan model yang cocok.


berikut hasil fitted dengan model aslinya

t <- ts((sales_data$Product_Sales), start = c(2019,1), frequency = 12) 
t2 <- ts(c(fitted.values(arima_model)), start = c(2019,1), 
                frequency = 12) 
t3 <- ts(c(predict(ridge_aic) ), start = c(2019,1), 
                frequency = 12) 
ridge.predict <- predict(ridge_aic)
r.p <- c(ridge.predict[1:2],NA,NA, ridge.predict[3:6], NA, ridge.predict[7:62])

r.p<- ts(r.p, start=c(2019,1), frequency=12)

hchart(t, type = "line", name="Actual")%>%
  hc_add_series(t2, type="line", name="ARIMA")%>%
  hc_add_series(r.p, type="line", name= "ridge regression")%>%
  hc_yAxis(opposite = TRUE,
         labels = list(format = "{value}"))%>%
  hc_title(    text = "Product Sales 2019-2024",
    margin = 20,
    align = "left",
    style = list(color = "#240750", useHTML = TRUE))

ARIMA dan Regresi RIdge memiliki nilai error yang berbeda sangat sedikit 0.01%, Pada prediksi regressi ridge tidak agak sulit untuk dipresentasikan. Regresi ridge memiliki oulier di data ke 3,4 dan 8. sehingga sebenarnya pada regresi ridge tidak bisa digunakan untuk memvisualisasikan data time series. Dari kedua hasil regresi dan ARIMA, kedua model sama-sama memiliki nilai MAPE (ts model) dan R-square (regression model) yang sangat baik.

Study Case

Collect dataset from BPS to analyze the relationship between monthly economic indicators, including inflation rate, GDP growth rate, unemployment rate, interest rate, consumer confidence, stock market index, and exchange rate, over a specified period. In this study case you need to apply econometric techniques to analyze the relationships between various economic indicators and gain insights into the dynamics of the economy. By conducting correlation analysis, time series analysis, and regression modeling, you can identify key factors driving economic trends and inform decision-making processes.

Data Collection

Obtain the extended dataset containing monthly observations of economic indicators from January 2010 to December 2021.

library(readxl)
setwd(getwd())

interest_rate <- read_excel("data_UAS/INTEREST_RATE.xlsx")
stock_index <- read.csv("data_uas/JKLQ45_uas.csv")
consumer_confidence <- read_excel("data_UAS/Book1.xlsx", sheet= "kepercayaan_konsumen")
inflation_rate <- read_excel("data_UAS/Book1.xlsx", sheet= "tingkat_inflasi")
Kurs_USD <- read_excel("data_UAS/Kurs_Transaksi_USD.xlsx")

Data yang telah diambil merupakan data yang berasal dari BPS, BI dan sumber internet.

see data

data.table(interest_rate)
data.table(stock_index)
data.table(consumer_confidence)
data.table(inflation_rate)
data.table(Kurs_USD)

Rebuild Data

library(reshape2)
library(dplyr)
interest_rate <- melt(interest_rate, "Bulan/Tahun") 
interest_rate <- interest_rate%>%
  arrange(`Bulan/Tahun`)


consumer_confidence <- melt(consumer_confidence, "kepercayaan konsumen")
inflation_rate <- melt(inflation_rate, "tingkat inflasi")

dc <- data.frame(
  dates = seq(as.Date("2010/1/1"),as.Date("2021/12/1"), "month"),
  interest_rate = interest_rate$value,
  stock_index = stock_index$Close,
  consumer_confidence = consumer_confidence$value,
  inflation_rate = inflation_rate$value
)
library(dplyr)
library(tidyr)
USD_buy <- Kurs_USD%>%
  separate(Tanggal, c("bulan","hari", "tahun"))%>%
  group_by(tahun,bulan)%>%
  summarise(BuyUSD = mean(`Kurs Jual`))%>%
  mutate_at(c("tahun","bulan"), as.numeric)%>%
  arrange(tahun,bulan)
dc$USD = USD_buy$BuyUSD

Manipulating Data

Data unemplyment rate dan GDP didapat bukan per-bulan. Data unemployment rate dapat data per-6 bulan (semester), dan data GDP dapat data per-3 bulan (triwulan). Maka dari itu, diperlukan manipulasi data untuk menentukan nilai per-bulannya. data tersebut akan dimanipulasi menggunakan rnorm() dengan mean merupakan data yang didapat, dan standar deviasi yang dihitung per-tahunnya berapa.

unemployment_rate <- read_excel("data_UAS/Pengangguran.xlsx", sheet="Sheet2")
data.table(unemployment_rate)
unemployment_rate <- read_excel("data_UAS/Pengangguran.xlsx", sheet="Sheet2")
ur <- unemployment_rate$unemployment_rate

for (j in c(1,3,5,7,9,11,13,15,17,19,21)) {
  print(sd(ur[j : (j+1)]))
}
## [1] 0.1909188
## [1] 0.3676955
## [1] 0.1697056
## [1] 0.205061
## [1] 0.1697056
## [1] 0.2616295
## [1] 0.07778175
## [1] 0.1202082
## [1] 0.1414214
## [1] 0.1767767
## [1] 1.506137
ali <- data.frame(
   0.1909188,
   0.3676955,
   0.1697056,
   0.205061,
   0.1697056,
   0.2616295,
   0.07778175,
   0.1202082,
   0.1414214,
   0.1767767,
   1.506137
)

set.seed(1993)

a= c(
rnorm(6, mean= ur[1] , sd = ali[1,1]),
rnorm(6, mean= ur[2] , sd = ali[1,1]),
rnorm(6, mean= ur[3] , sd = ali[1,2]),
rnorm(6, mean= ur[4] , sd = ali[1,2]),
rnorm(6, mean= ur[5] , sd = ali[1,3]),
rnorm(6, mean= ur[6] , sd = ali[1,3]),
rnorm(6, mean= ur[7] , sd = ali[1,4]),
rnorm(6, mean= ur[8] , sd = ali[1,4]),
rnorm(6, mean= ur[9] , sd = ali[1,5]),
rnorm(6, mean= ur[10], sd = ali[1,5]),
rnorm(6, mean= ur[11], sd = ali[1,6]),
rnorm(6, mean= ur[12], sd = ali[1,6]),
rnorm(6, mean= ur[13], sd = ali[1,7]),
rnorm(6, mean= ur[14], sd = ali[1,7]),
rnorm(6, mean= ur[15], sd = ali[1,8]),
rnorm(6, mean= ur[16], sd = ali[1,8]),
rnorm(6, mean= ur[17], sd = ali[1,9]),
rnorm(6, mean= ur[18], sd = ali[1,9]),
rnorm(6, mean= ur[19], sd = ali[1,10]),
rnorm(6, mean= ur[20], sd = ali[1,10]),
rnorm(6, mean= ur[21], sd = ali[1,11]),
rnorm(6, mean= ur[22], sd = ali[1,11])
)

c <- c(rnorm(12, mean= mean(ur), sd = sd(ur)) )

dc$unemployment_rate = c(a,c)
percentage <- c(5.69,6.17,5.8,6.9,
                  6.48,6.27,6.01,5.94,
                  6.11,6.21,5.94,5.87,
                  5.54,5.59,5.52,5.58,
                  5.12,4.94,4.93,5.05,
                  4.83,4.74,4.78,5.15,
                  4.94,5.21,5.03,4.94,
                  5.01,5.01,5.06,5.19,
                  5.06,5.27,5.17,5.18,
                  5.07,5.05,5.02,4.97,
                  2.97,-5.32,-3.49,-2.19,
                  -0.69,7.08,3.53,5.03)
sd_gdp <- c(
 sd(c(5.69,6.17,5.8,6.9)),
 sd(c(6.48,6.27,6.01,5.94)),
 sd(c(6.11,6.21,5.94,5.87)),
 sd(c(5.54,5.59,5.52,5.58)),
 sd(c(5.12,4.94,4.93,5.05)),
 sd(c(4.83,4.74,4.78,5.15)),
 sd(c(4.94,5.21,5.03,4.94)),
 sd(c(5.01,5.01,5.06,5.19)),
 sd(c(5.06,5.27,5.17,5.18)),
 sd(c(5.07,5.05,5.02,4.97)),
 sd(c(2.97,-5.32,-3.49,-2.19)),
 sd(c(-0.69,7.08,3.53,5.03)))
 
GDP_growth <- c(
rnorm(3, mean = percentage[1],  sd= sd_gdp[1]),
rnorm(3, mean = percentage[2],  sd= sd_gdp[1]),
rnorm(3, mean = percentage[3],  sd= sd_gdp[1]),
rnorm(3, mean = percentage[4],  sd= sd_gdp[1]),
rnorm(3, mean = percentage[5],  sd= sd_gdp[2]),
rnorm(3, mean = percentage[6],  sd= sd_gdp[2]),
rnorm(3, mean = percentage[7],  sd= sd_gdp[2]),
rnorm(3, mean = percentage[8],  sd= sd_gdp[2]),
rnorm(3, mean = percentage[9],  sd= sd_gdp[3]),
rnorm(3, mean = percentage[10], sd= sd_gdp[3]),
rnorm(3, mean = percentage[11], sd= sd_gdp[3]),
rnorm(3, mean = percentage[12], sd= sd_gdp[3]),
rnorm(3, mean = percentage[13], sd= sd_gdp[4]),
rnorm(3, mean = percentage[14], sd= sd_gdp[4]),
rnorm(3, mean = percentage[15], sd= sd_gdp[4]),
rnorm(3, mean = percentage[16], sd= sd_gdp[4]),
rnorm(3, mean = percentage[17], sd= sd_gdp[5]),
rnorm(3, mean = percentage[18], sd= sd_gdp[5]),
rnorm(3, mean = percentage[19], sd= sd_gdp[5]),
rnorm(3, mean = percentage[20], sd= sd_gdp[5]),
rnorm(3, mean = percentage[21], sd= sd_gdp[6]),
rnorm(3, mean = percentage[22], sd= sd_gdp[6]),
rnorm(3, mean = percentage[23], sd= sd_gdp[6]),
rnorm(3, mean = percentage[24], sd= sd_gdp[6]),
rnorm(3, mean = percentage[25], sd= sd_gdp[7]),
rnorm(3, mean = percentage[26], sd= sd_gdp[7]),
rnorm(3, mean = percentage[27], sd= sd_gdp[7]),
rnorm(3, mean = percentage[28], sd= sd_gdp[7]),
rnorm(3, mean = percentage[29], sd= sd_gdp[8]),
rnorm(3, mean = percentage[30], sd= sd_gdp[8]),
rnorm(3, mean = percentage[31], sd= sd_gdp[8]),
rnorm(3, mean = percentage[32], sd= sd_gdp[8]),
rnorm(3, mean = percentage[33], sd= sd_gdp[9]),
rnorm(3, mean = percentage[34], sd= sd_gdp[9]),
rnorm(3, mean = percentage[35], sd= sd_gdp[9]),
rnorm(3, mean = percentage[36], sd= sd_gdp[9]),
rnorm(3, mean = percentage[37], sd= sd_gdp[10]),
rnorm(3, mean = percentage[38], sd= sd_gdp[10]),
rnorm(3, mean = percentage[39], sd= sd_gdp[10]),
rnorm(3, mean = percentage[40], sd= sd_gdp[10]),
rnorm(3, mean = percentage[41], sd= sd_gdp[11]),
rnorm(3, mean = percentage[42], sd= sd_gdp[11]),
rnorm(3, mean = percentage[43], sd= sd_gdp[11]),
rnorm(3, mean = percentage[44], sd= sd_gdp[11]),
rnorm(3, mean = percentage[45], sd= sd_gdp[12]),
rnorm(3, mean = percentage[46], sd= sd_gdp[12]),
rnorm(3, mean = percentage[47], sd= sd_gdp[12]),
rnorm(3, mean = percentage[48], sd= sd_gdp[12])
)
dc$GDP_growth = GDP_growth
GDP_growth<- data.frame(GDP_growth)

Data Exploration

  • Load the dataset into your preferred statistical software (e.g., R, Python, or Excel).

  • Examine the structure of the dataset, including variable names and data types.

  • Explore summary statistics and visualize the distributions of each variable.

data.table(dc)
summary(dc)
##      dates            interest_rate    stock_index     consumer_confidence
##  Min.   :2010-01-01   Min.   :3.500   Min.   : 496.0   Min.   : 77.3      
##  1st Qu.:2012-12-24   1st Qu.:4.750   1st Qu.: 712.9   1st Qu.:108.6      
##  Median :2015-12-16   Median :5.750   Median : 839.3   Median :114.3      
##  Mean   :2015-12-16   Mean   :5.788   Mean   : 826.6   Mean   :113.1      
##  3rd Qu.:2018-12-08   3rd Qu.:6.750   3rd Qu.: 934.1   3rd Qu.:119.3      
##  Max.   :2021-12-01   Max.   :7.750   Max.   :1105.8   Max.   :128.0      
##  inflation_rate        USD        unemployment_rate   GDP_growth    
##  Min.   :-0.610   Min.   : 8575   Min.   :2.731     Min.   :-8.840  
##  1st Qu.: 0.890   1st Qu.: 9725   1st Qu.:5.427     1st Qu.: 4.917  
##  Median : 1.710   Median :13367   Median :5.843     Median : 5.100  
##  Mean   : 2.112   Mean   :12305   Mean   :5.949     Mean   : 4.483  
##  3rd Qu.: 2.663   3rd Qu.:14134   3rd Qu.:6.511     3rd Qu.: 5.642  
##  Max.   : 8.380   Max.   :15947   Max.   :8.525     Max.   :13.479
hchart(density(dc$interest_rate), type = "area", color= "#003285",name = "interest rate")%>%
  hc_add_series(density(dc$inflation_rate), type="area", color = "#FF7F3E", name="inflation rate")%>%
  hc_add_series(density(dc$unemployment_rate), type="area", color = "#94FFD8", name="unemployment rate")%>%
  hc_add_series(density(dc$GDP_growth), type="area", color = "#FF0080", name="GDP growth")%>%
  hc_title(text = "Density Plots of Rate",
           style = list(color = "#240750", useHTML = TRUE))
hchart(density(dc$stock_index), type="area", color="#A0DEFF", name="LQ45 Stock Index")%>%
  hc_title(text = "Density Plots of Stock Index LQ45",
           style = list(color = "#240750", useHTML = TRUE))
hchart(density(dc$USD), type="area", color="#A1DD70", name="USD to Rupiah") %>%
  hc_title(text = "Density Plots of USD",
           style = list(color = "#240750", useHTML = TRUE))

Dari semua plot density, bisa disimpulkan bahwa tidak ada variabel yang memiliki distribusi normal.

Correlation Analysis

  • Calculate the correlation matrix to examine the relationships between pairs of economic indicators.

  • Visualize the correlations using a heatmap or correlation plot.

cor(dc[,-1])*100 # correlation in percentage
##                     interest_rate stock_index consumer_confidence
## interest_rate          100.000000  -39.776376           -3.407922
## stock_index            -39.776376  100.000000           56.863719
## consumer_confidence     -3.407922   56.863719          100.000000
## inflation_rate          30.672687  -21.685176            1.808067
## USD                    -45.064439   73.409723           14.416465
## unemployment_rate       22.923106  -64.800773          -30.377054
## GDP_growth              39.114612   -6.063271           44.334155
##                     inflation_rate       USD unemployment_rate GDP_growth
## interest_rate            30.672687 -45.06444          22.92311  39.114612
## stock_index             -21.685176  73.40972         -64.80077  -6.063271
## consumer_confidence       1.808067  14.41646         -30.37705  44.334155
## inflation_rate          100.000000 -29.91644          19.96567  22.724255
## USD                     -29.916442 100.00000         -63.09429 -42.762542
## unemployment_rate        19.965674 -63.09429         100.00000  18.538719
## GDP_growth               22.724255 -42.76254          18.53872 100.000000

Jika dilihat dari hasil korelasi, kita bisa memilih variabel untuk dijadikan variabel dependen. Variabel dependen yang akan dipilih adalah variabel yang memiliki korelasi sedang pada variabel tersebut. Variabel yang bisa dipilih adalah rata-rata kurs dollar amerika. Kurs USD memiliki nilai korelasi dengan variabel dependen yang lumayan besar.

Time Series Analysis

  • Conduct time series analysis for each economic indicator, including:

  • Checking for stationarity using appropriate tests (e.g., Augmented Dickey-Fuller test).

  • Identifying trends, seasonality, and autocorrelation patterns.

  • Applying transformations if necessary to achieve stationarity (e.g., first differencing).

  • Building time series models (e.g., ARIMA, GARCH) to forecast future values of each indicator.

dc %>%
  hchart(., type = "line", 
         hcaes(x= dates, y= USD), name="Rupiah to US", color="#FEAE6F")%>%
  hc_yAxis(opposite = TRUE,
         labels = list(format = "{value}"))%>%
  hc_title(    text = "Rupiah to US Dollar Over The Years",
    margin = 20,
    align = "left",
    style = list(color = "#615EFC", useHTML = TRUE))

Stationarity

Untuk mengecek stasionaritas data dari time series, kita akan menggunakan Dicky-Fuller Test.
Hypothesis :

  • H0 = data tidak stasioner

  • H1 = data stasioner


kita akan menguji menggunakan coding adf.test() menggunakan library tseries.

library(forecast)
library(tseries)

adf.test(dc$USD) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dc$USD
## Dickey-Fuller = -1.4491, Lag order = 5, p-value = 0.8058
## alternative hypothesis: stationary

p-value dari adf test adalah 0.8058 yang berarti data tidak stasioner. Maka dari itu, diperlukan differencing agar stasioner.

adf.test(diff(dc$USD))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(dc$USD)
## Dickey-Fuller = -4.8683, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

P-value adf test sudah < , maka data sudah stasioner. Karena data membutuhkan differencing satu kali, I dari model akan 1.
Lalu, Lanjut kita melihat plot autokorelasinya. Plot yang biasa digunakan adalah ACF dan PACF. ACF adalah autocorrelation function plot, dan PACF adalah partial autocorrelarion function. PLot” ini membandingkan nilai autokorelasi dengan lag-lag sebelumnya. Suatu data sudha stasioner terhadap varians ketika model acf atau pacf plot akan membentuk white noise. white noise adalah ketika 95% dari nilai korelasi antar lag berada di antara critical value.

acf(diff(dc$USD), lag= 144, pl=T)

Dari hasil Plot ACF, kita bisa menentukan berapa banyak AR(p). Ketika dilihat, 4 ada korelasi yang lebih dari batasan critical value. Critical value adalah garis yang berwarna bitu putus”. Modelnya sudah white noise. Dari hasil ACF ini, disimpulkan datanya memiliki AR(p) dengan p=1-4.

Lalu, lanjut ke plot PACF. PACF bisa menentukan berapa banyak MA(q).

pacf(diff(dc$USD), lag =144,pl=T)

Dari hasil grafik PACF, ada 1 lonjakan korelasi yang lebih dari critical value. Maka, MA(q) yang cocok adalah q=1.

Hasil dari subbab stasionarity ini adalah adalah :

  1. p = 1-4

  2. d = 1

  3. q = 1


### decompose

ts_dc <- ts((dc$USD), start = c(2010,1), end = c(2021,12), frequency = 12)
decomp <- decompose(ts_dc)
plot(decomp)

Pada data, memiliki trend naik, dan tidak seasonal.

Box Jenkins Method

Box-Jenkins Method biasa disebut dengan model ARIMA. Sekarang kita membuat model dengan hasil stasionerity. Lalu kita juga akan membuat ARIMA dengan model auto dengan auto.arima().

dc_arima1 <- arima(dc$USD, order = c(1,1,1)) # ARIMA(1,1,1)
dc_arima2 <- arima(dc$USD, order = c(2,1,1)) # ARIMA(2,1,1)
dc_arima3 <- arima(dc$USD, order = c(3,1,1)) # ARIMA(3,1,1)
dc_arima4 <- arima(dc$USD, order = c(4,1,1)) # ARIMA(4,1,1)

auto_arima <- auto.arima(dc$USD, max.p = 4, d = 1, max.q = 1) 
auto_arima
## Series: dc$USD 
## ARIMA(2,1,0) with drift 
## 
## Coefficients:
##          ar1      ar2    drift
##       0.2764  -0.3590  35.9005
## s.e.  0.0778   0.0775  18.4398
## 
## sigma^2 = 57880:  log likelihood = -985.63
## AIC=1979.26   AICc=1979.55   BIC=1991.11

Dari hasil model auto arima, didapatkan ARIMA(2,1,0).

ts_error <-data.frame(
  actual = dc$USD
)
ts_error$"ARIMA(1,1,1)" =dc_arima1$residuals
ts_error$"ARIMA(2,1,1)" =dc_arima2$residuals
ts_error$"ARIMA(3,1,1)" =dc_arima3$residuals
ts_error$"ARIMA(4,1,1)" =dc_arima4$residuals
ts_error$"ARIMA(2,1,0)" =auto_arima$residuals

ARCH and GARCH

ARCH (autoregressive Conditional Heteroscedasity) dan GARCH (General Autoregressive Conditional Heteroscedasity). Model ARCH dan GARCH digunakan untuk model yang memiliki

library(FinTS)
ArchTest(dc$USD, d=1)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  dc$USD
## Chi-squared = 122.44, df = 12, p-value < 2.2e-16
library(rugarch)
# ARCH(1) model
arch1_spec <- ugarchspec(variance.model= list(model="sGARCH", garchOrder = c(1,0)), mean.model = list(armaOrder = c(0,0), include.mean=T))
arch1_fit <- ugarchfit(spec=arch1_spec, data=dc$USD)

ts_error$"ARCH(1)" = residuals(arch1_fit)
# GARCH(1) model
garch1_spec <- ugarchspec(variance.model= list(model="sGARCH", garchOrder = c(1,1)), mean.model = list(armaOrder = c(0,0), include.mean=T))
garch1_fit <- ugarchfit(spec=garch1_spec, data=dc$USD)

ts_error$"GARCH(1,1)" = residuals(garch1_fit)
PE_ts <-abs(ts_error)/ ts_error$actual
PE_ts_tes <-PE_ts[c(116:144),] 

acc_ts <- data.frame(
  Method = colnames(PE_ts_tes[,-1]),
  MAPE = c(mean(PE_ts_tes$`ARIMA(1,1,1)`) *100,
           mean(PE_ts_tes$`ARIMA(2,1,1)`) *100,
           mean(PE_ts_tes$`ARIMA(3,1,1)`) *100,
           mean(PE_ts_tes$`ARIMA(4,1,1)`) *100,
           mean(PE_ts_tes$`ARIMA(2,1,0)`) *100,
           mean(PE_ts_tes$`ARCH(1)`)*100,
           mean(PE_ts_tes$`GARCH(1,1)`*100)
  )
)

ts_rsq <- data.frame(
  actual = dc$USD
)
ts_rsq$"ARIMA(1,1,1)" = fitted.values(dc_arima1)
ts_rsq$"ARIMA(2,1,1)" = fitted.values(dc_arima2)
ts_rsq$"ARIMA(3,1,1)" = fitted.values(dc_arima3)
ts_rsq$"ARIMA(4,1,1)" = fitted.values(dc_arima4)
ts_rsq$"ARIMA(2,1,0)" = fitted.values(auto_arima)
ts_rsq$"ARCH(1)" = sigma(arch1_fit)
ts_rsq$"GARCH(1,1)"= sigma(garch1_fit)

ts_rsquare<-as.data.frame(cor(ts_rsq[c(116:144),-1], ts_rsq$actual[116:144])^2)
acc_ts$rsq<- ts_rsquare$V1

Holts Method

Holt’s Method adalah metode lanjutan dari eksponential smoothing. hanya Holt’s bisa memasukkan trend ke dalam model. Modelnya akan dibuat sebagai berikut :

library(fpp2)
holt_train <- window(dc$USD, end = 115)
holt_test <- window(dc$USD, start = 116)

beta <- seq(.0001, .5, by = .001)
RMSE <- NA
for(i in seq_along(beta)) {
  fit <- holt(holt_train,
              beta = beta[i], 
              h = 100)
  RMSE[i] <- accuracy(fit, 
                      holt_test)[2,2]
}

beta.fit <- data_frame(beta, RMSE)
beta.min <- filter(beta.fit, 
                   RMSE == min(RMSE))

holts <- holt(holt_train, h=100)
holts2 <- holt(holt_train, h=100, beta=beta.min$beta)

accuracy(holts , holt_test)
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set   6.225125 207.9699 144.3144 0.08898472 1.186428 0.9676789
## Test set     183.318828 473.7263 291.6796 1.18509676 1.957945 1.9558150
##                   ACF1 Theil's U
## Training set 0.1881446        NA
## Test set     0.5611894  1.137742
accuracy(holts2, holt_test)
##                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set 11.38941 208.2707 145.4785 0.1395392 1.197545 0.9754848 0.1986061
## Test set     82.38314 458.6395 306.1691 0.4868101 2.068370 2.0529723 0.5833957
##              Theil's U
## Training set        NA
## Test set      1.104712
(cor(holts$fitted[1:29], holt_test))^2
## [1] 0.006535874
(cor(holts2$fitted[1:29], holt_test))^2
## [1] 0.006257628
holts_acc <- data.frame(
  Method = c("Holt's b=0.0568", "Holt's b=0.0451"),
  MAPE = c(1.957945, 2.068370),
  rsq = c((cor(holts$fitted[1:29], holt_test))^2,(cor(holts2$fitted[1:29], holt_test))^2)
  )

acc_ts <- rbind(acc_ts, holts_acc)

Regression Analysis

  • Formulate a regression model with one or more economic indicators as dependent variables and other indicators as independent variables.

  • Check for multicollinearity among independent variables.

  • Estimate the regression model using appropriate techniques (e.g., Ordinary Least Squares).

  • Assess the significance and interpret the coefficients of the independent variables.

multicollinearity

1/(1-cor(dc[,-1])^2) # multicollinearity vif
##                     interest_rate stock_index consumer_confidence
## interest_rate                 Inf    1.187953            1.001163
## stock_index              1.187953         Inf            1.477865
## consumer_confidence      1.001163    1.477865                 Inf
## inflation_rate           1.103852    1.049345            1.000327
## USD                      1.254832    2.168721            1.021225
## unemployment_rate        1.055461    1.723882            1.101657
## GDP_growth               1.180631    1.003690            1.244635
##                     inflation_rate      USD unemployment_rate GDP_growth
## interest_rate             1.103852 1.254832          1.055461   1.180631
## stock_index               1.049345 2.168721          1.723882   1.003690
## consumer_confidence       1.000327 1.021225          1.101657   1.244635
## inflation_rate                 Inf 1.098297          1.041518   1.054451
## USD                       1.098297      Inf          1.661375   1.223786
## unemployment_rate         1.041518 1.661375               Inf   1.035592
## GDP_growth                1.054451 1.223786          1.035592        Inf

Semua nilai VIF di bawah 10, hal ini menunjukkan bahwa data tidak memiliki multikolinearitas antar variabel.

linear regression

with dates

# membuat regresi dengan semua variabel
lm_2.1 = lm(USD~., data=dc)
summary(lm_2.1) 
## 
## Call:
## lm(formula = USD ~ ., data = dc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1156.29  -503.91   -44.94   566.57  1372.06 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -2.104e+04  2.009e+03 -10.472  < 2e-16 ***
## dates                1.900e+00  9.556e-02  19.878  < 2e-16 ***
## interest_rate        5.478e+02  5.622e+01   9.743  < 2e-16 ***
## stock_index         -1.999e-01  8.514e-01  -0.235  0.81473    
## consumer_confidence -1.709e+00  8.704e+00  -0.196  0.84464    
## inflation_rate      -1.225e+01  2.906e+01  -0.421  0.67411    
## unemployment_rate   -2.197e+02  7.877e+01  -2.789  0.00604 ** 
## GDP_growth          -4.134e+00  2.244e+01  -0.184  0.85414    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 595.2 on 136 degrees of freedom
## Multiple R-squared:  0.9296, Adjusted R-squared:  0.926 
## F-statistic: 256.5 on 7 and 136 DF,  p-value: < 2.2e-16

dari hasil plot, variabel yang siknifikan adalah dates, interest rate, dan unemployment_rate. Lalu, kita akan coba membuat regrei haya dengan variabel tersebut

# membuat regresi linear baru hanya dengan variabele yang siknifikan 
lm_2.2 = lm(USD~dates+interest_rate+unemployment_rate, data=dc)
summary(lm_2.2)
## 
## Call:
## lm(formula = USD ~ dates + interest_rate + unemployment_rate, 
##     data = dc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1138.63  -520.03   -56.62   553.25  1431.57 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -2.143e+04  1.593e+03 -13.447   <2e-16 ***
## dates              1.895e+00  6.482e-02  29.238   <2e-16 ***
## interest_rate      5.416e+02  5.425e+01   9.982   <2e-16 ***
## unemployment_rate -2.039e+02  7.151e+01  -2.852    0.005 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 588.3 on 140 degrees of freedom
## Multiple R-squared:  0.9292, Adjusted R-squared:  0.9277 
## F-statistic: 612.5 on 3 and 140 DF,  p-value: < 2.2e-16

Dari hasil summary regresi di atas, semua p-value model kurang dari . Pada lm_2.1 masih terdapat variabel yang tidak siknifikan, lm_2.2 semua variabel siknifikan. Meski dari hal tersebut, nilai r-squarenya masih lebih besar lm_2.1. Meskipun lm_2.2semua variabelnya siknifikan, perbedaannya hanya sebesar 0.04%. R-square dari regresi tersebut adalah lm_2.1 dan lm_2.2. Kedua regresi memiliki R-square 92.96% dan 92.92%. Dari hasil kedua regresi ini kita bisa curiga, mengapa regresi bisa seakurat itu dengan variabel yang sangat sedikit? Kita harus coba cek lagi dengan tidak menggunakan variabel dates/waktu.

cor(dc[,-1])
##                     interest_rate stock_index consumer_confidence
## interest_rate          1.00000000 -0.39776376         -0.03407922
## stock_index           -0.39776376  1.00000000          0.56863719
## consumer_confidence   -0.03407922  0.56863719          1.00000000
## inflation_rate         0.30672687 -0.21685176          0.01808067
## USD                   -0.45064439  0.73409723          0.14416465
## unemployment_rate      0.22923106 -0.64800773         -0.30377054
## GDP_growth             0.39114612 -0.06063271          0.44334155
##                     inflation_rate        USD unemployment_rate  GDP_growth
## interest_rate           0.30672687 -0.4506444         0.2292311  0.39114612
## stock_index            -0.21685176  0.7340972        -0.6480077 -0.06063271
## consumer_confidence     0.01808067  0.1441646        -0.3037705  0.44334155
## inflation_rate          1.00000000 -0.2991644         0.1996567  0.22724255
## USD                    -0.29916442  1.0000000        -0.6309429 -0.42762542
## unemployment_rate       0.19965674 -0.6309429         1.0000000  0.18538719
## GDP_growth              0.22724255 -0.4276254         0.1853872  1.00000000
lm_2.nd = lm(USD~interest_rate+unemployment_rate, data=dc)
summary(lm_2.nd)
## 
## Call:
## lm(formula = USD ~ interest_rate + unemployment_rate, data = dc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2682.2  -907.2    14.1   875.4  5596.5 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        23870.0      989.7  24.117  < 2e-16 ***
## interest_rate       -550.2      104.5  -5.263 5.16e-07 ***
## unemployment_rate  -1408.7      155.2  -9.075 8.98e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1563 on 141 degrees of freedom
## Multiple R-squared:  0.4969, Adjusted R-squared:  0.4898 
## F-statistic: 69.64 on 2 and 141 DF,  p-value: < 2.2e-16

Setelah dicek, kita bisa lihat bahwa dates memiliki pengaruh yang sangat siknifikan terhadap USD, sedangkan jika untuk memprediksi waktu pastinya akan lebih cocok model time series. Untuk regresi, kita akan menggunakan vaiabel lain untuk mencari variabel lain yang memiliki pengaruh siknifikan terhadap variabel dependen.

summary(lm(interest_rate~., data=dc[,-1]))$r.squared
## [1] 0.3309123
summary(lm(stock_index~., data=dc[,-1]))$r.squared
## [1] 0.7834162
summary(lm(consumer_confidence~., data=dc[,-1]))$r.squared
## [1] 0.5807829
summary(lm(inflation_rate~., data=dc[,-1]))$r.squared
## [1] 0.1332784
summary(lm(unemployment_rate~., data=dc[,-1]))$r.squared
## [1] 0.4869431
summary(lm(GDP_growth~., data=dc[,-1]))$r.squared
## [1] 0.4953778
summary(lm(interest_rate~., data=dc[,-1]))$r.squared
## [1] 0.3309123
summary(lm(USD~., data=dc[,-1]))$r.squared
## [1] 0.7250227

Kita akan coba membuat 2 model, yaitu untuk memprediksi stock index dan USD.

Predict USD

Lalu apakah model terikat dengan tanggal, maka dari itu kita akan coba membuat regresi tanpa dates.

# regresi menggunakan semua variabel
lm_2.4 =lm(USD~., data=dc[,-1])
summary(lm_2.4)
## 
## Call:
## lm(formula = USD ~ ., data = dc[, -1])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2653.7  -652.9  -205.2   663.7  3028.7 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         12417.879   2160.442   5.748 5.62e-08 ***
## interest_rate         -55.482     93.187  -0.595  0.55257    
## stock_index            11.236      1.236   9.091 9.92e-16 ***
## consumer_confidence   -48.490     16.499  -2.939  0.00387 ** 
## inflation_rate        -52.939     57.072  -0.928  0.35526    
## unemployment_rate    -459.887    153.261  -3.001  0.00320 ** 
## GDP_growth           -166.909     41.147  -4.056 8.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1172 on 137 degrees of freedom
## Multiple R-squared:  0.725,  Adjusted R-squared:  0.713 
## F-statistic:  60.2 on 6 and 137 DF,  p-value: < 2.2e-16

p-value <, variabel yang siknifikan terhadap USD adalah stock_index, consumer_confidence, inflation rate, unemployment rate dan GDP growth. R-square model ini berbeda jauh dengan model dengan date. R-square hanya 72.5%, yang berarti model hanya dapat menjelaskan 72.5%. Masih ada 27.5% variabel yang velum masuk model. Sepertinya variansi tersebut disebabkan oleh variabel dates.
Mari kita coba bentuk model dengan variabel tersebut.

# regresi dengan variabel yang siknifikan
lm_2.5 =lm(USD~stock_index+consumer_confidence+unemployment_rate+GDP_growth, data=dc[,-1])
summary(lm_2.5)
## 
## Call:
## lm(formula = USD ~ stock_index + consumer_confidence + unemployment_rate + 
##     GDP_growth, data = dc[, -1])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2659.4  -710.6  -233.0   679.9  3285.4 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         11797.813   1963.266   6.009 1.55e-08 ***
## stock_index            11.657      1.153  10.108  < 2e-16 ***
## consumer_confidence   -49.717     16.434  -3.025  0.00296 ** 
## unemployment_rate    -453.455    151.338  -2.996  0.00324 ** 
## GDP_growth           -180.401     38.554  -4.679 6.76e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1170 on 139 degrees of freedom
## Multiple R-squared:  0.7221, Adjusted R-squared:  0.7141 
## F-statistic: 90.29 on 4 and 139 DF,  p-value: < 2.2e-16

Dari hasil variabel tersebut, didapatkan r-squarenya adalah 0.7221 dan semua variabel independen siknifikan terhadap US dollar. Lalu, apa yang harus dilakukan dengan 2 variabel yang tidak masuk kedalam model? Kita akan coba menambahkan dengan menggunakan interaksi dari variabel agar hampir semua variabel bisa masuk ke dalam model. Setelah itu, kita akan coba lihat lagi apakah terdapat variabel lain yang memiliki korelasi yang lumayan tinggi dan diinteraksikan dengan variabel yang sudah siknifikan terhadap USD. 

cor(dc[,-c(1,8)])
##                     interest_rate stock_index consumer_confidence
## interest_rate          1.00000000  -0.3977638         -0.03407922
## stock_index           -0.39776376   1.0000000          0.56863719
## consumer_confidence   -0.03407922   0.5686372          1.00000000
## inflation_rate         0.30672687  -0.2168518          0.01808067
## USD                   -0.45064439   0.7340972          0.14416465
## unemployment_rate      0.22923106  -0.6480077         -0.30377054
##                     inflation_rate        USD unemployment_rate
## interest_rate           0.30672687 -0.4506444         0.2292311
## stock_index            -0.21685176  0.7340972        -0.6480077
## consumer_confidence     0.01808067  0.1441646        -0.3037705
## inflation_rate          1.00000000 -0.2991644         0.1996567
## USD                    -0.29916442  1.0000000        -0.6309429
## unemployment_rate       0.19965674 -0.6309429         1.0000000

dilihat dari hasil korelasi, interest rate mempunyai korelasi yang besar pada stock index dan gdp growth sebesar 39%. Inflation rate memiliki korelasi dengan interest rate sebesar 30%. Karena korelasi inflation rate yang paling tinggi dengan interest rate dan hasil nilai korelasi cukup kecil, maka kita akan coba menggunakan interest saja yang mana interest rate sudah bisa mewakili sedikit inflation rate.

# mencoba dengan menambahkan interkasi interest rate dengan gdp growth
lm_2.6 =lm(USD~stock_index+consumer_confidence+unemployment_rate+interest_rate:GDP_growth, data=dc[,-1])
summary(lm_2.6)
## 
## Call:
## lm(formula = USD ~ stock_index + consumer_confidence + unemployment_rate + 
##     interest_rate:GDP_growth, data = dc[, -1])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2654.3  -660.7  -120.5   611.7  2957.7 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              12709.134   1847.737   6.878 1.89e-10 ***
## stock_index                 10.594      1.195   8.868 3.24e-15 ***
## consumer_confidence        -47.018     16.004  -2.938  0.00387 ** 
## unemployment_rate         -472.178    148.101  -3.188  0.00177 ** 
## interest_rate:GDP_growth   -37.603      7.217  -5.210 6.66e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1151 on 139 degrees of freedom
## Multiple R-squared:  0.7309, Adjusted R-squared:  0.7231 
## F-statistic: 94.37 on 4 and 139 DF,  p-value: < 2.2e-16

model tersebut menambah nilai r-square model dari 72.21% menjadi 73.09% ketika interaksi interest rate dan GDP growth ditambahkan. Interaksi tersebut memiliki dampak yang sangat siknifikan terhadap USD. Apakah yang akan terjadi ketika kita menginteraksikan interest ratae dengan stock index?

# mencoba dengan menambahkan interkasi interest rate dengan stock index
lm_2.7 =lm(USD~interest_rate:stock_index+consumer_confidence+unemployment_rate+GDP_growth, data=dc[,-1])
summary(lm_2.7)
## 
## Call:
## lm(formula = USD ~ interest_rate:stock_index + consumer_confidence + 
##     unemployment_rate + GDP_growth, data = dc[, -1])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2522.1 -1117.3   -59.7   796.9  5173.1 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.592e+04  2.579e+03   6.171 7.00e-09 ***
## consumer_confidence        4.134e+01  1.811e+01   2.282    0.024 *  
## unemployment_rate         -1.240e+03  1.745e+02  -7.103 5.77e-11 ***
## GDP_growth                -2.932e+02  5.075e+01  -5.777 4.78e-08 ***
## interest_rate:stock_index  8.553e-02  1.347e-01   0.635    0.527    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1539 on 139 degrees of freedom
## Multiple R-squared:  0.5192, Adjusted R-squared:  0.5053 
## F-statistic: 37.52 on 4 and 139 DF,  p-value: < 2.2e-16

Model interaksi interest rate dan stock index mempunyai nilai r-square yang anjlok menjadi 0.5192, yang berarti terdapat 48.08% bagian model yang hilang yang belum masuk ke dalam model.

Predict Stock Index

lm_2.8 <- lm(stock_index~., data=dc[,-1])
summary(lm_2.8)
## 
## Call:
## lm(formula = stock_index ~ ., data = dc[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -245.436  -38.593   -1.368   45.534  122.669 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -53.698355 131.326967  -0.409  0.68326    
## interest_rate       -13.506177   4.961862  -2.722  0.00733 ** 
## consumer_confidence   6.070463   0.770466   7.879 9.15e-13 ***
## inflation_rate       -0.439251   3.125487  -0.141  0.88844    
## USD                   0.033490   0.003684   9.091 9.92e-16 ***
## unemployment_rate   -25.845597   8.351073  -3.095  0.00239 ** 
## GDP_growth            3.261679   2.361192   1.381  0.16941    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 63.99 on 137 degrees of freedom
## Multiple R-squared:  0.7834, Adjusted R-squared:  0.7739 
## F-statistic: 82.59 on 6 and 137 DF,  p-value: < 2.2e-16

dari hasil model stock index, model memliki r-square sebesar 0.7834. dan variabel yang memiliki pengaruh siknifikan adalah interest rate, consumer confidence, unemployment rate dan USD. Model sudah memiliki r-square yang lebih tinggi dari memprediksi USD, tetapi intercept dari model ini tidak memiliki pengaruh siknifikan. Mari kita coba dengan variabel yang berpengaruh.

lm_2.9 <- lm(stock_index~interest_rate + consumer_confidence+unemployment_rate+USD, data=dc[,-1])
summary(lm_2.9)
## 
## Call:
## lm(formula = stock_index ~ interest_rate + consumer_confidence + 
##     unemployment_rate + USD, data = dc[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -239.391  -37.032   -1.212   40.220  125.021 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1.024e+02  1.264e+02  -0.810  0.41936    
## interest_rate       -1.177e+01  4.682e+00  -2.513  0.01311 *  
## consumer_confidence  6.671e+00  6.323e-01  10.550  < 2e-16 ***
## unemployment_rate   -2.512e+01  8.320e+00  -3.019  0.00302 ** 
## USD                  3.187e-02  3.452e-03   9.232    4e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 63.97 on 139 degrees of freedom
## Multiple R-squared:  0.7804, Adjusted R-squared:  0.7741 
## F-statistic: 123.5 on 4 and 139 DF,  p-value: < 2.2e-16

R-square model sebanyak 0.7804, jika dibandingkan dengan model sebelumnya, r-square berkurang sebanyak 0.003. Intercept masih tidak memiliki pengaruh siknifikan. Saya memutuskan untuk tidak menggunakan model ini karena intercept dari model ini tidak memiliki pengaruh yang siknifikan. Meskipun model ini memiliki R-square yang lebih bbesar dari model prediksi USD.

Summary

Model perbedaan r-square model bisa mencapai 49%. Maka bisa dikatakan bahwa dates memiliki pengaruh yang lumayan banyak dan dates merupakan faktor yang siknifikan.
Dari hasil prediksi USD dan stock index, model yang akan digunakan adalh model yang terbaik dari prediksi USD. Meskipun mdoel prediksi USD memiliki r-square yang lebih rendah, semua variabel termasuk intercept memiliki pengaruh siknifikan terhadap . Intercept model prediksi stock index tidak berpengaruh siknifikan. Saya memutuskan untuk tidak menggunakan model prediksi stock index.
Dari semua model regresi linear tanpa variabel dates, nilai r-square modelnya adalah 72.5% untuk lm_2.4, 72.21% untuk lm_2.5, dan 73.09% untuk lm_2.6. Model yang terbaik dari regresi linear adalah model lm_2.6 dengan nilai r-square sebesar 73.09%. Model yang memiliki interaksi inilah yang memiliki nilai R-square yang paling tinggi diantara model tanpa dates. Terdapat perbedaan ketika model dengan variabel dates dimasukkan.
Variabel yang memiliki pengaruh siknifikan terhadap USD adalah stock index, consumer confidence, unemployment rate dan interaksi dari interest rate dan GDP growth. Variabel inilah yang akan digunakan dalam pembentukan model regresi yang lain.

lm_2.6$coefficients
##              (Intercept)              stock_index      consumer_confidence 
##              12709.13429                 10.59359                -47.01849 
##        unemployment_rate interest_rate:GDP_growth 
##               -472.17784                -37.60342

Model Regresi Terbaik :
USD = 12709.13429 - 47.01849 * consumer_confidence - 472.17784 * unemployment_rate + 10.59359 * stock_index - 37.60342 * interest_rate:GDP_growth

Dari hasil regresi tersebut didapatkan hasil bahwa Kurs USD dipengaruhi oleh consumer confidence, unemployment rate, stock index dan interaksi antara interest rate dan gdp growth. Nilai awal kurs USD atau nilai minimum ketika tidak ada faktor apapun adalah sebesar Rp12.709 . stock index(LQ45) memiliki dampak positif sebesar Rp10 untuk satu satuan. Unemployment rate, consumer confidence dan interaksi interesest rate dan gdp growth memiliki dampak negatif pada USD. Unemployment rate sebesar Rp472 , dan consumer confidence sebesar Rp47. Interaksi interest rate dan gdp growth memiliki dampak sebesar Rp37.

Decision Tree Regression

Kita akan membuat model regresi decision tree sesuai dengan variabel yang sudah ditentikan di regresi linear. Variabel yang akan digunakan adalah stock_index, consumer_confidence, unempolyment_rate, GDP_growth dan interest rate. Pada decision tree tidak bisa menggunakan interaksi, sehingga intraksinya akandipisah.

library(rpart)
fit2 <- rpart(USD ~ stock_index + consumer_confidence + unemployment_rate + 
    GDP_growth,
    method = "anova", data = dc[,-1]) 
dt_r2.2 <- cor(dc$USD, predict(fit2))^2
dt_r2.2
## [1] 0.8331815

Dari hasil model decision tree pertama, didapatkan r-square sebesar 0.8331 yang berarti 83.31% bariabel sudah masuk ke dalam model. Model decision tree sudah bisa dikatakan lebih bagus modelnya dibandingkan dengan regresi linear. Ssekarang, kita akan coba dengan ditambah interest rate.

fit3 <- rpart(USD ~ stock_index + consumer_confidence + unemployment_rate + 
    GDP_growth+ interest_rate,
    method = "anova", data = dc[,-1]) 
dt_r2.3 <- cor(dc$USD, predict(fit3))^2
dt_r2.3
## [1] 0.8476713
dt_r2.3-dt_r2.2
## [1] 0.01448981

Dari hasil model yang ditambahkan interest rate, r-squarenya adalah 0.8476. Ada kenaikan r-square sebesar 0.014489 atau 1.44%. Model yang lebih bagus untuk decition tree adalah model dengan tambahan interest rate.

MAPE.dt3 <- mean(abs(predict(fit3)- dc$USD) / dc$USD);MAPE.dt3 
## [1] 0.04974933

MAPE yang dihasilkan adalah 0.0497 atau 4.97%.

Random Forest

Random Forest adalah metode dengan membuat beberapa decision trees untuk memprediksi output lebih akurat. Bisa dibilang decision tree beberapa decision tree. Metode ini bisa menghindari overfitting.
Pada metode ini, kita akan gunakan random forest regression untuk memprediksi kurs US dollar. library yang akan digunakan adalah randomForest.

library(randomForest)
rf2 <- randomForest(USD~., data=dc[,-1],
                   importance=T)
rf2
## 
## Call:
##  randomForest(formula = USD ~ ., data = dc[, -1], importance = T) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##           Mean of squared residuals: 335600.8
##                     % Var explained: 92.94

hasil random forest mengatakan bahwa r-square dari model adalah 92.94%. mari kita coba untuk melihat dengan interaksinya. apakah r-squarenya meningkayt atau tidak.

rf3 <- randomForest(USD~stock_index+consumer_confidence+unemployment_rate+interest_rate:GDP_growth, data=dc,
                   importance=T)
rf3
## 
## Call:
##  randomForest(formula = USD ~ stock_index + consumer_confidence +      unemployment_rate + interest_rate:GDP_growth, data = dc,      importance = T) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##           Mean of squared residuals: 379333
##                     % Var explained: 92.02

Setelah variabel interaksi dimasukkan (interest rate dan gdp growth), hasil r-square random forest yang dihasilkan adalah 92.02%.

Model Evaluation


Mari kita bandingkan hasil regresi dan model time series.

linear_error <- data.frame(
  method = c("Linear Regression", "Decision Tree", " Random Forest"),
  rsq= c(summary(lm_2.6)$r.squared,dt_r2.3,0.9294)
)



linear_error
# MAPE regresi
test20 = c(round(length(dc$dates) *0.8) : length(dc$dates))
res_linear = abs(lm_2.6$residuals/dc$USD)
res_linear = mean(res_linear[test20])
MAPE.rf <- mean(abs(predict(rf2)- dc$USD) / dc$USD);MAPE.rf
## [1] 0.03289613
linear_error$MAPE <- c(res_linear, MAPE.dt3, MAPE.rf)*100
linear_error

Dari berdasarkan r-square dari Regression, model yang terbaik adalah model dengan random forest. Linear regresi sendiri punya 26%an variansi yang tidak masuk model. Model regresi linear sudah lumayan oke walaupun untuk r-square yang lumayan oke ada di angka 80%. Tetapi nilai MAPE yang dihasilkan lumayan besar, yaitu 3-6%an. Model yang akan kita coba bandingkan adalah randomforest dengan MAPE sebesar 3.2% dan R-square sebesar 0.92.

Lalu mari kita lihat untuk error model time series dengan menggunakan metode box-jenkins, ARCH, GARCH dan Holt’s.

data.table(acc_ts) # tabel nilai MAPE model time-series

Semakin besar nilai R-square, semakin banyak variabel tersebut masuk kedalam model. Model yang terbaik untuk time series adalah model ARIMA(3,1,1) dengan MAPE sebesar 1.40% dan r-square sebesar 0.517. R-square sebesar 51.7% yang berarti hanya 48.3% variansi yang tidak masuk ke dalam model. tabel di atas juga didapatkan model ARCH(1), GARCH(1,1) dan Holt’s tidak cocok, maka perwakilan dari model time series yang akan digunakan adalah ARIMA(3,1,1).

Berikut plot Prediksi dengan semua model time series.

h1= predict(holts,2)
h2=  predict(holts2,2)

h1= data.frame(h1)
h2 = predict(h2)
h11 = ts(c(fitted.values(holts), h1$Point.Forecast[1:30]    ), start = c(2010,1), frequency = 12) 
h22 = ts(c(fitted.values(holts2), h2$Point.Forecast[1:30]), start = c(2010,1), frequency = 12) 

p1= predict(dc_arima1,2)
p2= predict(dc_arima2,2)
p3= predict(dc_arima3,2)
p4= predict(dc_arima4,2)

arima5 <- arima(dc$USD, order = c(2,1,0))
p5= predict(arima5,2)

ts <- ts(c(dc$USD), start = c(2010,1), frequency = 12) 

a1 = ts(c(fitted.values(dc_arima1) ,p1$pred), start = c(2010,1), frequency = 12) 
a2 = ts(c(fitted.values(dc_arima2) ,p2$pred), start = c(2010,1), frequency = 12) 
a3 = ts(c(fitted.values(dc_arima3) ,p3$pred), start = c(2010,1), frequency = 12) 
a4 = ts(c(fitted.values(dc_arima4) ,p4$pred), start = c(2010,1), frequency = 12) 
a5 = ts(c(fitted.values(arima5),p5$pred), start = c(2010,1), frequency = 12) 


hchart(ts, type="line", name="Actual US Dollar", color = "black")%>%
  hc_add_series(h11, type="line", name= "Holt's b=0.0568", color="#F9A828")%>%
  hc_add_series(a1, type="line", name= "ARIMA(1,1,1)", color="#03AED2")%>%
  hc_add_series(a2, type="line", name= "ARIMA(2,1,1)", color="#68D2E8")%>%
  hc_add_series(a3, type="line", name= "ARIMA(3,1,1)", color="#FDDE55")%>%
  hc_add_series(a4, type="line", name= "ARIMA(4,1,1)", color="#615EFC")%>%
  hc_add_series(a5, type="line", name= "ARIMA(2,1,0)", color="#A3D8FF")%>%  
  hc_add_series(h22, type="line", name= "Holt's b=0.0451", color="#F29191")%>%
  hc_title(text="Prediction Rp to USD Time Series Model")

Dari hasil nilai MAPE time series dengan menggunakan split 80-20, 80% initial 20% test data. nilai hasil MPAE yang terbaik adalah dengan mengggunakan model Box-Jenkins ARIMA(3,1,1) dengan MAPE 1.408%. Data kurang cocok untuk menggunakan ARCH dan GARCH. Untuk ARCH GARCH error 7% lumayan besar. ### Predict and compare Time series and regressions. Mari kita coba prediksi nilai dollar di tahun 2022 bulan januari dan febuari

predict(dc_arima3,2)
## $pred
## Time Series:
## Start = 145 
## End = 146 
## Frequency = 1 
## [1] 14380.68 14370.20
## 
## $se
## Time Series:
## Start = 145 
## End = 146 
## Frequency = 1 
## [1] 237.0551 397.2527
ts <- ts(c(dc$USD), start = c(2010,1), 
                frequency = 12) 
ts2 <- ts(c(fitted.values(dc_arima3),14380.68,14370.20), start = c(2010,1), 
                frequency = 12) 
ts.rf <-  ts(c(predict(rf2)), start = c(2010,1), frequency = 12)


hchart(ts2, type="line", name="Predicted US Dollar", color = "#FF9A00")%>%
  hc_add_series(ts, type="line", name= "Actual US Dollar", color="#615EFC")%>%
  hc_title(text="Actual vs Prediction Rp to USD")%>%
  hc_add_series(ts.rf, type="line", name="Random Forest Regression", color="#A1DD70")

Dari hasil perbandingan time series di atas, model ARIMA bisa dikatakan sudah bisa mengikuti pola dari data asli nilai kurs USD. Prediksi model Randomforest bisa dikatakan kurang akurat jika dibandingkan dengan data asli dari kurs USD. Model yang baik adalah model ARIMA(3,1,1).

Policy Implications

Nilai dari kurs USD memiliki trend naik, yang pasti ekonomi kita bisa semakin baik atau menurun. Ketika Kurs USD naik, maka mata uang Rupiah melemah dan ekonomi tidak terlalu stabil. Untuk meningkatkan Rupiah (Rp), Indonesia perlu untuk meningkatkan rate keoptimisan masyarakat terhadap pemerintah. Lalu interaksi antara suku bunga dan PDB perlu ditambahkan agar dollar menurun. Stock index memiliki pengaruh yang positif terhadap kurs dollar. Semakin tinggi stock indeks, berarti semakin tinggi transaksi dan ekonomi berjalan. Walaupun LQ45 memiliki nilai positif, nilai yang dihasilkan juga lumayan ketika penutupan dari LQ45 ada di rentang 500-900an. Jika ekonomis semakin membalik, pasti juga kepercayaan pada pemerintah, tingkat pengangguran dan interaksi suku bunga dan PBD bisa menurunkan dari hasil stock indeks yang meningkat.

Sedangkan prediksi untuk USD sendiri jika diprediksi menyentuh Rp14.380. Trens rata-rata USD juga sementara ada kenaikan tetapi belum naik dan siknifikan. Dari tahun 2021-2022 sepertinya ekonomi Indonesia stabil, USD tidak memiliki kenaikan trend yang siknifikan dan selama 2 tahun itu, USD lumayan stabil. diabndingkan dengan tahun-tahun sebelumnya. Menurut saya pemerintah perlu mengatur masalah-masalah yang ada di bidang ekonomi (suku bunga dan PDB), meningkatkan kepercayaan masyarakat dengan membuat kebijakan-kebijakan yang menguntungkan masyarakat, mengatasi masalah pengangguran dan mengstabilkan stock index agar pasar sektor ini.