#!/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[ ]: