The purpose of this post is to give a workflow that one could follow if your goal is to fit Non-Homogeneous Hidden Markov Models (NHMMs) and to simulate sequences of states. Julie Quinn wrote a series of great posts on fitting Hidden Markov Models (HMMs) here but the goal of this post is to discuss the non-homogeneous counterparts of HMM. NHMMs can be distinguished from the former because they involve non-stationary transition probabilities between states. These dynamic transition probability matrices are conditioned on one or more external covariates that influences transitions between states. One example of the use of NHMMs is to model and simulate precipitation. Precipitation models that simulate without using atmospheric information cannot be expected to perform well on conditions that deviate from those on which the model was fit. Thus using a covariate that shows changes in atmospheric circulation (such as geopotential height), can help capture some of the nonstationarity that a solely precipitation-based model alone could not capture.

I have had a lot of time to work with NHMMs in R; however at the beginning, I found overwhelmingly few examples to work off of, so this post is meant to outline a workflow that I have settled on that you can apply to your own application. First off, there are tons of packages available for HMMs, but very few that handle the dynamic transition matrices of NHMMs. My favorite is depmixS4 found here. This package has all the tools needed to fit your NHMM, but it takes some experimenting to understand the syntax and to piece together the functions to create a full workflow.

First we want to fit a depmix model. In the first line of the code, I am creating a model called *modNHMM*. The data that I am fitting the NHMM with are the first 9 principal components of geopotential height. It is important that you list these components in this syntax and they should be the column titles of your dataset, which is *synoptic.pcs* in my case. The number of states is how many states you want to fit your model with. For a precipitation/streamflow application, you could fit a two-state NHMM to represent wet/dry conditions, or if you are trying to identify multiple regimes, you may have more.

```
library(depmixs4)
modNHMM <- depmix(list(PC1~1,PC2~1,PC3~1, PC4~1, PC5~1, PC6~1, PC7~1, PC8~1,
PC9~1),
nstates = num.states,
family=list(gaussian(),gaussian(),gaussian(),gaussian(),gaussian(),gaussian(),gaussian(),gaussian(),
gaussian()),
ntimes = nrow(synoptic.pcs),
data = data.frame(synoptic.pcs),transition = ~predictor$PC1+predictor$PC2+predictor$PC3+predictor$PC4)
fit.modNHMM.depmix.paleo.check <- fit(modNHMM)
synoptic.state.assignments<- posterior(fit.modHMMs.depmix.paleo.check)$state # state sequence using the Viterbi algorithm
```

Then we choose a distribution for our responses, which I choose to be a Gaussian distribution. The argument, *ntimes*, refers to the length of our time series, which is the number of rows in our dataset. Next we specify the dataset that contains the PC data and the transition which is dictated by our external covariates. In this case, my covariates are in a dataframe called *predictor* and each column corresponds to the 4 covariates (which happen to be PCs of a different dataset) that will influence my transitions. Then we fit the model with our prescribed attributes using the *fit* function. Finally, we want to calculate the Viterbi sequence, which is the most probable path or sequence of states over the period that the model is fit. This step (last line of code) will return a vector of the length of the dataset with each day classified into one of *num.states*.

Now we need to locate the transition probability matrices. We use the following command:

`fit.modNHMM.depmix.paleo.check@transition`

If we were fitting an HMM, we would get one transition probability matrix. However, we get an output that looks like this:

Now instead of having constant values for transitions, we have equations of the form: *Intercept+b1*PC1+b2*PC2+b3*PC3+b4*PC4* where the *b* values are the coefficient values listed in the table. If were were to look at the first block [[1]], this block dictates transitions from State 1 to the other 5 states. The transitions from the states are as follows:

State 1 to State 1: 0+0+0+0+0 (State 1 is used as a reference variable in this case and the probability would be found by subtracting the other probabilities from 1 at the end)

State 1 to State 2: -3.11+-0.22*PC1+0.22*PC2+0.014*PC3-0.13*PC4

and so on. You have to remember to divide each value by the total across the row in order to return probabilities.

