```
import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
```


```
#matplotlib parameters
plt.rcParams["figure.figsize"] = (12, 8)
plt.rcParams.update({'font.size': 14})
```

# **What is a stochastic process?**
**A stochastic process is a collection of random
variables: $\{X_t,t \in I \},$ Where: $X_t$ is the set of random variables at time $t$, and $I$ is the index set of the process.**

**Discrete time stochastic processes which we will focus on in this tutorial are sequences of random variables.**

## **Simulation of Gambler's Ruin**

**We can use the famous example of gambler's ruin which is a stochastic process.**

**In our case, a gambler starts off with \$50. To keep things simple, they can only bet in increments of \$1. They can only win or lose \$1 per bet. They will keep gambling until they either lost all their money (leave with \$0) or if they win \$100. We can formalize this with the following notation:**

$\begin{equation}
  X_{g} =
    \begin{cases}
      \$+1, & \text{with a probability of 50%}\\
      \$-1, & \text{with a probability of 50%}\\
    \end{cases}       
\end{equation}$

**Where: $X_g$ is the gambling outcome in dollars.**

<sup>Source: [Introduction to Stochastic Processes with R](https://www.amazon.com/Introduction-Stochastic-Processes-Robert-Dobrow/dp/1118740653) by Robert P. Dobrow</sup>


```
def gamblers_ruin():
  gambling_money = 50
  gambling_goal = 100
  gambling_simulation = []

  while gambling_money in range(1,gambling_goal):
    bet_size = 1
    w_or_l = random.randrange(-1, 2, step = 2)
    gambling_money += bet_size * w_or_l
    gambling_simulation.append(gambling_money)
  return gambling_simulation
```


```
plt.plot(gamblers_ruin())
plt.yticks(np.arange(-20,120,10))
plt.axhline(y=0, color='r', linestyle='-')
plt.axhline(y=100, color='black', linestyle='-')
plt.xlabel('Number of bets')
plt.ylabel('Winnings')
plt.title('Gambling Simulation');
```


```
def prob_of_ruin(gambling_goal, initial_gambling_money):
  return (gambling_goal - initial_gambling_money)/gambling_goal
```


```
prob_of_ruin(100,50)
```




    0.5




```
sim_list = []

while len(sim_list) < 500:
  sim_list.append(gamblers_ruin()[-1])

np.mean(sim_list)
```




    54.2



<sup>Source: [Introduction to Stochastic Processes with R](https://www.amazon.com/Introduction-Stochastic-Processes-Robert-Dobrow/dp/1118740653) by Robert P. Dobrow</sup>

# **Markov Chains**

**A Markov chain is a type of stochastic process.**

**A Markov chain is a collection of random variables ($X_t$) where the future states ($j$) only depend on the current state ($i$). Markov chains can be either discrete or continuous.**

**For a Markov chain transition matrix (denoted as $P$):**

 - **Each row must add to one, where: $\sum\limits_{j}P_{ij} = 1$.**
 - **The probabilities must be non-negative where: $P_{ij} \geq 0 \ \ \forall \ \ i,j$**

<sup>Source: [Introduction to Stochastic Processes with R](https://www.amazon.com/Introduction-Stochastic-Processes-Robert-Dobrow/dp/1118740653) by Robert P. Dobrow</sup>

<sup>Source: [Markov Chain](https://mathworld.wolfram.com/MarkovChain.html?utm_source=twitterfeed&utm_medium=twitter) from Wolfram MathWorld</sup>

## **Markov Chain Simulation**

**A Markov chain can be simulated from an initial distribution and transition matrix. In our case, the initial state is New York City. From the initial state we can travel to: Paris, Cairo, Seoul or even within New York City.**

**The transition matrix contains the one step transition probabilities of
moving from state to state.**

<sup>Source: [Introduction to Stochastic Processes with R](https://www.amazon.com/Introduction-Stochastic-Processes-Robert-Dobrow/dp/1118740653) by Robert P. Dobrow</sup>


```
mc_example = {'NYC': [.25,0,.75,1],
              'Paris': [.25,.25,0,0],
              'Cairo': [.25,.25,.25,0],
              'Seoul': [.25,.5,0,0]}

mc = pd.DataFrame(data = mc_example, index = ['NYC', 'Paris', 'Cairo', 'Seoul'])
```

### **Markov Transition Graph**



**We can formalize the movement for the Markov chain at the initial starting point (New York City) with the following notation:**

$P(X_1 = \text{New York}|X_0 = \text{New York})=P(X_1 = \text{Paris}|X_0 = \text{New York})=P(X_1 = \text{Cairo}|X_0 = \text{New York}) = P(X_1 = \text{Seoul}|X_0 = \text{New York}) = 25\%$


```
travel_sim = []
travel_sim.append(mc.iloc[0].index[0])
city = np.random.choice(mc.iloc[0].index, p = mc.iloc[0])
travel_sim.append(city)

while len(travel_sim) < 25:
  city = np.random.choice(mc.iloc[mc.index.get_loc(city)].index, p = mc.iloc[mc.index.get_loc(city)])
  travel_sim.append(city)
```


```
travel_sim
```




    ['NYC',
     'Paris',
     'Paris',
     'Seoul',
     'NYC',
     'Seoul',
     'NYC',
     'Seoul',
     'NYC',
     'Cairo',
     'NYC',
     'Seoul',
     'NYC',
     'Seoul',
     'NYC',
     'Seoul',
     'NYC',
     'Paris',
     'Seoul',
     'NYC',
     'NYC',
     'Cairo',
     'NYC',
     'Seoul',
     'NYC']




```
mc
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>NYC</th>
      <th>Paris</th>
      <th>Cairo</th>
      <th>Seoul</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>NYC</th>
      <td>0.25</td>
      <td>0.25</td>
      <td>0.25</td>
      <td>0.25</td>
    </tr>
    <tr>
      <th>Paris</th>
      <td>0.00</td>
      <td>0.25</td>
      <td>0.25</td>
      <td>0.50</td>
    </tr>
    <tr>
      <th>Cairo</th>
      <td>0.75</td>
      <td>0.00</td>
      <td>0.25</td>
      <td>0.00</td>
    </tr>
    <tr>
      <th>Seoul</th>
      <td>1.00</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.00</td>
    </tr>
  </tbody>
</table>
</div>



### **Memorylessness: The future is independent of the past given the present**

**An important feature of Markov chains is that they are memorylessness. If we look at our example, the only state that matters is the current state. If our salesperson started in New York City ($X_0 = NYC$) and traveled to Paris ($X_1 = Paris$), then the movement to the next city ($X_2$) only depends on the probabilities of travel from Paris. The fact that the salesperson originally started in New York does not effect the movement to $X_2.$** 

## **$n$ Step Transition Matrix**


```
mc.to_numpy()
```




    array([[0.25, 0.25, 0.25, 0.25],
           [0.  , 0.25, 0.25, 0.5 ],
           [0.75, 0.  , 0.25, 0.  ],
           [1.  , 0.  , 0.  , 0.  ]])




```
def matrix_power(matrix, power):
  if power == 0:
    return np.identity(len(matrix))
  elif power == 1:
    return matrix
  else:
    return np.dot(matrix, matrix_power(matrix, power-1))
```

<sup>Source: [Introduction to Stochastic Processes with R](https://www.amazon.com/Introduction-Stochastic-Processes-Robert-Dobrow/dp/1118740653) by Robert P. Dobrow</sup>


```
matrix_power(mc.to_numpy(), 2)
```




    array([[0.5   , 0.125 , 0.1875, 0.1875],
           [0.6875, 0.0625, 0.125 , 0.125 ],
           [0.375 , 0.1875, 0.25  , 0.1875],
           [0.25  , 0.25  , 0.25  , 0.25  ]])




```
np.dot(np.dot(mc.to_numpy(), mc.to_numpy()), mc.to_numpy())
```




    array([[0.453125, 0.15625 , 0.203125, 0.1875  ],
           [0.390625, 0.1875  , 0.21875 , 0.203125],
           [0.46875 , 0.140625, 0.203125, 0.1875  ],
           [0.5     , 0.125   , 0.1875  , 0.1875  ]])




```
np.dot(mc.to_numpy(), mc.to_numpy())
```




    array([[0.5   , 0.125 , 0.1875, 0.1875],
           [0.6875, 0.0625, 0.125 , 0.125 ],
           [0.375 , 0.1875, 0.25  , 0.1875],
           [0.25  , 0.25  , 0.25  , 0.25  ]])




```
matrix_power(mc.to_numpy(), 3).sum(axis=1)
```




    array([1., 1., 1., 1.])



**What the above 2-step transition matrix is telling us is if we look at state $i$ for New York, in two steps there is a 50% probability that the salesperson ends up back in New York City, a 12.5% probability that they go onto Paris and an 18.75% chance that they end up in Cairo or Seoul.**


```
for i in range(1,10,1):
  print(f'n Step Transition Matrix at the nth power {i}\n', matrix_power(mc.to_numpy(), i),'\n')
```

    n Step Transition Matrix at the nth power 1
     [[0.25 0.25 0.25 0.25]
     [0.   0.25 0.25 0.5 ]
     [0.75 0.   0.25 0.  ]
     [1.   0.   0.   0.  ]] 
    
    n Step Transition Matrix at the nth power 2
     [[0.5    0.125  0.1875 0.1875]
     [0.6875 0.0625 0.125  0.125 ]
     [0.375  0.1875 0.25   0.1875]
     [0.25   0.25   0.25   0.25  ]] 
    
    n Step Transition Matrix at the nth power 3
     [[0.453125 0.15625  0.203125 0.1875  ]
     [0.390625 0.1875   0.21875  0.203125]
     [0.46875  0.140625 0.203125 0.1875  ]
     [0.5      0.125    0.1875   0.1875  ]] 
    
    n Step Transition Matrix at the nth power 4
     [[0.453125   0.15234375 0.203125   0.19140625]
     [0.46484375 0.14453125 0.19921875 0.19140625]
     [0.45703125 0.15234375 0.203125   0.1875    ]
     [0.453125   0.15625    0.203125   0.1875    ]] 
    
    n Step Transition Matrix at the nth power 5
     [[0.45703125 0.15136719 0.20214844 0.18945312]
     [0.45703125 0.15234375 0.20214844 0.18847656]
     [0.45410156 0.15234375 0.203125   0.19042969]
     [0.453125   0.15234375 0.203125   0.19140625]] 
    
    n Step Transition Matrix at the nth power 6
     [[0.45532227 0.15209961 0.20263672 0.18994141]
     [0.4543457  0.15234375 0.20288086 0.19042969]
     [0.45629883 0.15161133 0.20239258 0.18969727]
     [0.45703125 0.15136719 0.20214844 0.18945312]] 
    
    n Step Transition Matrix at the nth power 7
     [[0.45574951 0.15185547 0.20251465 0.18988037]
     [0.45617676 0.15167236 0.20239258 0.1897583 ]
     [0.45556641 0.15197754 0.20257568 0.18988037]
     [0.45532227 0.15209961 0.20263672 0.18994141]] 
    
    n Step Transition Matrix at the nth power 8
     [[0.45570374 0.15190125 0.20252991 0.18986511]
     [0.45559692 0.15196228 0.20256042 0.18988037]
     [0.45570374 0.15188599 0.20252991 0.18988037]
     [0.45574951 0.15185547 0.20251465 0.18988037]] 
    
    n Step Transition Matrix at the nth power 9
     [[0.45568848 0.15190125 0.20253372 0.18987656]
     [0.45569992 0.1518898  0.20252991 0.18988037]
     [0.45570374 0.15189743 0.20252991 0.18986893]
     [0.45570374 0.15190125 0.20252991 0.18986511]] 
    


**We can see above that as the number of steps increase, the probabilities converge to what is called a stationary distribution (also known as a steady state vector).**

**Let's say that this time our salesperson is starting their trip from Seoul. In this case, we can denote the initial distribution as: $\alpha = [0,0,0,1]$. Now let's find what the probability of ending up back in Seoul is two trips from now.**


```
initial_dist = np.asarray([0,0,0,1])

mc_p2 = matrix_power(mc.to_numpy(),2)

np.dot(initial_dist,mc_p2)
```




    array([0.25, 0.25, 0.25, 0.25])



# **References and Additional Learning**

## **Online Course**

- **[Introduction to Probability](https://ocw.mit.edu/resources/res-6-012-introduction-to-probability-spring-2018/) from MIT OpenCourseWare**

## **Textbook**
- **[Introduction to Stochastic Processes with R](https://www.amazon.com/Introduction-Stochastic-Processes-Robert-Dobrow/dp/1118740653) by Robert P. Dobrow**

## **Videos**
- **[Markov Chains](https://www.youtube.com/watch?v=uvYTGEZQTEs) by patrickJMT**

- **[Markov Steady-State Vectors](https://www.youtube.com/watch?v=D2wvhwOTWIo) by  Brandon Foltz**

- **[Origin of Markov chains](https://www.youtube.com/watch?v=Ws63I3F7Moc) from Khan Academy Labs**

# **Connect**

- **Feel free to connect with Adrian on [YouTube](https://www.youtube.com/channel/UCPuDxI3xb_ryUUMfkm0jsRA), [LinkedIn](https://www.linkedin.com/in/adrian-dolinay-frm-96a289106/), [Twitter](https://twitter.com/DolinayG) and [GitHub](https://github.com/ad17171717). Happy coding!**

683eeccf9c271e312f33125edf583004e7a4428b
ipynb
Jupyter Notebook
Python YouTube Tutorials/An_Intro_to_Markov_chains_with_Python!.ipynb
tudev/Workshops-2020-2021
e708f50cc6620bb8903943588260c9d125a9bafb
MIT
2020-09-01T00:19:22.000Z
2022-03-31T05:38:28.000Z
Python YouTube Tutorials/An_Intro_to_Markov_chains_with_Python!.ipynb
ad17171717/Workshops-2020-2021
f0bcd92c9171ccc43a1a7ccb5ec645c4e6bffe01
MIT
Python YouTube Tutorials/An_Intro_to_Markov_chains_with_Python!.ipynb
ad17171717/Workshops-2020-2021
f0bcd92c9171ccc43a1a7ccb5ec645c4e6bffe01
MIT
2020-08-28T02:20:39.000Z
2022-03-31T05:38:30.000Z
Qwen/Qwen-72B
1. YES
2. YES