# Time inhomogeneous Markov Chain in Mathematica

I would like to create a discrete 2-state Markov process, where the switching probabilities in the transition matrix vary with time.

I can currently do the following, which creates a process with fixed Transition matrix, and then simulates, and plots, a short time series:

model = DiscreteMarkovProcess[{1, 0}, {{1 - p1, p1}, {p2, 1 - p2}}];
ListLinePlot[RandomFunction[model /. {p1 -> 0.1, p2 -> 0.5}, {0, 100}]]


I would like the transition probabilities, $p1$ and $p2$ to vary over time according to sine waves. Specifically, I would like the resultant probabilities to vary according to, for example:

$p1 = 0.1[0.5 + 0.5sine(\frac{2\pi t}{10})]$

$p2 = 0.5[0.5 + 0.5sine(\frac{2\pi t}{10})]$

I have tried a number of ways of implementing this, although, to be honest none of them are getting particularly close to emulating the behaviour that I would like.

Does anyone have any idea how I could implement this?

Best,

Ben

• I haven't time to work this out myself but something involving FoldList and an Array of the relevant transition matrices should work. The main catch is that RandomFunction outputs TemporalData, so you might need to use ["Paths"] to strip it back to a plain list. – Verbeia Feb 12 '15 at 2:44
• @Verbeia - thanks for your message. I have implemented a version of what you said below. Let me know if it is what you had in mind! – ben18785 Feb 12 '15 at 14:18

Thanks to @Verbeia, I have managed to find an answer, although it is a bit inelegant, and involves construction of my own 'custom' Discrete Markov chain. Please see this below:

First select relevant probability from matrix:

SelectProb[aVPast_, TrM_] := Module[{tempProb}, tempProb = Switch[aVPast, {1, 0}, TrM[[1, 1]], {0, 1}, TrM[[2, 2]]]];


Then, probabilistically choose the next state in the 1st order Markov chain:

NextStep[aVPast_, TrM_] := Module[{aRand, aProb}, aProb = SelectProb[aVPast, TrM];
aRand = RandomReal[]; If[aRand >= aProb, {1, 1} - aVPast, aVPast]];


Now create a function which generates a sinusoidally-varying transition matrix:

TransitionSineMatrix[t_, pdw_, pwd_] := {{1 - pdw (0.5 + 0.5 Sin[2 π t/365]),
pdw (0.5 + 0.5 Sin[2 π t/365])}, {pwd (0.5 + 0.5 Sin[2 π t/365]), 1 - pwd (0.5 + 0.5 Sin[2 π t/365])}}


Put it all together in a function:

DiscreteMarkovProcessTimeInHomogeneous[aInitialState_, tMax_] := FoldList[NextStep, aInitialState,
Array[TransitionSineMatrix[#, 0.5, 0.5] &, tMax]]


However, I would be most appreciative if anyone comes up with a nicer option!

Best,

Ben

Here's an answer which will accept an arbitrary transition matrix, with possible time-dependence. It is similar in principle to Ben's answer in that it's a custom routine to simulate a Markov chain. I believe that this is unavoidable as long as DiscreteMarkovProcess does not allow explicit manipulation of the time parameter.

# Code

inhomogeneousdiscretemarkovprocess[initialtime_, finaltime_,
timestep_, initialstate_, transitionmatrix_, pathnumbers_] :=
Module[{unsortedpaths, currentstate = initialstate, currenttime = initialtime, t},
unsortedpaths = Reap[
Do[
Do[
t = currenttime;
Sow[{currenttime, currentstate}];
currentstate =
RandomChoice[
Normal[transitionmatrix[t][[currentstate]]] ->
Range[Length[transitionmatrix[t]]]];
currenttime = currenttime + timestep;
, {Floor[finaltime - initialtime + 1]}];
currenttime = initialtime; currentstate = initialstate;
, {pathnumbers}]
];
If[pathnumbers =!= 1,
Partition[unsortedpaths[[2, 1]],
Floor[finaltime - initialtime + 1]], unsortedpaths[[2, 1]]]
]


# Usage

For the discrete case, the function takes 6 parameters:

• initial time of the simulation
• final time
• time step
• initial state of the system (numbered 1 to n)
• the transition matrix
• the number of realisations of the simulation (if you want multiple trajectories)

The important thing to note is that the transition matrix should be defined as a function of t, so that it can be written naturally with the explicit time-dependence. So in this case, the transition matrix would be

transitionmatrix[t_] := {{1 - 0.1 (0.5 + 0.5*Sin[2*Pi*t/10]),
0.1 (0.5 + 0.5*Sin[2*Pi*t/10])}, {0.3 (0.5 + 0.5*Sin[2*Pi*t/10]),
1 - 0.3 (0.5 + 0.5*Sin[2*Pi*t/10])}}


Finally, the output is a list containing the couples {time,state} for easy plotting, with each trajectory in a separate list.

# Example

Here's the function simulating the system in the question:

inhomogeneousdiscretemarkovprocess[1, 100, 1, 1, transitionmatrix, 1]

(*output->{{1, 1},{2, 1},{3, 1},{4, 1},{5, 1},{6, 1},{7, 1},{8,1},
{9, 1},{10, 1}}*)


And a plot:

ListPlot@inhomogeneousdiscretemarkovprocess[1,100,1,1,transitionmatrix,1]


Multiple trajectories:

 ListPlot@inhomogeneousdiscretemarkovprocess[1,100,1,1,transitionmatrix,3]


# What about continuous-time Markov chains?

Well, for the time-homogeneous case, you could extend the first function without too much problem. Actually, there is an algorithm out there that does just this called the Stochastic Simulation Algorithm (SSA) or Gillespie algorithm (see original or quick guide).

However, for time-inhomogeneous systems, things become tricky because expected waiting times for the next transition events change as you wait for the next transition. For that, you would need more sophisticated algorithms (such as this one). But that's another story...