
12:48 PM (7 hours ago)



Dear Professor Donoho,
Continuoustime Fourier inversion is the critical step of the applications I have been working on. Fast and accurate computations for these need more careful attention. I follow the attached papers to produce reliable code:
Thus I need to code up usable approximations of inverse continous Fourier transforms. The fractional diffusion Green’s functions and most Levy distributions have computable formulae only for Fourier transforms, so this is a necessary deviation from examination of Hurst exponents.

2:47 PM (5 hours ago)



A quick examination of literature on evaluation of oscilatory integrals shows that the star is a method developed in 1928 by Filon, which essentially calculates the integral
\int_a^b f(x) sin(\omega x) dx
by approximating f(x) in a polynomial and uses exact formulae for \int_a^b x^i sin(\omega x) dx.
Here is a formal method to calculate the inverse Fourier transform that I explore for implementation:
We want to compute the inverse Fourier transform of f(k), and we assume analytic f = \int_{j=0}^\infty a_j k^j. Then
InverseFourier(f)(x) = \int_{\infty}^\infty f(k) e^{i x k} dk
Now let u = ixk so that dk = du/(ix) and k^j = u^k/(ix)^j, so the sum becomes a series in (ix)^{1} powers and Gamma integrals since
Gamma(t) = \int_0^\infty u^{t1} e^{u} du
This same trick allows us to compute fractional power Fourier inverses, so the problem of calculating Fourier inversion reduces to a Taylor expansion problem since translations in frequency are multiplications by a phase factor.
This works when the function has a global analytic expression since we have the Gamma function. Unfortunately, we will have functions that are not globally expressible as power series and have to implement this Filon’s method. I hope to produce some working code for this next.

3:37 PM (4 hours ago)



Finon’s method is available in Macsyma, so this is a translation task (http://www.cs.berkeley.edu/~fateman/papers/oscillate.pdf)
The following code is the Finon algorithm accepting a function object ‘f’ integrating between 0,1 f(x)sin(mx). It’s described in (http://www.cs.berkeley.edu/~fateman/papers/oscillate.pdf). We’re interested in the inverse fourier transforms
\int_{\infty}^\infty f(k) (cos(x k) dk + i sin(x k)) dk
which modifies the interpolation points and we need the cosine version, both of which are relatively simple from here. The application is to Green’s function of fractional diffusions which are theoretically computed as inverse Fourier transforms of some computable quantities. Then we sample from these and evaluate the Hurst exponents of such sampling in order to assess whether financial time series can be productively described by fractional diffusions.
def finons0(f,m,p):
h=1/(2*p)
k=m
theta=k*h
s2p = 0
for i in range(p):
s2p += f(2*h*i)*sin(2*k*h*i)
s2p = 1/2 * f(1)*sin(k)
s2pm1 = 0
for i in range(p):
s2pm1 += f(h*(2*i1))*sin(k*h*(2*i1))
alf = 1/theta + sin(2*theta)/2/theta**2 – sin(2*theta)/theta**3
bet = 2*((1+cos(theta)**2)/theta**2 – sin(2*theta)/theta**3)
gam = 4*(sin(theta)/theta**3 – cos(theta)/theta**2)
return h*(alf*(f(0)f(1)*cos(k)+bet*s2p+gam*s2pm1)

7:00 PM (1 hour ago)



The FINON METHOD for general inverse Fourier transform is the following:
def finon_zero( x, xi1, xi2):
return 1/(1j*x)*(exp( 1j*x*xi2) – exp(1j*x*xi1))
def finon_one(x, xi1, xi2):
u1= 1j*x*xi1
u2 = 1j*x*xi2
return (1/x**2)*((exp(u2)*u2 – exp(u1)*u1 – exp(u2) + exp(u1))
def finon_two(x,xi1,xi2):
u1= 1j*x*xi1
u2 = 1j*x*xi2
return (1/(1j*x))**3*(u2**2*exp(u2) – u2*exp(u2)+exp(u2) – u1**2*exp(u1)u1*exp(u1) + exp(u1))
def finon_integral( f1, f2, f3, xi1, xi2, xi3, x ):
# solve a 3×3 linear system for coefficients
# for the quadratic passing through f1, f2, f3
A=np.ndarrary([ x1**2, x1, 1, xi2**2, xi2,1, xi3**2, xi3, 1],shape=(3,3))
quad = np.linalg.solve(A,np.array([f1,f2,f3])
res = 0
res += quad[0]*finon_two(x, xi1,xi3)
res += quad[1]*finon_one(x,xi1,xi3)
res += quad[2]*finon_three(x,xi1,xi3)
return res
def inverse_fourier( f, x, xi_start, xi_end, N):
xivals = np.linspace(xi_start, xi_end, N)
ifval = 0
for j in range(N1):
xi1 = xivals[j]
xi3 = xivals[j1]
xi2 = (xi1+xi3)/2.0
f1 = f(xi1)
f2 = f(xi2)
f3 = f(xi3)
ifval += finon_integral(f1, f2, f3, xi1, xi2, xi3, x)
return(ifval)

8:04 PM (16 minutes ago)



Now we return the Green function for fractional diffusions of Mainardi before returning to debug the inverse_fourier implementation. If we take a fractional diffusion equation with timederivative beta and space operator RieszFeller with symbol psi(xi) = xi^alpha exp( i (sign(xi))*theta*pi/2, and let L and F stand for Fourier and Laplace transforms, then the Green’s function satisfies
LF (G)( xi, s) = s^{beta1}/(s^beta + psi(xi))
So we use the fourier_inverse above and the attached Laplace inversion code to obtain G(t,x).
Example:
from talbot import *
N=200
xvals = np.linspace( 1.0, 1.0, num=N)
xivals = np.linspace(10.0, 10.0, num=N)
svals = np.linspace(0.0, 100.0, num=N)
def LFG( xi, s, alf, bet):
return s**(bet1)/(s**bet + psi(xi, alf, theta))
So we need to put all the code together to take inverse Laplace and Fourier transforms of this LFG and then we can sample from the distributions and check their Hurst exponents.

8:08 PM (13 minutes ago)



This entire exercise is useful because of the major hypothesis (intuitively obvious but no consensus at all at least among Wall Street and academic finance workers and scholars) that even though it’s possible to have Hurst exponents deviate from 1/2 for Markov processes, the natural interpretation of such things is that there’s a fractional diffusion (nonlinear or not) governing the underlying continous time process.
Attachments area
Preview attachment talbot.py
talbot.py
Advertisements
Leave a Reply