In this notebook we discuss a mathematical model for stock prices. We present some numerical approximations and illustrate these with simple plots.

Consider the general ordinary differential equation for unknown $U : (0,\infty) \to \mathbb{R}$,

$$ \tag{ODE} \frac{dU(t)}{dt} = f(t,U(t)), \quad t > 0, $$

with a given initial value $U(0) = U_0$ and function $f : (0,\infty) \times \mathbb{R} \to \mathbb{R}$.

By integrating, we see that it is equivalent with the integral equation

$$ \tag{IE} U(t) = U_0 + \int_0^t f(s,U(s)) \, ds , \quad t > 0. $$

For example, we could be modelling our bank account balance $B$ by

$$ \frac{dB(t)}{dt} = \mu B(t) + F(t), \quad t > 0, $$

where $\mu$ is the interest rate and $F$ stands for, say, your salary or expenses.

It is easy to see that the solution to this equation is given by

$$ B(t) = e^{\mu t}B_0 + \int_0^t e^{\mu (t-s)} F(s) \, ds , \quad t > 0, $$

where $B_0$ is the initial balance.

Remark: For a bank account, a discrete model would probably be a more realistic choice, but we wish to emphasize the analogy between this and the continuous time model for stock price.

In [1]:

```
# We inspect our model over a time frame
import pandas as pd
import numpy as np
start_date = '2017-01-01'
end_date = '2017-12-31'
date_range = pd.date_range(start=start_date, end=end_date)
# We assume that the daily earnings and expenses vary somewhat randomly
daily_earnings_avg = 10.0 # average daily earnings
daily_expenses_scale = 50.0 # scale of unexpected expenses
F_df = pd.DataFrame(data=np.random.normal(loc=daily_earnings_avg,
scale=daily_expenses_scale,
size=len(date_range)),
index=date_range)
def F(s):
return F_df.loc[s].values
B_0 = 10000 # Initial bank balance at 10000 euros
mu = 0.01 / 365 # 1 percent annual interest rate transformed into daily interest rate
def integrand(t,s):
return np.exp(mu * (t - s).days) * F(s)
def integral(t):
return sum([integrand(t,s) for s in pd.date_range(start=start_date, end=t)])
def B(t):
return np.exp(mu * (t - pd.to_datetime(start_date)).days) * B_0 + integral(t)
bank_balance = pd.DataFrame(data=[np.round(B(t),2) for t in date_range], index=date_range, columns=['bank_balance'])
```

In [2]:

```
# We plot the bank balance over our chosen time frame
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=(15,5))
plt.xlabel('Time')
plt.ylabel('Bank balance')
sns.lineplot(data=bank_balance)
```

Out[2]:

What if we wanted to invest in some stocks intead? It is generally agreed that stock values fluctuate quite randomly. We could therefore try to model the stock price $S_t$ at time $t$ by

$$ \tag{1} \frac{dS_t}{dt} = \mu S_t + \sigma S_t \xi_t , \quad t > 0, $$

where $\mu$ and $\sigma$ stand for expected rate of return and volatility, respectively, and $\xi_t$ denotes random "noise".

It is not immediately clear how random noise should be desrcibed mathematically. The generally accepted way is to view it formally as a time derivative of a *standard Wiener process* $W$, that is

$$ \xi_t = \frac{dW_t}{dt} . $$

By definition, the standard Wiener process

- has continuous paths,
- $W_0 = 0$ almost surely,
- has independent increments so that $W_t - W_s$ is independent from $W_u - W_v$ whenever $t > s > u > v$,
- $W_t - W_s$ is normally distributed with mean $0$ and variance $t-s$, whenever $t > s$

The last property has the effect that the displacement of Wiener process scales, not linearly with time, but with its square root, that is

$$ dW_t \sim (dt)^{1/2}. $$

In [3]:

```
# Here we plot a few sample paths of a discretized Wiener process over (0,1)
N = 100
K = 5
dt = 1 / N
dW = np.random.normal(loc=0.0, scale=np.sqrt(dt), size=(N,K))
W = np.cumsum(dW, axis=0)
W_df = pd.DataFrame(data=W, index=np.array(range(100)) / N)
plt.figure(figsize=(15,5))
sns.lineplot(data=W_df,
palette=sns.light_palette("blue", n_colors=K),
legend=False)
```

