rtime-seriespredictionhidden-markov-models

How to predict out-of-sample observations with depmixS4 package in R?

I have a series of univariate data and I want to fit a Hidden Markov Model on it using the depmixS4 package on R. My final goal is to predict the next k observations (let's say k = 10) for the data series. I am not really interested in predicting the new state (that is important, but not my final goal), but I want to predict the next values for the data series.

It is a snippet of code:

# My series
data = rnorm(10000)
df_1_col = data.frame(data)
colnames(df_1_col) <- c('obs')

# Model
mod <- depmix(obs ~ 1, data = draws, nstates = n_state)
fit.mod <- fit(mod)

At this point I don't know how to predict the next out-of-samples values. I would like something similar to the forecast function in the forecast package.

I try using the following code:

state_ests <- posterior(fit.mod)
pred_resp <- matrix(0, ncol = n_state, nrow = 10)

for(i in 1:n_state) {
  pred_resp[,i] <- predict(fit.mod@response[[i]][[1]])
}

Using this code the predict function generates a number of predicted values that is equal to the number of observations into data, so it is not right.

How can I do this quite basic operations? I am new to HMM, but I already tried to look into many resources and I did not find any information. Thanks :)

Solution

  • A hidden Markov model models the observed variables conditional upon the hidden states. Forecasting the observed variables therefore requires an intermediate step of forecasting the hidden states. Once you have the forecasted probabilities of the hidden states, you can then forecast the observed variables from the marginal distribution of the observed variables, e.g.

    P(Y[T+k]|Y[1:T]) = \sum_i P(Y[T+k]|S[T+k] = i) * P(S[T+k] = i|Y[1:T])

    You can obtain the predicted state distribution by multiplying P(S[T]|Y[1:T]) and the state-transition matrix.

    library(depmixS4)

n_state <- 2

# My series
draws <- data.frame(obs=rnorm(10000))

# Model
mod <- depmix(obs ~ 1, data = draws, nstates = n_state, stationary=TRUE)
fit.mod <- fit(mod)

# extract the state-transition matrix
transition_mat <- rbind(getpars(getmodel(fit.mod,"transition",1)),getpars(getmodel(fit.mod,"transition",2)))

# extract the probability of the states at the final time point in the data (t=T)
# this will act as a "prior" to compute the forecasted state distributions
prior_vec <- as.numeric(posterior(fit.mod)[1000,-1])

# state-wise predictions for the observed variables
pred_r_by_state <- c(getpars(getmodel(fit.mod,"response",1))[1],
                     getpars(getmodel(fit.mod,"response",2))[1])

# for T + 1
# the forecasted state distribution is (prior_vec %*% transition_mat)
# so hence the prediction of the observed variable is
sum(pred_r_by_state * (prior_vec %*% transition_mat))

# for T + 2
# the forecasted state distribution is (prior_vec %*% transition_mat %*% transition_mat)
# so hence the prediction of the observed variable is
sum(pred_r_by_state * (prior_vec %*% transition_mat %*% transition_mat))

# for T + 3
sum(pred_r_by_state * (prior_vec %*% transition_mat %*% transition_mat %*% transition_mat))

# etc

    You will probably want to use the expm package, which contains the %^% operator, so you can use

    transition_mat %^% 3

    instead of

    transition_mat %*% transition_mat %*% transition_mat

    If the model contains covariates in the models of the observed predictors, you would also need to take these into account, i.e. try to somehow predict the values of these when computing pred_r_by_state.