The problem of underlying laws of financial time series remains open still. Markovian stochastic processes cannot model these; otherwise the Hurst exponents would not have such a wide distribution away from .

The figure above shows the Hurst exponent distribution of 1900 stocks daily returns over many years. Now the Hurst exponent can be interpreted in terms of Holder continuity of the series. So a natural idea is to attempt to model the underlying processes in terms of a fractional diffusion with a Gaussian initial, i.e. with initially with $u(0,x)=npdf(0,1)$.

The Green’s function for the fractional diffusion-wave equation was obtained by F. Mainardi in 1996 using Laplace transforms (green-function-fracdiff-mainardi). Recall that a fractional derivative of a causal function is defined for as

and the Laplace transform formula for the fractional derivative proved by M. Caputo in 1969 is

This Laplace transform formula allows us to solve the simple fractional diffusion easily taking Laplace transforms noting :

or

(*)

where and

So now we have the ordinary differential equation if the initial function is a Gaussian or Meixner distribution or any other univariate distribution that might model individual returns. Now recall from elementary differential equation theory that the homogeneous component of (*) is solved easily via the roots of the characteristic polynomial which are 0 and . We can use the following approach to solve this, which is also standard.

:Let and let and and , and rewrite (*) as the matrix first order system . The homogeneous equation produces a matrix , and then we get the particular solution via .

The general solution has matrix and using this, we find that the particular solution with the inhomogeneity is:

So the general solution is

import numpy as np

from numpy import exp

import scipy.integrate as SI

import scipy.special as SS

from talbot import *

from meixnermle import *

def meixnerpdf( x, a=0.1, b=0.1, m=1.1, d=1.0):

G1 = np.abs(SS.gamma(d + (m+1j*x))/a)**2

G2 = SS.gamma(2*d)

A = (2*np.cos(b/2))**(2*d)/(2*a*np.pi*G2)

B = np.exp(b*(x-m)/a)

return A*B*G1

def particular_solution_frac_diffusion( f, s, x, gamma ):

def exp_times_f( s, x, gamma):

return exp( s**(gamma+1) * x ) * f(x)

return exp(-s**(gamma+1)*x)*SI.quad(exp_times_x,0,x)

class function_object:

def __init__(self, xv, yv):

self.xv=xv

self.yv=yv

self.N = len(xv)

def __call__(self,x):

for k in range(1,self.N):

if np.real(x)<np.real(self.xv[0]):

return self.yv[0]

if np.real(self.xv[k-1])<=np.real(x) and np.real(self.xv[k])>np.real(x):

return self.yv[k-1]

if np.real(x)>=np.real(self.xv[self.N-1]):

return self.yv[self.N-1]

N=500

M=50

gamma=0.7

S = np.linspace(0.0,100.0,num=N)

X = np.linspace(-1.0,1.0,num=M)

V = np.zeros((M,N))

W = np.zeros((M,N))

for j in range(M):

for k in range(N):

V[j,k] = particular_solution( meixnerpdf, S[j], x[j], gamma)

for j in range(M):

for k in range(N):

W[j,k] = Talbot(function_object(S,V[j,:]))(-S[j])

# Now W should contain the Laplace inverse of the solution of

# the FDE with a Meixner PDF initial distribution

## Leave a Reply