R/markov_step.R

In ricky-kotecha/rkHMM:

markov_step <- function(aFit,
                        last_change = c(),
                        p_previous_state = c(),
                        density = NULL) {

  if (is.null(aFit$transition_matrix)) {
    stop("Missing transition matrix")
  }

  if (length(last_change) == 0) {
    last_change = aFit$regimeTimeSeries$Yield_diff[ length(aFit$regimeTimeSeries$Yield_diff)]
  }

  if (length(p_previous_state) == 0) {
    previousState1 = aFit$regimeTimeSeries$S1[length(aFit$regimeTimeSeries$S1) -1 ]
    previousState2 = aFit$regimeTimeSeries$S2[length(aFit$regimeTimeSeries$S2) -1 ]
  } else {
    previousState1 = p_previous_state[1]
    previousState2 = p_previous_state[2]
  }

  toState1FromState1 = previousState1 * aFit$transition_matrix[1,1]
  toState1FromState2 = previousState2 * aFit$transition_matrix[2,1]
  p_state1 = max(toState1FromState1, toState1FromState2)

  toState2FromState1 = previousState1 * aFit$transition_matrix[1,2]
  toState2FromState2 = previousState2 * aFit$transition_matrix[2,2]
  p_state2 = max(toState2FromState1, toState2FromState2)

  p_transition = append(p_state1, p_state2)

  if (is.null(density)) {
    parameters = dplyr::bind_rows(aFit$regime_density_parameters)
    parameters = dplyr::select(aFit$states, mean = coefficients, sd)
    density = dnorm(last_change, parameters$mean, parameters$sd)
  }

  p_weighted_density = as.numeric(p_transition) * as.numeric(density)

  p_state = p_weighted_density / sum(p_weighted_density)
  p_state
}