Feeds:
Posts

FRACTIONAL NONLINEAR DIFFUSION EVOLUTION OF GAUSSIAN

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 $H=1/2$.

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. $((d/dt)^\gamma + A)u=F(u)$ with $F=0$ 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 $n < \gamma as

$(d/dt)^\gamma f(t) = \frac{1}{\Gamma(n-\gamma)}\int_0^t \frac{f^{(n)}(s)}{(t-s)^{(\gamma+1-n)}} ds$

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

$L( (d/dt)^\gamma f)(s) = s^\gamma Lf(s) + \sum_{k=0}^{n-1} s^{\gamma-1-k}f^{(k)}(0^+)$

This Laplace transform formula allows us to solve the simple fractional diffusion $((d/dt)^\gamma + d^2/dx^2) u = 0$ easily taking Laplace transforms noting $L(1)(s)= 1/s$:

$s^{\gamma} Lu(x,s) - s^{\gamma-1} u(x,0^+) + s^{-1} L\frac{d^2}{dx^2}u(x,s) = 0$

or

$g_{xx}(x,s) + s^{\gamma+1} g(x,s) = s^{\gamma} f(x)$ (*)

where $g(x,s) = Lu(x,s)$ and $f(x) = s^{\gamma} u(x,0^+)$

So now we have the ordinary differential equation if the initial function $u(x,0^+)$ 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 $s^{\gamma+1}$.  We can use the following approach to solve this, which is also standard.

:Let $A = \diag( s^{\gamma+1}, 1)$ and let $x_1(s) = g(s)$ and $x_2(s) = g's(s)$ and $\Phi = [ s^\gamma f(x), 0]^T$, and rewrite (*) as the matrix first order system $X' + AX = \Phi$.  The homogeneous equation produces a matrix $Z$, and then we get the particular solution via $Z\int Z^{-1} \Phi$.

The general solution has matrix $Z = \diag( \exp(-s^{\gamma+1}x), 1)$ and using this, we find that the particular solution with the inhomogeneity is:

$g_p(x,s) = \exp( -s^{\gamma+1}x) \int_0^x \exp(s^{\gamma+1} y) s^{\gamma} u(y,0^+) dy$

So the general solution is

$g(x,s) = c_1 + c_2 e^{-s^{\gamma+1} x} + g_p(x,s)$

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)

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