Indonesian GDP Forecasting Using QRNN-MIDAS
Wisnowan Hendy Saputra
LIBRARY
library(readxl)
library(qrnn)
library(ggplot2)
library(forecast)
library(scales)
library(quantreg)
library(ggcorrplot)
library(neuralnet)DATA AND FUNCTION
This research applies the QRNN-MIDAS model to forecast Indonesia’s GDP data. Forecasting is based on GDP growth data which is modeled using the Financial Condition Index (FCI) predictor and includes the Covid-19 dummy variable. Complete data in this publication can be downloaded at the following link:
getwd()## [1] "D:/WORKSHOP/FORECAST PDB/SKENARIO 6"setwd("D:/WORKSHOP/FORECAST PDB/SKENARIO 6")
source("Function.R")
#Import training data
train0 <- read_excel("data6.xlsx",sheet = 1)
train <- data.frame(train0)
head(train0)## # A tibble: 6 × 10
## Kuartal GDP `GDP Growth` Bulan Rate Exc JCI CCI Inflasi Covid
## <chr> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2011Q1 1748731. 6.48 2011-01-01 6.5 9057 3409. 106. 0.89 0
## 2 2011Q2 1816268. 6.27 2011-02-01 6.75 8823 3470. 106. 0.13 0
## 3 2011Q3 1881850. 6.01 2011-03-01 6.75 8709 3679. 107. -0.32 0
## 4 2011Q4 1840786. 5.94 2011-04-01 6.75 8574 3820. 107. -0.31 0
## 5 2012Q1 1855580. 6.11 2011-05-01 6.75 8537 3837. 105. 0.12 0
## 6 2012Q2 1929019. 6.21 2011-06-01 6.75 8597 3889. 109 0.55 0#Import testing data
test0 <- read_excel("data6.xlsx", sheet = 2)
test <- data.frame(test0)
ytest <- test$GDP.Growth[1:4];ytest## [1] 5.044077 5.174244 4.936915 5.034001#Variable
tr.PDB <- tr.y <- na.omit(train[,2])
tr.Pertumbuhan_PDB <- na.omit(train[,3])
tr.Eksogen <- train[,5:9]CREATION OF EXOGENOUS PREDICTORS
Financial Condition Index (FCI) Variable
#Standardization of Exogenous Variables
tr.X <- scale(data.frame(tr.Eksogen))
#Principal Component Analysis
data.pca <- prcomp(tr.X)
summary(data.pca) #chosen 2 PC (>75%)## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.5076 1.1066 0.9412 0.7181 0.31824
## Proportion of Variance 0.4546 0.2449 0.1772 0.1031 0.02026
## Cumulative Proportion 0.4546 0.6995 0.8766 0.9797 1.00000#PC matrix rotation
rot.pca <- prcomp(scale(tr.X))$rotation
#Get FCI1 and FCI2 Variables
FCI1 <- FCI(1, tr.X)
FCI2 <- FCI(2, tr.X)
#Store FCI as x1 and x2
x1 <- matrix(rep(0,length(FCI1)), ncol = 3)
x2 <- matrix(rep(0,length(FCI2)), ncol = 3)
for (i in 1:(length(FCI1)/3)) {
x1[i,3] <- FCI1[(i-1)*3+1]
x1[i,2] <- FCI1[(i-1)*3+2]
x1[i,1] <- FCI1[(i-1)*3+3]
x2[i,3] <- FCI2[(i-1)*3+1]
x2[i,2] <- FCI2[(i-1)*3+2]
x2[i,1] <- FCI2[(i-1)*3+3]
}Covid-19 Dummy Variable
#Get Dummy Variable
X3 <- train$Covid
#Store Dummy as x3
x3 <- matrix(rep(0,length(X3)), ncol = 3)
for (i in 1:(length(FCI1)/3)) {
x3[i,3] <- X3[(i-1)*3+1]
x3[i,2] <- X3[(i-1)*3+2]
x3[i,1] <- X3[(i-1)*3+3]
}CREATION OF RESPONSE VARIABLE
y <- tr.Pertumbuhan_PDBTRAINING PREDICTION PROCESS
Initialization of Training Data Information
N <- length(y) # number of training sample
M <- 4 # maximum lag of the response variable (number of testing sample)Save training & testing data
#training
tr <- data.frame(Y = y[-c(1:M)],
X11 = x1[-c(45:48),1],
X12 = x1[-c(45:48),2],
X13 = x1[-c(45:48),3],
X21 = x2[-c(45:48),1],
X22 = x2[-c(45:48),2],
X23 = x2[-c(45:48),3],
x31 = x3[5:48,1],
x32 = x3[5:48,2],
x33 = x3[5:48,3],
l2 = y[-c(1:2,47:48)],
l4 = y[-c(45:48)])
#testing
te <- data.frame(X11 = x1[c(45:48),1],
X12 = x1[c(45:48),2],
X13 = x1[c(45:48),3],
X21 = x2[c(45:48),1],
X22 = x2[c(45:48),2],
X23 = x2[c(45:48),3],
x31 = test$Covid[c(1,4,7,10)],
x32 = test$Covid[c(2,5,8,11)],
x33 = test$Covid[c(3,6,9,12)],
l2 = c(y[c(47:48)],NA,NA),
l4 = y[c(45:48)])Initialize Experiment Settings (Hyperparameters)
cek.smape <- 200 # max of smape value
smape.final <- 200
kombinasi <- rep(1,5) # save the best model combination
tau <- seq(0.05,0.95,0.05) # 19 quantile values
elap.p <- matrix(rep(0,32),ncol = 2) # 16 parameter combinations of ELAP
# ELAP: Exponential Lag Almon Polynomial
elap.p[,1] <- rep(c(-0.2,-0.1,0.2,0.1),each=4)
elap.p[,2] <- rep(c(-0.2,-0.1,0.2,0.1),4)
eps <- 2^seq(-8,-16,-1) # 9 huber parameter values
neu <- c(6:12) # 7 number of nodes
maxit <- 10 # 10 initial value trials
# Best hyperparameter (after tunning)
i <- 9
j <- J <- 3
k <- 1
l <- 1
m <- 10Final QRNN-MIDAS Model
#Final Model
set.seed(m) # setting the initial value of hyperparameters
mo.final <-
qrnn.fit(x = as.matrix(cbind( # predictor Variables (Input)
x1 = elap(tr[,2:4],elap.p[j,1],elap.p[j,2]), # exogenous high frequency predictor (X1)
x2 = elap(tr[,5:7],elap.p[j,1],elap.p[j,2]), # exogenous high frequency predictor (X2)
x3 = elap(tr[,8:10],elap.p[j,1],elap.p[j,2]), # exogenous high frequency predictor (X3)
l2 = tr[,11], # endogenous low frequency predictor (2nd lag)
l4 = tr[,12])), # endogenous low frequency predictor (4th lag)
y = as.matrix(tr[,1]), # response variable (Target)
n.hidden = neu[l], # architectural scheme in hidden layer (number of nodes)
tau = tau[i], # quantile value (tau)
iter.max = 1000, # number of iterations (quasi-newton)
n.trials = 1, # number of trial modeling Neural Network (NN)
Th = sigmoid, # activation function
eps.seq = eps[k]) # huber approximation parameter values## tau = 0.45
## 1/1
## 1 0.005222152
## * 0.005222152Training Data Predictions
j <- J
pred.final.tr <- qrnn.predict(
x = as.matrix(cbind(x1 = elap(tr[,2:4],elap.p[j,1],elap.p[j,2]),
x2 = elap(tr[,5:7],elap.p[j,1],elap.p[j,2]),
x3 = elap(tr[,8:10],elap.p[j,1],elap.p[j,2]),
l2 = tr[,11],
l4 = tr[,12])),
parms = mo.final)[,1]FORECAST TESTING DATA
N2 <- 4 # predicted number of test samples
tes <- te
ypred <- y[-c(1:4)]
pred.final.te <- rep(0,N2)
for (i in 1:N2) {
j <- J
if (is.na(tes[i,10])) {
tes[i,10] <- ypred[length(ypred)]
pred.final.te[i] <- qrnn.predict(
x = as.matrix(cbind(x1 = elap(tes[i,1:3],elap.p[j,1],elap.p[j,2]),
x2 = elap(tes[i,4:6],elap.p[j,1],elap.p[j,2]),
x3 = elap(tes[i,7:9],elap.p[j,1],elap.p[j,2]),
l2 = tes[i,10],
l4 = tes[i,11])),
parms = mo.final)[,1]
ypred <- c(ypred,pred.final.te[i])
}else{
pred.final.te[i] <- qrnn.predict(
x = as.matrix(cbind(x1 = elap(tes[i,1:3],elap.p[j,1],elap.p[j,2]),
x2 = elap(tes[i,4:6],elap.p[j,1],elap.p[j,2]),
x3 = elap(tes[i,7:9],elap.p[j,1],elap.p[j,2]),
l2 = tes[i,10],
l4 = tes[i,11])),
parms = mo.final)[,1]
ypred <- c(ypred,pred.final.te[i])
}
}
PREDICTION <- rep(0,N-4)
for (i in 1:length(PREDICTION)) {
PREDICTION[i] <- pred.final.tr[i]/100*tr.PDB[i]+tr.PDB[i]
}
te.PDB <- c(tail(tr.PDB,4))
FORECAST <- rep(0,N2)
for (i in 1:length(FORECAST)) {
FORECAST[i] <- pred.final.te[i]/100*te.PDB[i]+te.PDB[i]
}FORECAST FUTURE SAMPLES (Q1-Q4 2024)
Initialization of Forecast Data Information
N3 <- 4
te.Eksogen <- test[,5:9]Creation of FCI and dummy (forecast)
# standarization
te.X <- scale(data.frame(te.Eksogen), center = sapply(tr.Eksogen, mean),
scale = sapply(tr.Eksogen,sd))
# get FCI1, FCI2, and dummy
FCI1 <- FCI(1, te.X)
FCI2 <- FCI(2, te.X)
X3f <- test$Covid
# save FCI as x1 and x2
x1 <- matrix(rep(0,length(FCI1)), ncol = 3)
x2 <- matrix(rep(0,length(FCI2)), ncol = 3)
x3 <- matrix(rep(0,length(X3f)), ncol = 3)
for (i in 1:(length(FCI1)/3)) {
x1[i,3] <- FCI1[(i-1)*3+1]
x1[i,2] <- FCI1[(i-1)*3+2]
x1[i,1] <- FCI1[(i-1)*3+3]
x2[i,3] <- FCI2[(i-1)*3+1]
x2[i,2] <- FCI2[(i-1)*3+2]
x2[i,1] <- FCI2[(i-1)*3+3]
x3[i,3] <- X3f[(i-1)*3+1]
x3[i,2] <- X3f[(i-1)*3+2]
x3[i,1] <- X3f[(i-1)*3+3]
}
y <- na.omit(test$GDP.Growth)
# save training dataset
fore <- data.frame(X11 = x1[,1],
X12 = x1[,2],
X13 = x1[,3],
X21 = x2[,1],
X22 = x2[,2],
X23 = x2[,3],
X31 = x3[,1],
X32 = x3[,2],
X33 = x3[,3],
l2 = c(y[c(3:4)],NA,NA),
l4 = y[c(1:4)])
# Forecasting Calculation
yfore <- c(0)
fore.final <- rep(0,N3)
for (i in 1:N3) {
j <- J
if (is.na(fore[i,10])) {
fore[i,10] <- yfore[length(yfore)]
fore.final[i] <- qrnn.predict(
x = as.matrix(cbind(x1 = elap(fore[i,1:3],elap.p[j,1],elap.p[j,2]),
x2 = elap(fore[i,4:6],elap.p[j,1],elap.p[j,2]),
x3 = elap(fore[i,7:9],elap.p[j,1],elap.p[j,2]),
l2 = fore[i,10],
l4 = fore[i,11])),
parms = mo.final)[,1]
yfore <- c(yfore,fore.final[i])
}else{
fore.final[i] <- qrnn.predict(
x = as.matrix(cbind(x1 = elap(fore[i,1:3],elap.p[j,1],elap.p[j,2]),
x2 = elap(fore[i,4:6],elap.p[j,1],elap.p[j,2]),
x3 = elap(fore[i,7:9],elap.p[j,1],elap.p[j,2]),
l2 = fore[i,10],
l4 = fore[i,11])),
parms = mo.final)[,1]
yfore <- c(yfore,fore.final[i])
}
}
fore## X11 X12 X13 X21 X22 X23 X31 X32 X33
## 1 1.606750 1.690406 1.571104 -0.7933052 -0.7149091 -0.9687749 0 0 0
## 2 1.529466 1.510744 1.571933 -0.9680019 -0.9844955 -1.1868838 1 0 0
## 3 1.813351 1.800403 1.688060 -0.7211626 -0.7197450 -0.8709680 1 1 1
## 4 1.865708 1.739057 1.723182 -1.2344174 -1.1413720 -0.8568368 1 1 1
## l2 l4
## 1 4.936915 5.044077
## 2 5.034001 5.174244
## 3 4.112099 4.936915
## 4 4.988000 5.034001Forecast <- rep(0,N3)
for (i in 1:length(Forecast)) {
Forecast[i] <- fore.final[i]/100*test$GDP[i]+test$GDP[i]
}FINAL ACCURACY
MAPE.PREDICTION <- round(accuracy(PREDICTION, tr.PDB[-c(1:4)])[5],4)
MAPE.FORECAST <- round(accuracy(FORECAST, te.PDB)[5],4)
data.frame(Sample = c("Training", "Testing"),
`Accuracy (MAPE)` = c(paste(MAPE.PREDICTION,"%"),paste(MAPE.FORECAST,"%")))## Sample Accuracy..MAPE.
## 1 Training 0.0282 %
## 2 Testing 5.0688 %FINAL VISUALIZATION
Actual2 <- c(tr.PDB,test$GDP[1:4],rep(NA,4))
Training.Prediction <- c(rep(NA,4),PREDICTION,FORECAST[1],rep(NA,7))
Testing.Prediction <- c(rep(NA,N),FORECAST,Forecast[1],rep(NA,3))
Forecast2 <- c(rep(NA,N+N2),Forecast)
fore.Bulan <- as.character(seq.Date(from = as.Date("2024-01-01"),
to = as.Date("2024-12-01"), by = "1 months"))
Dates <- c(train$Bulan,test$Bulan,fore.Bulan)
dfplot <- data.frame(Q = Dates[1:56*3],
Actual = Actual2,
Training.Prediction = Training.Prediction,
Testing.Prediction = Testing.Prediction,
Forecast = Forecast2)
colorplot <- c("Actual" = "black",
"Training.Prediction" = "red",
"Testing.Prediction" = "blue",
"Forecast" = "brown")
g1 <- ggplot(dfplot,aes(x = as.Date(Q)))+
geom_line(aes(y = Actual,color="Actual"))+
geom_line(aes(y = Training.Prediction,color="Training.Prediction"))+
geom_line(aes(y = Testing.Prediction,color="Testing.Prediction"))+
geom_line(aes(y = Forecast,color="Forecast"))+
labs(x = "Time (Quarter)",y = "Gross Domestic Product (Billion Rupiah)",color = "")+
scale_color_manual(values = colorplot)+
scale_y_continuous(labels = dollar_format(prefix="Rp"))+
scale_x_date(labels = function(x) zoo::format.yearqtr(x, "%YQ%q"),date_breaks="year")+
theme(
plot.title = element_text(color = "#0099f9", size = 20, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 13, face = "bold", hjust = 0.5),
plot.caption = element_text(face = "italic",hjust = 0),
axis.text.x = element_text(angle=45),legend.position = "bottom")+
geom_rect(aes(xmin = as.Date("2024-03-01"),xmax = as.Date("2024-12-01"),
ymin = -Inf, ymax = Inf),fill="yellow",alpha = .005)
g1
cat("Future Forecast for (Q1-Q4 2024)","\n",na.omit(dfplot$Forecast))## Future Forecast for (Q1-Q4 2024)
## 3111292 3202256 3280868 3294147