Once you have these symbolic equations for the transition probability matrices, you can create a list of matrices which will allow you to simulate sequences of states for new sets of the PC1, PC2, PC3, PC4 covariates. You can get a sense how you might create *n* different transition matrices if you have times series of length *n* of the covariates. Below I am just representing those symbolic equations in code using the *getpars* function to acquire the coefficients and store the resulting daily matrices in a list called *mm*. Depending on the number of covariates or states, you will need to adjust the indices accordingly.

```
n=dim(df)[[1]]
mm<-matrix(list(), 1,n)
for (j in 1:n){
transition_matrix=matrix(, nrow = 5, ncol = 5)
for (i in 6:10){
transition_matrix[1,i-5]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[1,]))
for (i in 6:10){
transition_matrix[1,i-5]=exp(transition_matrix[1,i-5])/denominator
}
for (i in 31:35){
transition_matrix[2,i-30]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[2,]))
for (i in 31:35){
transition_matrix[2,i-30]=exp(transition_matrix[2,i-30])/denominator
}
for (i in 56:60){
transition_matrix[3,i-55]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[3,]))
for (i in 56:60){
transition_matrix[3,i-55]=exp(transition_matrix[3,i-55])/denominator
}
for (i in 81:85){
transition_matrix[4,i-80]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[4,]))
for (i in 81:85){
transition_matrix[4,i-80]=exp(transition_matrix[4,i-80])/denominator
}
for (i in 106:110){
transition_matrix[5,i-105]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[5,]))
for (i in 106:110){
transition_matrix[5,i-105]=exp(transition_matrix[5,i-105])/denominator
}
mm[[j]]=transition_matrix
}
```

Once we have these matrices, we can then simulate state sequences that can result from the chain of transition matrices. For this part, we need to create markov lists with our transition matrices:

```
library(markovchain)
mcObject <- mclapply(X=1:iter,mc.preschedule=TRUE,mc.cores=1,FUN=function(j){
mcObject.time.varying <- mclapply(X=1:n.sim,mc.preschedule=TRUE,mc.cores=1,FUN=function(t){
tr.prob=as.matrix(mm[[t]])
mcObject.time.varying.out <- new("markovchain", states = c("1","2","3","4","5"),
transitionMatrix = tr.prob, name = paste("McObject",t,sep=""))
return(McObject.time.varying.out)
}
)
mcObject.final <- new("markovchainList",markovchains = mcObject.time.varying, name = "mcObject.nh")
return(
mcObject.final
)
}
```

Finally we simulate using the following:

```
simulate.mc <- function(mcObject,num.states,dates.sim,last.month,last.day,n.sim,iter) {
#this function will simulate the Markov chain iter times
#Arguments:
#mcObject = a Markov chain object from the markovchain package
#num.states = the number of states
#date.sim = a time series of dates for the simulation period
#last.month = last month of the season
#last.day = last day of the last month of the season
#iter = the number of iterations
#day and month sequences for the simulation period
days.sim <- as.numeric(format(dates.sim,"%d"))
months.sim <- as.numeric(format(dates.sim,"%m"))
n.sim <- length(dates.sim) #length of simulation
final.mc.sim <- mclapply(X=1:iter,mc.preschedule=TRUE,mc.cores=1,FUN=function(i){
mc.sim <- as.numeric(rmarkovchain(n=1,object=mcObject[[i]][[1]]))
end.state <- paste(mc.sim[length(mc.sim)])
for (t in 1:n.sim) {
mc.sim <- c(mc.sim,as.numeric(rmarkovchain(n=1,object=mcObject[[i]][[t]],t0=end.state)))
end.state <- paste(mc.sim[length(mc.sim)])
if(months.sim[t]==last.month & days.sim[t]==last.day) {end.state <- paste(sample(1:num.states,size=1))}
}
#here is the final mc simulation
final.mc.sim.iter <- mc.sim[2:(n.sim+1)]
return(final.mc.sim.iter)
}
)
return(final.mc.sim)
}
simulations=matrix(list(), 1,1000)
for (i in 1:1000){
simulations[[i]] <- simulate.mc(mcObject=mcWeather.Regime,num.states=num.states,
dates.sim=dates.sim,last.month=last.month,last.day=last.day,iter=iter)
}
```

And that’s it! You can simulate for many different iterations (1000 in my case) and you will be returned a large list with your 1000 sequence of states over the simulation period.