Grey-Markov method in R

Question

In R, I have loaded the built-in time series: AirPassengers and split it in train- and testdata like this:

rm(list = ls())
data = AirPassengers

traindata = ts(data[1:(0.75*length(data))], frequency = 12)
testdata = ts(data[((0.75*length(data))+1):length(data)], frequency = 12)

from here I want to estimate future values of a time series with the traindata using the Grey-Markov method. I know the Grey-Markov method consist of a Grey GM(1, 1) forecasting model followed by a Markov chain forecasting model refinement. But is there a function in R that performs this Grey-Markov method on its own, just like, for example, the auto.arima function?

Answer 1

No one has answered my question, so I thought I'll try to code a solution based on an example and found some minor different results but I don't think that my code is wrong. I hope that this method will help someone in the future. Also, I'd like to hear if someone finds errors in my code, than I'll change this.

The example I used is the example out of the following article:
Zhan-Li, M., & Jin-Hua, S. (2011b). Application of Grey-Markov Model in Forecasting Fire Accidents. Procedia Engineering, 11, 314–318. https://doi.org/10.1016/j.proeng.2011.04.663

# Clean workspace
rm(list = ls())

##### Grey GM(1,1) Forecasting Method #####

# The original data series
x0 = c(8000, 5261, 4657, 4728, 4325)

# The accumulated generating operation (AGO)
x_ago = c()
for(i in 1:length(x0)) {
  if (i==1) {
    x_ago[i] = x0[i]
  }
  else {
    x_ago[i] = x_ago[i-1] + x0[i]
  }
}

# Determine estimate for the developing coefficient (notation: a)
# Determine estimate for the grey input (notation: b)
# (a, b)^T = (B^T B)^-1 B^T Y
B = matrix(nrow = length(x_ago)-1, ncol = 2)
B[,2] = 1
for (i in 1:nrow(B)) {
  B[i,1] = -((x_ago[i+1]+x_ago[i])/2)
}

Y = matrix(nrow = length(x0)-1, ncol = 1)
for (i in 1:nrow(Y)) {
  Y[i,1] = x0[i+1]
}

ab = (solve(t(B)%*%B))%*%(t(B)%*%Y)
a = ab[1, 1]
b = ab[2, 1]

# Applying the inverse accumulated generating operation (IAGO)
# x_iago(k) = (1 - e^a)(x0(1) - b/a)e^-a(k-1)
x_iago = c()
for(i in 1:length(x0)) {
  if (i==1){
    x_iago[i] = x0[i]
  }
  else {
    x_iago[i] = (1-exp(a))*(x0[1]-(b/a))*exp(-a*(i-1))
  }
}

##### Markov Forecasting Method #####
# Calculating the relative error (notation: re)
# re(k) = ((x0(k) - x_iago(k))/x0(k))*100%
rel_error = c()
for(i in 1:length(x0)) {
  rel_error[i] = ((x0[i]-x_iago[i])/x0[i])*100
}

# We have n = 5 observations, so we create n - 1 = 4 states based on the relative error
# Highest re = 2.797% (round = 3), lowest re = -4.597% (round = -5). Difference high - low = 8%
# State_width = Difference high - low / n - 1 = 8%/4 = 2%.
state_width = (round(max(rel_error))-round(min(rel_error)))/(length(x0)-1)
# So state 1 = (-5, -3], state 2 = (-3, -1], state 3 = (-1, 1] and state 4 = (1, 3]
initial_state = c()
for(i in 1:length(x0)) {
  for (k in 1: length(x0)-1){
    if(rel_error[i]> round(min(rel_error)) && rel_error[i]<= round(min(rel_error))+(state_width*k)){
      initial_state[i] = k
      break
    }
  }
}

# In addition, it is necessary to go in transition in k-steps from the initial state
steps = c(0)
for(i in 2:length(x0)){
  steps[i] = i-1
}

# Determine whether we can get from the initial state to the transition state in k-steps
# Normally you determine P(1), P(2), ... , P(n) separately and grab one row from every P(k)
# To make one statematrix. In this example those steps are done together 
statematrix = matrix(0, nrow = length(x0), ncol = length(x0)-1)
statematrix[1,] = c(0, 0, 0, 0)
Pk = c()
Pk[1:length(x0)-1] = 0

for(i in 1:(length(x0)-1)){
  for(k in 1:(length(x0)-i)) {
    if (initial_state[k] == initial_state[i+1]) {
      Pk[initial_state[k+i]] = 1
    }
  }
  statematrix[i+1,] = Pk
  Pk = c()
  Pk[1:length(x0)-1] = 0
}

# To make it transparent, I make a dataframe with all vectors in it
df = data.frame(x0, x_iago, rel_error, initial_state, steps, statematrix)
df

# Now we take the column of the statematrix with the highest sum, which is column 3 with boundaries (-1, 1]
# To forecast the next x_iago-value, which is x_iago(6), we use the grey-function
x_iago[6] = (1-exp(a))*(x0[1]-(b/a))*exp(-a*(6-1))

# Now we add the Markov refinement for the forecast-value (notation: y)
# y(6) = x_iago(6)*(1+(delta_boundaries/2)*100%)
y = x_iago[6]*1.01
y