A Novel Optimized Nonlinear Grey Bernoulli Model for Forecasting China’s GDP

Algorithm 1

	#R program of simulation and prediction in the original NGBM(1,1) model
	raw_data <- c(53.8580, 59.2963, 64.1280, 68.5992, 74.0060, 85.0754, 90.0309)
	x <- raw_data [1:5] #training data
	y <- x[−1] #Y
	m <- length (x)
	ago <- cumsum(x) #1-AGO sequence
	B <- matrix(nrow = m − 1, ncol = 2)
	n <- seq(−1, 0.99, 0.01) #Given range of parameter n
	#compute all errors from different n
	err <- NULL
	for (i in 1:length(n)){
	for (j in 1:dim(B) [1]){
	B[j, 1] = −0.5 ∗ (ago[j] + ago[j + 1])
	B[j, 2] = (0.5 ∗ (ago[j] + ago[j + 1])) ^ n[i]
	}#construct matrix B
	a <- (solve(t(B)% ∗ %B)% ∗ %t(B)% ∗ %y) [1].
	b <- (solve(t(B)% ∗ %B)% ∗ %t(B)% ∗ %y) [2].
	x_simu <- NULL
	for (k in 1:m){
	x_simu[k] = (b/a + (x [1] ^ (1 − n[i]) − b/a) ∗ exp(−(1 − n[i]) ∗
	a ∗ (k − 1))) ^ (1/(1 − n[i]))
	}
	#Compute the sum of absolute percentage errors
	err[i] <- sum(abs(diff(x_simu) − x[−1])/x[−1])
	}
	n = n[which.min(err)] #Optimal parameter n
	for (i in 1:dim(B) [1]){
	B[i, 1] <- −0.5 ∗ (ago[j] + ago[j + 1])
	B[i, 2] <- (0.5 ∗ (ago[j] + ago[j + 1]))^n
	}
	a <- (solve(t(B)% ∗ %B)% ∗ %t(B)% ∗ %y) [1] #parameter a
	b <- (solve(t(B)% ∗ %B)% ∗ %t(B)% ∗ %y) [2]#parameter b
	for (k in 1:(m + 2)){
	x_simu[k] = (b/a + (x[1] ^ (1 − n) − b/a) ∗ exp(−(1 − n) ∗ a ∗ (k − 1))) ^ (1/(1 − n))
	}
	x_pre <- diff(x_simu) #predicted values
	#compute error
	RPE <- (x_pre-raw_data[−1])/raw_data[−1]
	MAPE <- sum(abs(x_pre-raw_data[−1])/raw_data[−1])/(length(raw_data) − 1)
	#R program of the optimized NGBM(1,1) model
	raw_data <- c(53.8580, 59.2963, 64.1280, 68.5992, 74.0060, 85.0754, 90.0309)
	x <- raw_data[1:5] #training data & first rolling
	#x <- raw_data[2:6] #training data & second rolling
	y <- x[−1] #Y
	m <- length(x)
	ago <- cumsum(x) #1-AGO sequence
	B <- matrix(nrow = m − 1, ncol = 2)
	n <- seq(−1, 0.99, 0.01) #Given range of parameter n
	alpha <- seq(0, 1, 0.01) #Given range of parameter alpha
	#compute all errors from different n and alpha
	matri_err < - list()
	for (i in 1:length(n)){
	err <- NULL
	for (j in 1:length(alpha)){
	for (k in 1:dim(B) [1]){
	B[k,1] = −alpha[j] ∗ (ago[k] + ago[k + 1])
	B[k,2] = (alpha[j] ∗ (ago[k] + ago[k + 1]))^n[i]
	}#construct matrix B
	a <- (solve(t(B)% ∗ %B)% ∗ %t(B)% ∗ %y) [1].
	b <- (solve(t(B)% ∗ %B)% ∗ %t(B)% ∗ %y) [2].
	x_simu <- NULL
	for (l in 1:m){
	x_simu[l] = (b/a + (s[m] ^ (1 − n[i]) −b/a) ∗ exp(−(1 − n[i])∗
	a ∗ (l − m))) ^ (1/(1 − n[i]))
	}
	#Compute the sum of absolute percentage errors
	err[j] <- sum(abs(diff(x_simu) − x[−1])/x[−1])
	}
	matri_err[[i]] <- err
	}
	matri <- do.call(cbind, matri_err)
	dimo <- arrayInd(sort.list(matri) [1],dim(matri))
	n = r[dimo [1]] #Optimal parameter n
	alpha = alpha[dimo [2]] #Optimal parameter alpha
	for (k in 1:dim(B) [1]){
	B[k, 1] = −alpha ∗ (ago[k] + ago[k + 1])
	B[k, 2] = (alpha ∗ (ago[k] + ago[k + 1]))^n
	}#construct matrix B
	a <- (solve(t(B)% ∗ %B)% ∗ %t(B)% ∗ %y) [1]#parameter a
	b <- (solve(t(B)% ∗ %B)% ∗ %t(B)% ∗ %y) [2]#parameter b
	for (l in 1:(m + 2)){
	x_simu[l] = (b/a + (ago[m]^(1 − n) − b/a) ∗ exp(−(1 − n) ∗ a ∗ (l − m))) (1/(1 − n))
	}
	x_pre <- diff(x_simu) #predicted values
	#compute first rolling error
	RPE <- (x_pre-raw_data[−1])/raw_data[−1]
	MAPE <- sum(abs(x_pre-raw_data[−1])/raw_data[−1])/(length(raw_data) −1)