# Simulating a Markov chain

382 views (last 30 days)
John on 4 Jan 2013
Commented: Nazer Hdaifeh on 27 Aug 2020
Hello,
Would anybody be able to help me simulate a discrete time markov chain in Matlab?
I have a transition probability matrix with 100 states (100x100) and I'd like to simulate 1000 steps with the initial state as 1.
I'd appreciate any help as I've been trying to do this myself all week with no success and have not been able to source suitable code online.
Kind Regards
John

John D'Errico on 20 Apr 2019
Edited: John D'Errico on 21 Apr 2019
There seems to be many followup questions, it may be worth discussing the problem in some depth, how you might attack it in MATLAB. So lets start out with a discussion of such a Markov process, and how we would work with it. First, create a simple markov process. I'm not feeling terribly creative right now, so lets just pick something simple, thus a 5x5 transition matrix.
T = triu(rand(5,5),-1);
T = T./sum(T,2)
T =
0.17362 0.0029508 0.33788 0.19802 0.28752
0.059036 0.16812 0.29036 0.36644 0.11604
0 0.15184 0.30054 0.275 0.27262
0 0 0.42829 0.30672 0.26499
0 0 0 0.28417 0.71583
I've arbitrarily chosen a matrix with no absorbing state, but one where at each step, the process will tend to push the state to a higher number, but some chance you can climb backwards too.
If we look at the matrix above, if you are in state 5, with probability 0.71583 you will stay in state 5, but 28% of the time, you will drop back to state 4, etc. Next, consider a vector that describes your current state. Suppose we start out in state 1 (thus initially, 100% of the time, we are in state 1.)
X = [1 0 0 0 0];
After one step of the process, we don't know what state we will be in, but we can determine the probability that we will lie in any given state. That is found simply using a matrix multiplication.
% After one time step
X = X*T
X =
0.17362 0.0029508 0.33788 0.19802 0.28752
% After two time steps
>> X = X*T
X =
0.030319 0.052314 0.24588 0.27082 0.40066
% After three time steps
>> X = X*T
X =
0.0083525 0.04622 0.21532 0.28971 0.44039
% After four time steps
>> X = X*T
X =
0.0041788 0.040491 0.20504 0.29181 0.45848
Gradually, the probability grows that we will lie in state 5 MOST of the time. So after 100 time steps, the probability keeps growing that we lie in state 5.
X100 = [1 0 0 0 0]*T^100
X100 =
0.0025378 0.035524 0.19457 0.29167 0.4757
Does a steady state prediction of the long term state of this process exist? Yes. We can get that from the eigenvector of T', corresponding to the eigenvalue 1. That would be the vector V1, such that V1*T = V1.
eig(T')
ans =
1 + 0i
0.48873 + 0i
0.16513 + 0i
0.0054905 + 0.13467i
0.0054905 - 0.13467i
[V,D] = eig(T');
V1 = V(:,1)';
V1 = V1/sum(V1)
V1 =
0.0025378 0.035524 0.19457 0.29167 0.4757
That last step was to normalize the eigenvector to sum to 1, so they could be viewed as probabilities. Remember that eig normalizes its vectors to have unit norm.
So over the long term, the probability is 47.57% that the process lies in state 5.
This is what we can learn about the long term behavior of that system. But how about simulating the process? So, instead of thinking about where we will be as this process goes to infinity, can we simulate a SINGLE instance of such a Markov chain? This is a very different thing, since it does not rely on eigenvalues, matrix multiplication, etc.
Now, we need to look at the rows of T. I'll build this as a random walk that runs for many time steps. At any time, I'll simulate the progress of the random walk as it proceeds from one state to the next state. Then I'll simulate 1000 such random walks, and see where we ended at the final step of that process.
The trick is to use the cumulative sum of T, along the rows of T. What does that do for us?
CT = cumsum(T,2)
CT =
0.17362 0.17657 0.51445 0.71248 1
0.059036 0.22716 0.51752 0.88396 1
0 0.15184 0.45239 0.72738 1
0 0 0.42829 0.73501 1
0 0 0 0.28417 1
Suppose we start out in state 1. The first row of CT is pertinent here. Generate a single random number (using rand). If that number is less than 0.17362, then we started in state 1, and will remain in state 1. If the number lies between 0.51445 and 0.71248, then we will have moved to state 4, etc. We should see how to simulate this process now. I'll simulate 10000 such random Markov processes, each for 1000 steps. I'll record the final state of each of those parallel simulations. (If I had the parallel tolbox, I suppose I could do this using a parfor loop, but not gonna happen for me.)
I'll use histcounts, but older MATLAB releases need to use histc. I'll append a column of zeros at the beginning of CT to make histcounts return the proper bin index.
CT = [zeros(size(T,1),1),cumsum(T,2)];
numberOfSims = 10000;
simTime = 1000;
initialState = 1;
currentState = repmat(initialState,[1,numberOfSims]);
for n = 1:simTime
r = rand(1,numberOfSims);
for m = 1:numberOfSims
% for each sim, where does r(m) lie, relative to the boundaries
% in the corresponding row of CT?
[~,~,currentState(m)] = histcounts(r(m),CT(currentState(m),:));
end
end
finalState = currentState;
This took a minute or so to run, but that was a fair amount of work. The first 10 such parallel simulations ended in these states:
finalState(1:10)
ans =
2 4 5 3 5 5 5 5 4 5
Now, how often did our process end in any given state? We can count that using accumarray.
finalStateFreq = accumarray(finalState',ones(numberOfSims,1))/numberOfSims
finalStateFreq =
0.0029
0.0336
0.1957
0.2891
0.4787
If I did a good job here, this should mesh well with the steady state predictions from before. Was the simulation time long enough? Surely yes. This is a small Markov process, with only 5 states. And 10K total sims will be entirely adequate to predict if the result matches the steady state predictions.
V1 =
0.0025378 0.035524 0.19457 0.29167 0.4757
That I got a consistent answer from the simulations with the steady state predictions suggests I did a good job in the simulation. As you would expect from a random simulation, it will only approach the theoretical frequency as the number of simulations goes to infinity. With a little extra effort, I suppose could even have estimated the variance of the counts in each bin, based on the variance of a binomial distribution.
Nazer Hdaifeh on 27 Aug 2020
Thank you for the great discussion. Can you add the case for Continuous time Markove chain with a five states. Thanks

mona faraji on 1 Jun 2013
Edited: Walter Roberson on 20 Apr 2019
chack the following for a 2*2 transition matrix and for 1000 states begining at 1:
transition_probabilities = [0.1 0.9;0.8 0.2]; starting_value = 1; chain_length = 1000;
chain = zeros(1,chain_length);
chain(1)=starting_value;
for i=2:chain_length
this_step_distribution = transition_probabilities(chain(i-1),:);
cumulative_distribution = cumsum(this_step_distribution);
r = rand();
chain(i) = find(cumulative_distribution>r,1);
end
% provides chain = 1 2 1 2 1 2 1 2 1 1 2 1 2 1 2....
Walter Roberson on 8 Jun 2020
Gn Gnk: Could you confirm that you copied mona's code from https://www.mathworks.com/matlabcentral/answers/57961-simulating-a-markov-chain#answer_87358 and pasted it into your command line, but you got an error? What was the error message?

Sean de Wolski on 4 Jan 2013
Edited: Sean de Wolski on 4 Jan 2013
If you also have the emissions matrix, you can use hmmgenerate()
More
Pseudo-ish-code (from my understanding, (disclosure: not a Markov Model expert by any means))
Use a for-loop to loop n times for length you want. S
transC = [zeros(size(trans,1),1), cumsum(trans,2)]; %cumulative sum of rows, we will use this to decide on the next step.
n = 10;
states = zeros(1,n); %storage of states
states(1) = 1; %start at state 1 (or whatever)
for ii = 2:n
%Generate a random number and figure out where it falls in the cumulative sum of that state's trasition matrix row
[~,states(ii)] = histc(rand,transC(states(ii-1),:));
end
Jay Hanuman on 22 Oct 2016
this is first order markov chain or not

Paul Fackler on 21 Aug 2013
You can simulate a Markov chain using the function ddpsimul in my CompEcon toolbox available at www4.ncsu.edu/~pfackler/compecon

Ragini Gupta on 10 Nov 2017
Hey there, I am using Markov Chain model to generate synthetic data. However, I get this error when I simulate it to find the next markov state.
*Edge vector must be monotonically non-decreasing.
Error in MarkovChain2 (line 53) [~,states(ii)] = histc(rand,transC(states(ii-1),:));*
Any ideas where I could be going wrong:?
filename = 'newTestingExcel.xlsx';
edges = linspace(min(Furnace),max(Furnace),8); [counts,bins] = histc(Furnace, edges); [counts,bins] = histc(Furnace, edges);
%to calculate the pdf of each state' f ' is the pdf of each state: [f,xi,bw] = ksdensity(Furnace,edges);
alpha = 2.43; beta = 1; pd = @(f)gampdf(f,alpha,beta); % Target distribution
%# transition matrix trans = full(sparse(bins(1:end-1), bins(2:end), 1)); trans = bsxfun(@rdivide, trans, sum(trans,2));
transCum = [zeros(size(trans,1),1), cumsum(trans,8)]; n = 8; states = zeros(1,n); states(1) = 1; for length = 2:n
[~,states(length)] = histc(rand,transCum(states(ii-1),:));
end

Christine Tobler on 22 Apr 2019
There is a specific class that represents a discrete-time Markov chain in the Econometrics toolbox: dtmc.