And a tutorial on how to simulate a discrete time Markov process using Python
A Discrete Time Markov Chain can be used to describe the behavior of a system that jumps from one state to another state with a certain probability, and this probability of transition to the next state depends only on what state the system is in currently, i.e. it does not depend on which states the system was in prior to the current state.
The above definition pretty much sums up the long and short of what is a Discrete Time Markov Process or a Markov Chain as it’s often called.
But the above description is deceptively simple. As we’ll see, even the simplest of Markov Chains can create richly interesting patterns of system behavior.
We’ll use the terms Markov process and Markov chain interchangeably.
A simple 2-state Markov chain

The Markov chain shown above has two states, or regimes as they are sometimes called: +1 and -1. There are four types of state transitions possible between the two states:
- State +1 to state +1: This transition happens with probability p_11
- State +1 to State -1 with transition probability p_12
- State -1 to State +1 with transition probability p_21
- State -1 to State -1 with transition probability p_22
The above figure is known as the state transition diagram of the Markov process.
The following diagram shows another way to represent this Markov process. The X-axis is the time axis. The bubbles represent the different possible states that the process could be in at each time step.

In reality, at each time step, the process will be in exactly one of the many possible states. In the following diagram, the black bubbles depict one such realization of the Markov process across 4 successive time steps:

There are many such realizations possible. In a 2-state Markov process, there are 2^N possible realizations of the Markov chain over N time steps.
By illustrating the march of a Markov process along the time axis, we glean the following important property of a Markov process:
A realization of a Markov chain along the time dimension is a time series.
The state transition matrix
In a 2-state Markov chain, there are four possible state transitions and the corresponding transition probabilities. We can represent them in a state transition matrix P as follows:

For a Markov chain over n-states (1,2,3,…,n), the transition matrix looks like this:

The state transition matrix has the following two important properties:
- Since each element p_ij is a probability, 0 ≤ p_ij ≤ 1
- Each row of P sums to 1.0 i.e.

This is because the row index represents the source state at time t and the column index represents the destination state at time (t+1). If the process is in source state i at time t, at (t+1), it has to be in one of the allowed set of states (1,2,3,…,n).
Thus, we can restate the transition matrix of a 2-state Markov process as follows:

Introducing the Markov distributed random variable
We will now introduce a random variable X_t. The suffix t in X_t denotes the time step. At each time step t, X_t takes a value from the state space [1,2,3,…,n] as per some probability distribution. One possible sequence of values that X_t takes is {X_0=1, X_1=3, X_2=4, X_3=0,…,X_t=k}.
Transition probabilities as conditional probabilities
At time t, the probability that X_t takes some value j given that X_t had taken a value i at the previous time step (t-1) is given by the following conditional probability:

Indeed, this is the transition probability p_ij.
The Markov property
The Markov property states that p_ij is independent of the state in which the system was at times (t-2), (t-3),…,0. The Markov property is stated as follows:

n-step transition probabilities
The state transition matrix P has this nice property by which if you multiply it with itself k times, then the matrix P^k represents the probability that system will be in state j after hopping through k number of transitions starting from state i.
For example, consider the following transition probabilities matrix for a 2-step Markov process:

The 2-step transition probabilities are calculated as follows:

In P², p_11=0.625 is the probability of returning to state 1 after having traversed through two states starting from state 1. p_12=0.375 is the probability of reaching state 2 in exactly two time steps starting from state 1. And so one.
The state probability distribution
We could continue multiplying P with itself forever to see how the n-step probabilities change over time. What would be more interesting is if we could know what is the unconditional probability distribution of X_t at each time step t.
For example, in our 2-step Markov process example, at each time step t, what is the probability with which X_t could be +1 or -1? These probabilities constitute what is known as the state probability distribution of the Markov variable X_t at time t and it is denoted by π_t, (or δ_t). For example, for the 2-state Markov process consisting of states +1 and -1, the state distribution is given as follows:

An n-state Markov process that operates over the states [1,2,3,…,n] could be in any one of those n states at time t, so π_t is a vector of length n as follows:

