mc_fbm.py

import numpy as np
import pandas as pd
import nolds
from fbm import FBM
import matplotlib.pyplot as plt

# number of time steps (days) to predict ahead
n = 30
# number of FBMs to realize
c = 1000

# read in the csv
df = pd.read_csv("eth-cad-max.csv")
print(df.head())
# get the prices col slice from df
prices = np.array(df['price'])
# scale prices
SCALER = 1000
prices = prices/SCALER
start = 0
subset = prices[start:]

# compute the (per timestep) volatility and drift of the time series by log returns
log_ret = []
foo = prices[0]
for i in range(len(prices)-1):
    ret = prices[i+1]/foo
    log_ret.append(np.log(ret))
    foo = prices[i+1]
log_ret = np.array(log_ret)
sigma = np.std(log_ret)
mu = np.average(log_ret)

# compute hurst exponent
H = nolds.hurst_rs(subset)
print("Calculated hurst exponent is " + str(H))

# realize the FBMs
FBMs = []
for i in range(c):
    f = FBM(n=n, hurst=H, length=1, method='hosking')
    f_sample = f.fbm()
    FBMs.append(f_sample)
FBMs = np.array(FBMs)
W_0 = prices[len(prices)-1]
WD = []
means = []
sds = []
# compute the mean and sd's at each time step across all realizations
for i in range(c):
    foo = []
    for t in range(n+1):
        if(t==0):
            continue
        # apply the Wiener process with drift solution to the FBM realizations
        w = W_0*np.exp(mu * t + sigma* FBMs[i,t])
        foo.append(w)
    foo = np.array(foo)
    WD.append(foo)
WD = np.array(WD)
for t in range(n):
    m = WD[:,t]
    foo = np.average(m)
    means.append(foo)
    foo = np.std(m)
    sds.append(foo)

means = np.array(means)
sds = np.array(sds)
# plot predicted means and 95% confidence region (2 sd's)
x = np.linspace(1,n,n)
plt.figure(figsize=(20,10))
# randomly plot every 11th realization
mod = 11
for i in range(c):
    if (i % 11 == 0):
        plt.plot(x,WD[i],linewidth=0.5, color = 'black')
plt.plot(x,means,  linewidth=3, color='r', label='Predicted mean')
plt.fill_between(x, means-2*(sds), means+2*(sds), color = 'grey', label = "95% confidence region")
plt.xlabel('Time step (days)')
plt.ylabel('CAD$ (x1000)')
plt.legend(loc = 2)
plt.savefig("test.png",bbox_inches='tight')
plt.show()