Out[3]:

Proceeding formally, we may then multiply both sides of equation (1) by $dt$ and arrive at

$$ \tag{SDE} dS_t = \mu S_t \, dt + \sigma S_t \, dW_t, \quad t > 0, $$

which can be interpreted rigorously as the integral equation

$$ S_t = S_0 + \mu \int_0^t S_s \, ds + \sigma \int_0^t S_s \, dW_s , \quad t > 0. $$

It remains to be defined what we mean by the stochastic integral on the right-hand side. We will interpret it as an *ItÃ´ integral*, that is, as the limit of approximating sums

$$ \int_0^t S_s \, dW_s \sim \sum_{k=0}^N S_{s_k}(W_{s_{k+1}} - W_{s_k}) $$

with increasingly fine-grained partitions $0 = s_0 < \cdots < s_k < s_{k+1} < \cdots < s_N = t$.

Recall the Taylor expansion of a smooth function $f$ around $(t',x')$:

$$ f(t,x)-f(t',x')=\frac{\partial f}{\partial t}(t',x')(t-t')+\frac{\partial f}{\partial x}(t',x')(x-x')+\frac{1}{2}\frac{\partial^2 f}{\partial t^2}(t',x')(t-t')^2+\frac{1}{2}\frac{\partial^2 f}{\partial x^2}(t',x')(x-x')^2+\cdots $$

As $(t,x)$ approaches $(t',x')$ we arrive at the first order approximation

$$ df(t,x) = \frac{\partial f}{\partial t} (t,x) \, dt + \frac{\partial f}{\partial x}(t,x) \, dx $$

A stochastic version of this arises from substituting $W_t$ for $x$. Now, since $dW_t \sim (dt)^{1/2}$, also the second order terms contribute to the first order approximation:,

$$ \frac{1}{2}\frac{\partial^2 f}{\partial t^2}(t,W_t)\,dt^2+\frac{1}{2}\frac{\partial^2 f}{\partial x^2}(t,W_t)\,dW_t^2\sim\frac{1}{2}\frac{\partial^2 f}{\partial x^2}(t,W_t)\,dt . $$

This leads us to ItÃ´'s formula:

$$ df(t,W_t) = \frac{\partial f}{\partial t}(t,W_t) \, dt + \frac{\partial f}{\partial x} (t,W_t) \, dW_t + \frac{1}{2} \frac{\partial^2 f}{\partial x^2} (t,W_t) \, dt $$

We are now ready to face the SDE modelling our stock prices. Indeed, let us use ItÃ´'s formula to compute the differential of

$$ S_t = e^{W_t} . $$

Without difficulty, we see that

$$ de^{W_t} = e^{W_t} \, dW_t + \frac{1}{2} e^{W_t} \, dt $$

so that we've arrived to the solution of

$$ dS_t = \frac{1}{2} S_t \, dt + S_t \, dW_t . $$

This is a special case of our equation with

$$ \begin{cases} \mu = 1/2, \\ \sigma = 1, \\ S_0 = 1. \end{cases} $$

With some thought, we notice that the general case is solved by

$$ S_t = S_0 \exp((\mu - \sigma^2 / 2)t + \sigma W_t) . $$

This process is known as *Geometric Brownian motion*.

In [4]:

```
# Here we plot a few sample paths of the solution S_t
N = len(date_range)
K = 5
dt = 1 / N
dW = np.random.normal(loc=0.0, scale=np.sqrt(dt), size=(N,K))
W_df = pd.DataFrame(data=np.cumsum(dW, axis=0), index=date_range)
def W(t):
return W_df.loc[t]
# Try out different values for the parameters
S_0 = 100 # Initial stock value at 100 euros
sigma = 0.1
mu = 0.0055
def S(t, mu, sigma, W):
return S_0 * np.exp((mu - sigma**2 / 2) * (t - pd.to_datetime(start_date)).days + sigma * W(t))
stock_price = pd.DataFrame(data=[np.round(S(t, mu, sigma, W),2) for t in date_range],
index=date_range)
plt.figure(figsize=(15,5))
plt.xlabel('Time')
plt.ylabel('Stock value')
sns.lineplot(data=stock_price,
palette=sns.light_palette("blue", n_colors=K),
legend=False)
```

Out[4]: