CODEWITHSUNDEEP
rtime-seriespredictionhidden-markov-models
A1010
A1010

Reputation: 360

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

Upvotes: 2

Views: 1867

Answers (2)

mspeek
mspeek

Reputation: 176

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.

Upvotes: 4

batlike
batlike

Reputation: 698

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.

Upvotes: 0

Related Questions