Is this a job for mixture of experts regression or semi-hidden markov models or something else? Data
I have several thousand timeseries each comprising around 365 data points. Browsing through a few of them, it looks like each timeseries consists of several regimes (different number f regimes per series). I think each regime could be adequately modelled using a GLM, i.e. within each regime, the data looks like it is either normally or poisson distributed around a mean which is itself a linear function of time. 
Question
I am looking for a technique that can automatically run across all thousand timeseries and identify the best start and end points for each regime, as well as fit a GLM to each regime. If there was no trend in each regime, I could do this using a semi-hidden markov model. I'm wondering if there is a technique that extends semi-hidden markov models so that the emission distribution accepts covariates? And if so, is there an R package that implements this? Or is this a mixture of experts regression? I've had a look at {mixtools} and it doesn't seem encouraging.
I've looked into Bai & Parson's work from econometrics, but this seems like it will only work if the regime can be modeled with a GLM with gaussian link, and I don't love that this procedure requires a 'run in' period within each regime - i'd like to be able to discover arbitrarily short regimes.
Software
I'm working in R, but can probably use Python to do this in a pinch. Worst case, MATLAB.
 A: I think the depmixS4 package may do the job for you.  Here's a simple example:
# Create an HMM with covariates
x1 <- rnorm(1000)
x2 <- rnorm(1000)

p11 <- 0.6
p22 <- 0.8

state <- rep(1,1000)
u <- runif(1000)
for (i in 2:1000) {
  if (state[i-1] == 1) {
    if (u[i] < p11) state[i] <- 1 else state[i] <- 2
  } else {
    if (u[i] < p22) state[i] <- 2 else state[i] <- 1
  }
}

y <- rnorm(1000)
y[state == 1] <- y[state == 1] + x1[state == 1]
y[state == 2] <- y[state == 2] + x2[state == 2]

# Estimate the simple model
s1 <- depmix(response=y~x1+x2, nstates=2, trstart=runif(4), ntimes=1000)
s2 <- fit(s1, emc=em.control(rand=FALSE))
summary(s2)

with output:
Transition matrix 
            toS1      toS2
fromS1 0.6888864 0.3111136
fromS2 0.2130215 0.7869785

Response parameters 
Resp 1 : gaussian 
    Re1.(Intercept)      Re1.x1      Re1.x2    Re1.sd
St1      0.03555102  0.89855381 -0.02027341 1.1037094
St2     -0.01743131 -0.01589832  1.00617375 0.9434336

which is reasonably accurate.  
You are able to specify a GLM family, which I have not done in the above example.  
I have to say that if your series are short, the different states are going to have to have quite different generating processes in order to be able to identify them as different states with any accuracy, and if there are a lot of states, e.g., 15 different ones observed over 365 data points, you are going to have a hard time getting accurate parameter estimates - harder than normal for a series of the given length, because of the added complexity of estimating the state transition points.
