Rawr does not currently consider interaction between combustion and hot streak. What Rawr does is compute the average number of casts it takes to consume the 3 charges of combustion given a starting crit rate using dynamic programming and then sets the crit rate to 3 divided by this average number of casts needed to consume 3 charges. The actual duration of combustion is then determined based on the cast rate of the cycle.

Sorry for slipping off the radar for a week and a half. At the risk of necro, I'll throw in my experience with Markov chains, since you guys are spending too much effort re-inventing the wheel.

A markov chain is a stochastic process with the markov property. A stochastic process is a series of states with probability distributions describing the movement in time from one state to the next. The markov property means that the probability distribution of your next transition depends only on your current state, not on any other previous states. Requiring knowledge of a finite number of previous states can be modeled as a markov chain as well, by appending the previous state-information to your current state, although this makes your chain considerably more complex (needing the last k transitions in n possible states is simulated by having n^k states).

If you have a finite number of states, then you can describe the probability distribution of going from each state as a simple vector, avoiding all the Hilbert Space machinery and landing you nicely in the soft, cushy place of linear algebra. Each state has a vector describing possible future states, with each entry in the vector corresponding to the probability of transitioning to the corresponding state. The entries of the vector sum to 1. If you append these vectors to each other, then you get a transition matrix. If you stretch your head a little bit, you can see that if your current state is a probability distribution vector, then your future state distribution is given by a simple matrix multiplication.

If you remember eigendata from linear algebra, then you remember or can easily be convinced that the steady-state distribution, ie what you end up with after transitioning a bunch which is the same as repeatedly multiplying by the transition matrix, is the eigenvector corresponding to the largest eigenvalue of the transition matrix. Since all of our transition vectors are nice probability vectors that sum to 1, the largest eigenvalue is 1 (this is a special case of a more general result about upper and lower bounds of largest eigenvalues). This is nice, because it means that the problem is solveable in a finite amount of time; finding one eigenvalue given the other is easy but finding either from scratch cannot be done in exact arithmatic for matrices larger than 4x4, as a result of Abel's theorem (insolvability of quintic equations). Instead of an insolvable equation we get Ax=x, something linearish.

It turns out finite-precision mathematics on computers gives an even faster way of finding eigenvectors from eigenvalues: v = inv(A-cI)x for (almost)* any x. An explanation follows, but you can take my word for it. Just remember two things: a, eigenvalues can be scaled and be the same eigenvalue, so normalize it; b, using floating point math, not fixnum, for the love of god.
* Almost in the mathematical sense meaning that the probability of chosing an x where this is not the case is 0. Not actually impossible, but probability 0.

In exact math (A-cI) is noninvertible for eigenvalue c, but finite precision makes it huge rather than nonexistant (I recomend floating point, not fixnum, computations). Then we get some crazy things going on. Inv(A) has the same eigenvectors as A, with the eigenvalues being inverted. A-I has the same eigenvectors with eigenvalues reduced by 1. And repeatedly multiplying x by A approaches an eigenvector, with the speed of approach determined by the ratio of the largest and second-argest eigenvalues. So if your second-largest eigenvalue of A is 0.99, then the corresponding eigenvalue of inv(A-I) is 100; the one corresponding to 1 is the inverse of floating-point underflow, which means around 10^14 or so. Multiply anything by this matrix, and the eigenvector corresponding to 1 is now 12 orders of magnitude larger than anything from any other eigenvector, totally drowning it out. Multiply it twice and everything but your eigenvector is swamped well below the limits of your computer to tell the difference.

There are extensions of Markov chains to continuous time, and to infinite states, but I can't think of anything where we have to model those. We want caster state at times when we cast something, which is discrete time, and the various states are finite and very granular (if we hold gear constant). Combinatorics tends to give us a large number of states, but they're fairly regular since most effects are independent.

For the combusion problem, the matrix is pretty regular-looking. You have 40 (or less) states representing the various combinations of your two counters (several of which won't exist for various reasons). Any particular column has only two possible outcomes with non-zero probability: you crit or you don't, moving your crit-counter up by 1 and either decrementing or not decrementing your charge counter by 1, with distribution matching your crit chance. Your matrix will be mostly 0's with what few numbers you do have lying on two minor diagonals. For the sake of actually getting information out of this, there are 10 out-of-combution states that always cycle to themselves with probability 1 that correspond to how much crit you had when combutions ran out. Your steady state should be some distribution of these telling you the average effectiveness of combustion... I think.

If we want to use markov processes for combustion, don't we have to link the last state to the first with a transition ratio that represents combustion's cooldown? Then the steady state would give you your average crit chance for the entire fight when you have combustion talented based on your base crit chance.

What the steady-state gives you is a measure of how many more crits your combustion button gives you, and in what timeframe. You can average this out over the fight length by the standard technique of dividing by cooldown. What you describe would not be inaccurate, but it would be both unnecessary, and unduly sensitive to environmental influences, while the markov whatsit I described would be a pretty reliable measure.

If you don't have some kind of a loop, won't the steady state be 100% time spent with combustion not active? After all that's what actually happens if you only use combustion once and then keep DPSing to infinity...

You are misunderstanding what the model is measuring. The end result is what state you were in when you used the last combustion charge, ie how many spells it took to get the third crit. That gives a notion of how effective your combustion was. Trying to find DPS from via the markov chain is putting the carriage in front of the horse. You get the average number of casts from the chain, you figure the average number of bonus crits, you calculate the raw damage (*not* DPS) that this would net you when all your other cooldowns blown, and you divide by the cooldown to get the DPS. Everything else is just multiplication, the markov chain is the hard part.