#!/usr/bin/env python
# coding: utf-8
# # Stock price modelling for the mathematically curious
#
# In this notebook we discuss a mathematical model for stock prices. We present some numerical approximations and illustrate these with simple plots.
# ### Bank balance as an Ordinary Differential Equation
#
# 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
get_ipython().run_line_magic('matplotlib', 'inline')
plt.figure(figsize=(15,5))
plt.xlabel('Time')
plt.ylabel('Bank balance')
sns.lineplot(data=bank_balance)
# ### Fluctuations as a standard Wiener process
#
# 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)
# ### Stock Price as a Stochastic Differential Equation
#
# 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$.
#
#
# ### Itô's formula
#
# 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
# $$
#
#
# ### Solution to our SDE
#
# 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)
# In[ ]: