Question to a Non-markovian process model Suppose that a Piece of chocolate is chewed and as broken down it distributes in the mouth. The mouth is assumed as a 2-dimensional $N \times M$ grid with discrete cells $(i,j)$ in which one fragment of chocolate can stay. The probability for a chocolate Fragment jumping from Position $(i,j)$ to $(i',j')$ within one chew movement is denoted by $p_{i,j;i',j'}$. As an assumption, $p_{i,j;i',j'}$ obeyes an uniform Distribution and the entries are $\frac{1}{9}$ for $i=i',j=j'$ and for the $(i',j')$ that neighboring on $(i,j)$; all other entries are zero. But because it is a non-markovian process (the Person who chews can feel the chocolate in the mouth and recognizes where the chocolate is located), the Transition Matrix $p_{i,j;i',j'}$ changes after every chew movement.
The probability for jumping from $(i,j)$ to $(i',j')$ under the condition that $(i,j)$ was reached from $(i'',j'')$ is denoted by $p_{i'',j'';i,j;i',j'}$. It holds $p_{i'',j'';i,j;i'',j''}= \frac{1}{9}+q$ and $p_{i'',j'';i,j;i',j'}= \frac{1}{9}- \frac{q}{8}$ for $(i',j') \neq (i'',j'')$ neighboring $(i,j)$ and $p_{i'',j'';i,j;i',j'}=0$ otherwise. Here, $q$ is a non-markovian Parameter for the chewing process.
Higher-order "memory" can be neglected in this model. Therefore, the probability for a chocolate Fragment hopping from $(i,j)$ to $(i^{(n)},j^{(n)})$ within $n$ chewing cycles according to this model is given by ($I$ denotes all in-between grid places):
$p_{i,j;(i^{(n)},j^{(n)})} = \sum_I p_{i,j;i^{(1)},j^{(1)}}p_{i,j;i^{(1)},j^{(1)};i^{(2)},j^{(2)}}p_{i^{(1)},j^{(1)};i^{(2)},j^{(2)};i^{(3)},j^{(3)}}...p_{i^{(n-2)},j^{(n-2)};i^{(n-1)},j^{(n-1)};i^{(n)},j^{(n)}}$.
QUESTION: Is there a General formula for $p_{i,j;(i^{(n)},j^{(n)})}$ in dependence on $q$? Multiplying all These matrices together is very troublesome; how can I do this?
I'm stuck here in the Matrix multiplication; every hints will be greatly appreciated.
 A: This can actually be turned into a Markov process. Instead of $N\times M$ states, you have somewhat fewer than $9\times N\times M$ states. The extra factor represents one of eight directions, or staying put, encoding the small history you need. For example, the state $(0,0)$ can be reached by $\uparrow, \nwarrow, \leftarrow$ or $\star$ (staying put). So you represent this with the states $(\uparrow,0,0),(\nwarrow,0,0),(\leftarrow,0,0),(\star,0,0)$, any of which can transition to $(\star,0,0), (\rightarrow,0,1), (\searrow,1,1),$ or $(\downarrow,1,0)$.
Because of the edges, the number of states will be $9NM - 6(N+M) + 4$.
It wasn't clear to me what $q$ depends on, but probably no more detail is needed to incorporate it into the transition matrix, as it simply increases or decreases the probability of staying put at any position, given what happened immediately before. I have to admit however that if $MN$ is large I would probably program an equivalent simulation directly rather than create a giant sparse matrix.
