```
%matplotlib inline
```

```
import numpy as np
import matplotlib.pyplot as plt
import random as rm
from scipy.stats import norm
import itertools
import matplotlib.patches as mpatches
```

# Stochastic Processes and Applications¶

## Author: Stefan Panev¶

### Abstract¶

This notebook is a basic introduction into Stochastic Processes. It is meant for the general reader that is not very math savvy, like the course participants in the Math Concepts for Developers in SoftUni.

There is a basic definition. Some examples of the most popular types of processes like Random Walk, Brownian Motion or Weiner Process, Poisson Process and Markov chains have been given. Their basic characteristics and examples for some possible applications are stated. For all the examples there are simulations in Python, some are visualized.

The following packages have been used:

- nympy
- matplotlib.pyplot
- random
- scipy.stats
- itertools
- matplotlib.patches

### Introduction¶

A Stochastic Process is a collection of random variables over a period of time. It can be either discrete or continuous type. Techniques from calculus and probability theory are used to study the processes. There are numerous kids of stochastic processes. Several kinds of Stochastic processes are going to be given as an example in this notebook.

Lets define the stochastic process as a collection of random variables in defined on a common probability space $$(\Omega, \mathcal{F}, P)$$
,where omega is the the sample space or all the possible outcomes, mathematical F is the sigma-algebra where each set contains zero or several possible outcomes and P is the probability for the realization of the outcomes.

We define a function with an index t (representing the time) that maps a variable from the set T to a random variable from the state space S or: $$X_{t} : T \to S$$

### Random Walk¶

The first basic process we will look at is the random walk. It is defined as the path created by the steps which have no rule determining their development. The term 'random walk' was coined by the mathematician Karl Pearson (1857 - 1936) in 1905. The process can either be represented as one or several dimensions.