In general, each element of π_t can be referenced using the notation π_jt. π_jt is indexed using two variables j and t, indicating the unconditional probability π of the process being in state j at time t:
![The vector π contains the unconditional probabilities of the Markov process being in state j=[1,2,…n] at time t](https://timeseriesreasoning.files.wordpress.com/2021/11/2b442-175w1snzan5a_pjh2rl48ba.png)
Since π_t is a probability distribution, its elements always sum up to 1.0:

How to calculate π_t
Typically, one assumes a certain value for π_0 which is the probability vector for state 0. And given π_0, it can be shown that π_t can be calculated as follows:

The idea behind computing π_t is to matrix-multiply the state transition matrix with itself t number of times to get the t-step transition probabilities, and multiply the t-step transition probabilities with the unconditional probability distribution π_0 at t=0.
Notice the following two things about the above equation:
The state probability distribution of a Markov process at time t depends on:
- The initial probability distribution π_0 at time t=0, and
- It depends on time step t. Thus, the probability distribution evolves over time.
Markov Process Example
Let’s work through all the concepts learnt so far using an example.
Suppose the stock of Acme Corporation increases or decreases as per the following rules:
- As compared to the previous day’s closing price, the current day’s closing price is either higher or lower by some percentage points. For, e.g. on some day, Acme’s stock closes 2.3% higher than the previous day. On other days, it might close 0.8% lower than the previous day’s closing price.
- Given any sequence of 3 days denoted by ‘day before’, ‘yesterday’ and ‘today’, if Acme has closed higher yesterday than the day before, it will close higher today than yesterday with a probability p_11. Consequently, it would close lower today than yesterday with a probability (1 — p_11).
- Let’s also assume that if Acme has closed lower yesterday than the day before, the probability that it would once again close lower today is p_22. And therefore, the probability that it would close higher today than it did yesterday is (1 — p22).
You may have guessed that Acme’s stock price can be modeled using a 2-state Markov process. It’s state transition diagram looks like this:

And it’s transition probabilities matrix P looks like this:

Now let’s pin down some sample values for these probabilities. Suppose the probability of Acme’s price’s moving up two days in a row is 0.6 and the chances of the price moving down for two days in a row is 0.25. Thus,
p_11=0.6 and p_22=0.25
And therefore, we have:
p_12=1 — p_11 = 0.4, and, p_21=1 — p22 = 0.75
The state transition matrix for Acme’s stock is as follows:

The state transition diagram of the Markov process model for Acme’s stock is as follows:

We would now like to calculate the following unconditional probability distribution:

Let’s assume an initial value of π_0=[0.5,0.5], i.e. an equal probability of it’s closing higher (+1 state) or lower ( — 1 state) on listing day, as compared to the IPO price. If we repeatedly apply the formula for π_t=π_0*P^t, we will get the probability distribution vector π_t for each time step t.
The following Python code calculates π_t:
import numpy as np
from matplotlib import pyplot as plt
#initialize the transition matrix P
P=np.array([[0.6,0.4],[0.75,0.25]])
#initialize pi_0
pi_0=np.array([0.5, 0.5])
#set up the array to accumulate the state probabilities at times t=1 to 10
pi=[]
pi.append(pi_0)
P_mul=P.copy()
#calculate the state probability for each t and store it away
for i in range(10):
P_mul=np.matmul(P_mul,P)
pi_t = np.matmul(pi_0,P_mul)
pi.append(pi_t)
pi = np.array(pi)
The above code stores all the calculated π_t vectors in the array π. Let’s separately plot the components π_1t and π_2t of π_t= [ π_1t, π_2t] for t=1 through 10.
#plot pi_2t = P(X_t = -1) versus t
fig = plt.figure()
fig.suptitle('Probability of closing lower than previous day\'s close')
plt.plot(range(len(pi)), pi[:,1])
plt.show()
#plot pi_1t = P(X_t = +1) versus t
fig = plt.figure()
fig.suptitle('Probability of closing higher than previous day\'s close')
plt.plot(range(len(pi)), pi[:,0])
plt.show()
We get the following two plots:

We see that in just a few time steps, the unconditional probabilities vector π_t= [ π_1t, π_2t] settles down into the steady state value [0.65217391, 0.34782609].
In fact, it can be shown that if a two state Markov process with the following transition matrix is run for ‘long time’:

Then the state probability distribution π_t steady-states to the following constant (limiting) probability distribution that is independent of both t and the initial probability distribution π_0:

Simulating Acme’s stock price movement using a Markov process
Let’s simulate Acme Corp’s stock price movement at each time step t using the probability distribution π_t at time t. The simulation procedure goes like this:
- Assume Acme’s IPO price to be $100.
- Set the initial probability distribution:
π_0=[P(X_t = +1)=0.5, P(X_t = -1) = 0.5]. - Map the +1 Markov state to the action of increasing Acme’s previous day’s closing price by a random percentage in the interval [0.0%, 2.0%].
- Map the -1 Markov state to reducing Acme’s previous day’s closing price by a similar random percentage in the interval [0.0%, 2.0%].
- At each time step t, calculate the state probability distribution π_t=[P(X_t = +1), P(X_t = -1) = 0.5], by using the formula π_t = π_0*P^t.
- Generate a uniformly distributed random number in the interval [0, 1.0]. If this number is less than or equal to π_1t = P(X_t = +1), increase the closing price from the previous time step by a random percentage in the interval [0.0%, 2.0%], else decrease the previous closing price by the same random percentage. Repeat this procedure for as many time steps as needed.
Here is the plot of Acme’s stock price changes over 365 trading days:

Here is the source code for generating the above plot:
#Simulate the closing price of a company
closing_price = 100.0
#initialize pi_0
pi_0=np.array([0.5, 0.5])
#create a random delta in the range [0, 2.0]
delta = random.random() * 2
#generate a random number in the range [0.0, 1.0]
r = random.random()
#if r <= P(X_t = +1), increase the closing price by delta,
#else decrease the closing price by delta
if r <= pi_0[0]:
closing_price = closing_price*(100+delta)/100
else:
closing_price = math.max(closing_price*(100-delta)/100,1.0)
#accumulate the new closing price
closing_prices = [closing_price]
P_mul=P.copy()
T=365
#now repeat this procedure 365 times
for i in range(T):
#calculate the i-step transition matrix P^i
P_mul=np.matmul(P_mul,P)
#multiply it by pi_0 to get the state probability for time i
pi_t = np.matmul(pi_0,P_mul)
# create a random delta in the range [0, 2.0]
delta = random.random() * 2
# generate a random number in the range [0.0, 1.0]
r = random.random()
# if r <= P(X_t = +1), increase the closing price by delta,
# else decrease the closing price by delta
if r <= pi_t[0]:
closing_price = math.max(closing_price*(100+delta)/100,1.0)
else:
closing_price = closing_price*(100-delta)/100
# accumulate the new closing price
closing_prices.append(closing_price)
#plot all the accumulated closing prices
fig = plt.figure()
fig.suptitle('Closing price of Acme corporation simulated using a Markov process')
plt.xlabel('time t')
plt.ylabel('Closing price')
plt.plot(range(T+1), closing_prices)
plt.show()
Citations and Copyrights
Images
All images are copyright Sachin Date under CC-BY-NC-SA, unless a different source and copyright are mentioned underneath the image.
PREVIOUS: A Deeper Dive Into The Poisson Probability Distribution Function
NEXT: Hidden Markov Models