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
sums to 1.0 i.e.*P*

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:*

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