There are numerous applications in finance for modeling stocks and price movements (for further reading https://www.investopedia.com/terms/r/randomwalktheory.asp), the net worth of a gambler, movement of people in the market (such as some agent based modeling), movements of molecules and particles or changes in gene in the genome.

There are different kinds of random walks, if the steps follow the normal distribution it is said that it is a gaussian type. Other variants of a random walk are self-interacting walks, correlated walks, maximal entropy random walk and others.

```
# Initialization
np.random.seed(99)
# Initialize all_walks
all_walks = []
# Simulate random walk 10 times
for i in range(10) :
random_walk = [0]
for x in range(100) :
step = random_walk[-1]
dice = np.random.randint(1,7)
if dice <= 2:
step = step - 1
elif dice <= 5:
step = step + 1
else:
step = step + np.random.randint(1,5)
random_walk.append(step)
# Append random_walk to all_walks
all_walks.append(random_walk)
#print(all_walks)
# Plot random_walk
plt.plot(random_walk)
# Show the plot
plt.show()
```

We generate a simulate 100 random walks. The np.random.seed of 99 is being used that let the reader reproduce the same results. Because we are using the gaussian random generator of integers the walk will be gaussian in nature.

To make thing more interesting we create a basic algorithm, where we though a dice and depending on the result we move with a particular number of steps. We get the number of step from the previous iteration of the for cycle. The steps are decreased with 1 if the dice is less or equal then 2, if there are between 3 and 5 we increase the number of steps with 1, and if we get 6 we increase the number of steps with a random number between 1 and 4. Finally the values of the random walk are printed (optionally, it is commented out) and plotted.

We seem mostly a increasing trend.

```
np.random.seed(99)
all_walks = []
for i in range(20) :
random_walk = [0]
for x in range(100) :
step = random_walk[-1]
dice = np.random.randint(1,7)
if dice <= 2:
step = step - 1
elif dice <= 5:
step = step + 1
else:
step = step + np.random.randint(1,5)
if np.random.rand() <= 0.001 :
step = 0
random_walk.append(step)
all_walks.append(random_walk)
# Convert all_walks to Numpy array: np_aw
np_aw = np.array(all_walks)
# Plot np_aw and show
plt.plot(np_aw)
plt.show()
# Clear the figure
plt.clf()
# Transpose np_aw: np_aw_t
np_aw_t = np.transpose(np_aw)
# Plot np_aw_t and show
plt.plot(np_aw_t)
plt.show()
```

We start with 20 walks. The graphs are much clearer then with 500 walks below. Most paths follow and increasing trend.

```
np.random.seed(99)
all_walks = []
for i in range(500) :
random_walk = [0]
for x in range(100) :
step = random_walk[-1]
dice = np.random.randint(1,7)
if dice <= 2:
step = step - 1
elif dice <= 5:
step = step + 1
else:
step = step + np.random.randint(1,5)
if np.random.rand() <= 0.001 :
step = 0
random_walk.append(step)
all_walks.append(random_walk)
# Convert all_walks to Numpy array: np_aw
np_aw = np.array(all_walks)
# Plot np_aw and show
plt.plot(np_aw)
plt.show()
# Clear the figure
plt.clf()
# Transpose np_aw: np_aw_t
np_aw_t = np.transpose(np_aw)
# Plot np_aw_t and show
plt.plot(np_aw_t)
plt.show()
```

Now we will scale the simulation increasing the number of walks to 500. The walks have been plotted. With so many simulations so the graph is very dense and not very informative, but the Law of Large Numbers is set to apply now.

The trends are again mostly increasing.

```
# Select last row from np_aw_t: ends
ends = np_aw_t[-1]
# Plot histogram of ends, display plot
plt.hist(ends)
plt.show()
#lets caclulate the odds of the ending point being above 40
bool = ends >= 40
greater = sum(bool)
print("The odds of the ending point being above 40 is " + "{:0.2%}".format(greater / 500))
```

We define and ending point of the process, when all the walks are over with the variable ends. It has been plotted and is showing close to the normal distribution, although it seem some asymmetric to the left. At the end we calculate the odds of having and ending point more then 90 which is 20,40%. (this can be done for other percentages).

### Brownian Motion or Wiener Process¶

Both names are used to describe the process. It was first discovered by biologist Robert Brown (1773 - 1858) studying the particle movements under a microscope. However it was mathematically describer by Norbert Wiener (1894 - 1964). The process is a continuous-time stochastic process and a subset of the Levy Process. Martingales are derived from the process.

Lets represent the process with $W_{t}$, with the t index as the time.

- It starts with 0 or $W_{0}$ = 0.
- The increments are independent from the past values $W_{t + u} - W_{t}$ for u > s > 0 are independent form $W_{s}$.
- The increments have normal distribution $W_{t + u} - W_{t}$ $\sim N(O,u)$
- W has a.s. continuous path.

Some of the basic properties are:

The mathematical expectation is 0 or $E(W_{t}) = 0$.

The variance is 0 or $Var(W_{t}) = t$

The covariance is $cov(W_{s}, W_{t}) = min(s, t)$

The correlation is $corr(W_{s}, W_{t}) = \frac{min(s, t)}{\sqrt{st}}$

It can be further altered if we add an additional drift to show changes in a certain direction.

It is used to model white noise, also has major influence for fluid mechanics, quantum mechanics and the Black-Scholes option pricing model (which is not correct in every case).

```
# Brownian Motion steps
for i in range(0, 50):
xt = 0 + norm.rvs(scale=1**2* 4)
#print(xt)
```

Now we generate 50 values for $x_{t}$. It is comment out. If the reader wants to see the print statement can uncomment it.

```
def generate_brownian(x0, n, dt, delta, output=None):
x0 = np.asarray(x0)
r = norm.rvs(size=x0.shape + (n,), scale=delta* np.sqrt(dt))
if output is None:
output = np.empty(r.shape)
np.cumsum(r, axis=-1, out = output)
output += np.expand_dims(x0, axis=-1)
return output
```

The Generate*Brownian Function is defined, with inputs $x*{0}$, $n$, $d*{t}$ and $delta$. The $x*{0}$ = 0 and the formula is: $$X(t+dt)=X(t) + N(0,(delta)^2 dt;t,t+dt) $$

```
x = np.empty((50,101))
x[:, 0] = 100
generate_brownian(x[:, 0], 100, 1, 4, output=x[:,1:])
time = np.linspace(0, 100, 101)
for k in range(50):
plt.plot(time, x[k])
plt.title('Brownian Motion Simulation Multiple')
plt.xlabel('time')
plt.ylabel('x_values')
plt.grid(True)
plt.show()
```

Finally we plot the 50 simulations, we start with $x_{0}$ = 100 and go on for 100 periods. Most simulations stay close to the baseline of 100, but some stray further away.

### Poisson Process¶

The Poisson Process is named after the French mathematician Simon Poisson (1781 - 1840). We define a Poisson Process as random process when it satisfies the following conditions:

- The first value is 0 or $X(0) = 0$.
- The increments are independent. This is also named lack of memory or Markov Property, meaning the the Markov Process is more general.
- The probabilities of x = k follows the poisson distribution (a type of discrete probability distribution) with a probability mass function for k = 0, 1, 2 ...: $$P(X = k) = \frac{\lambda ^ k e^{-\lambda}}{k!}$$ The distribution can be visualized in the following picture. Lambda is also known as an intensity. The higher the lambda the more symmetric the distribution.

Some examples of applications are queueing theory for arrival of customers at a shop, trees in forests, subatomic particles smashing in a particle accelerator or getting calls from a telephone network.

Now we define the poisson*process function with 2 arguments lambdas and number. The lambdas show different intensities, while the number is the time periods that are generated. We will check different outcomes. The output is generated following this function: $$X(t) = X*{0} + \sum*{k = 1}^{t} X*{k}$$

```
def poisson_process(lambdas, number):
X_T = np.random.poisson(lambdas, size=number)
S = [np.sum(X_T[0:i]) for i in range(number)]
X = np.linspace(0, number, number)
graphs = [plt.step(X, S, label="Lambda = %d"%lambdas)[0] for i in range(lambdas)]
graph = plt.step(X, S, label="Lambda = %d"%lambdas)
plt.legend(handles=graph, loc=2)
plt.title('Poisson Process')
plt.xlabel('time')
plt.ylabel('intensity')
#The legend does not display perfectly so the parameters have been split.
```

```
poisson_process(1, 20)
plt.show()
```

```
poisson_process(2, 20)
plt.show()
```

```
poisson_process(5, 20)
plt.show()
```

```
poisson_process(10, 20)
plt.show()
```