2

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 :)

A1010
  • 360
  • 5
  • 18

2 Answers2

4

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.

Dharman
  • 30,962
  • 25
  • 85
  • 135
mspeek
  • 176
  • 5
  • That's great, thank you very much! When you call the `pred_by_r_state` you are basically using a very simple model. Is there a way to use something like an Arima model inside each of the states? – A1010 Dec 18 '20 at 15:58
0

HMMs routinely used, like the library you invoked, are often one-to-one. By one-to-one, what I mean is: you've noticed that the prediction sequence length is the same as the input (observation) length.

Figure by Andrej Karpathy

For one-to-many, you may want to try LSTMs, which is amenable to one[input]-to-many[output], many[input]-to-[many]output, etc. This should allow you to do some short-term forecasting. There are many answers (like this one) that give an intuitive illustration for how this might work. If you would prefer not to work with deep learning models, maybe you can take a look at Kalman filters.

batlike
  • 668
  • 1
  • 7
  • 19
  • 1
    I am not interested in working with Deep Learning models. My question was more focused on how to predict values in a time-series scenario using "classical" models – A1010 Dec 16 '20 at 08:48
  • 1
    I understand. HMMs may not be the best model for your forecasting. Kalyan filters may be a better fit. – batlike Dec 16 '20 at 09:54
  • Kalman filters are *also* completely unrelated to the question being asked – Taylor Jul 07 '23 at 15:58