In [ ]:
import numpy as np
from fgn import fbm, fgn
from hermite import hermite
from main import path_analyse
from crossing_tree import xtree_build
In [ ]:
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
# generator = fbm( N = 2**20+1, H = .75 )
generator = hermite( N = 2**18+1, d = 2, H = .9, K = 2**4, time = True )
generator.set_rnd( np.random.RandomState( ) )
In [ ]:
# X = np.cumsum( np.random.randn( len( T ) ) )
T, X = generator( )
In [ ]:
import numpy as np
from main import path_analyse
from crossing_tree import xtree_build
import scipy.io as io
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
mat = io.loadmat('./output/HRM_3_20-32_0.9000_20150525-144821_AC8A.mat')
X = mat['X'][0]
T = np.arange( len(X), dtype = np.float64 )/ (len(X)-1)
In [ ]:
from numpy.polynomial.hermite_e import hermeval
In [ ]:
gfgn = fgn( 2**18+1, H = (0.5+1)/2, sigma = 1.0 )
gfgn.set_rnd( np.random.RandomState( ) )
In [ ]:
gfgn.reset()
Z = gfgn( )
print Z.mean()
X = np.cumsum( hermeval( Z, c = [0,0,1] ) )[ 15::16 ]
T = np.arange( len(X), dtype = np.float64 )/ (len(X)-1)
In [ ]:
# T, X = generator( )
print np.median( np.abs( np.diff( X ) ) )
print np.std( np.diff( X ) )
print np.subtract( *np.percentile( np.diff( X ), [75, 25] ) )
In [ ]:
delta = np.median( np.abs( np.diff( X ) ) )
Tnk, Xnk, Znk, Vnde, Wnk = xtree_build( T, X, delta = delta )
l = len( Tnk )-5
plt.figure( figsize = ( 16, 9 ) )
plt.plot( T, X/np.max(np.abs( X )), '-', color = 'lightgray' )
# plt.plot( T[mask], X[mask], '-', color = 'orange', linewidth = 2 )
plt.step( Tnk[l-3], Xnk[l-3]/np.max( np.abs( X ) ), '-y>', where = 'pre' )
plt.step( Tnk[l-2], Xnk[l-2]/np.max( np.abs( X ) ), '-g>', where = 'pre' )
plt.step( Tnk[l-1], Xnk[l-1]/np.max( np.abs( X ) ), '-b>', where = 'pre' )
plt.step( Tnk[l-0], Xnk[l-0]/np.max( np.abs( X ) ), '-r>', where = 'pre' )
# plt.yticks( np.arange( -512, 256+1, 64 ) )
plt.grid( which = 'both' )
In [ ]:
Nn, Dnk, Vnde, _ = path_analyse( T, X, delta = delta, max_levels = 16, max_crossings = 18 )
In [ ]:
Dnk
In [ ]:
Djnk = np.reshape( Dnk, (1,)+Dnk.shape)
Njn = np.reshape( Nn, (1,)+Nn.shape)
import self_similarity as ss
ss.pooled( Djnk, Njn, H = None, p = 6, q = 9, h = 4 )
In [ ]:
Z = X - X.mean( )
G_s = np.correlate( Z, Z, mode = 'full' )[ -len( Z ): ] / len( Z )
acf = G_s / G_s[0]
plt.plot( acf )
In [ ]:
%lprun -f xtree_build for i in range( 10 ) : T1nk, X1nk, Z1nk, V1nde, W1nk = xtree_build( T, X, delta = delta )
In [ ]:
# %lprun -f xtree_build
# %timeit -n10
T2nk, X2nk, Z2nk, V2nde, W2nk = xtree_build( T, X, delta = delta )
In [ ]:
print [np.allclose( x1, x2 ) for x1, x2 in zip( T1nk, T2nk ) ]
print [np.allclose( x1, x2 ) for x1, x2 in zip( X1nk, X2nk ) ]
print [np.allclose( x1, x2 ) for x1, x2 in zip( Z1nk, Z2nk ) ]
print [np.allclose( x1, x2 ) for x1, x2 in zip( V1nde, V2nde ) ]
print [np.allclose( x1, x2 ) for x1, x2 in zip( W1nk, W2nk ) ]
In [ ]:
# mask = np.random.uniform( size = len( T ) ) > .99
# mask &= ( ( T > .125 ) & ( T < .475 ) | ( T > .625 ) & ( T < .975 ) )
mask = np.ones( len( T ), np.bool )
In [ ]:
# %timeit -n10
delta = np.median( np.abs( np.diff( X[mask] ) ) )
T1nk, X1nk, Z1nk, V1nde, W1nk = xtree_build_new( T[mask], X[mask],
shift = .15, delta = delta )
T2nk, X2nk, Z2nk, V2nde, W2nk = xtree_build_new( T[mask], X[mask],
shift = -.15, delta = delta )
In [ ]:
print np.unique( Z1nk[6], return_counts = True )
print np.unique( Z2nk[6], return_counts = True )
In [ ]:
plt.figure( figsize = ( 16, 9 ) )
plt.plot( T, X, '-', color = 'lightgray' )
plt.plot( T[mask], X[mask], '-', color = 'orange', linewidth = 2 )
plt.step( T1nk[10], X1nk[10], '-b', where = 'post' )
plt.step( T1nk[12], X1nk[12], '-r', where = 'post' )
plt.step( T2nk[10], X2nk[10], '-g', where = 'post' )
plt.step( T2nk[12], X2nk[12], '-k', where = 'post' )
# plt.yticks( np.arange( -512, 256+1, 64 ) )
plt.grid( which = 'both' )
In [ ]:
max_levels, max_crossings = 7, 20
In [ ]:
T, X = generator( )
# delta = np.std( np.diff( X ) )
In [ ]:
delta = np.subtract( *np.percentile( np.diff( X ), [ 75, 25 ] ) )
%lprun -f xtree_integer_crossings for i in range( 10 ) : Tnk, XTnk, Znk, Vnk, Wnk = xtree_build( T, X, delta = delta )
In [ ]:
# Nn = np.array( [ len( Tk ) - 1 for Tk in Tnk ],
# dtype = np.float ).reshape( ( len( Tnk ), -1 ) )
Nn = np.zeros( ( max_levels + 1, 1 ), dtype = np.int )
for n, Tk in enumerate( Tnk, 0 ) :
n = max_levels if n > max_levels else n
Nn[ n ] += len( Tk ) - 1
Dnk = np.zeros( ( max_levels + 1, max_crossings // 2 ), dtype = np.int )
for n, Zk in enumerate( Znk[ 1: ], 0 ) :
n = max_levels if n > max_levels else n
Z_count, Z_freq = np.unique( Zk, return_counts = True )
Z_count = np.minimum( Z_count, max_crossings )
mask = ( Z_count < max_crossings )
Dnk[ n, Z_count[ mask ] // 2 - 1 ] += Z_freq[ mask ]
Dnk[ n, max_crossings // 2 - 1 ] += np.sum( Z_freq[ ~mask ] )
Vnde = np.zeros( ( max_levels + 1, 2, 2 ), dtype = np.int )
for n, Vk in enumerate( Vnk[ 1: ], 0 ) :
n = max_levels if n > max_levels else n
Vnde[ n, 0 ] += np.sum( Vk[ Vk[ :, 2 ] < 0 ], axis = 0 )[:2]
Vnde[ n, 1 ] += np.sum( Vk[ Vk[ :, 2 ] > 0 ], axis = 0 )[:2]
## Reshape properly
Djnk = np.reshape( Dnk, (1,)+Dnk.shape)
Njn = np.reshape( Nn, (1,)+Nn.shape)
In [ ]:
Djnk
In [ ]:
print Njn[:2,1:,0]
print np.sum( Djnk[:2], axis = 2 )
In [ ]:
import self_similarity as ss
ss.pooled( Djnk, Njn, H = None, p = 4, q = 12, h = 5 )
In [ ]:
np.sum( Djnk, axis = ( 0, ) )
In [ ]:
print Dnk[0,:]#.dot( np.arange( 2, max_crossings+2, 2 ) )
print Nn
print Znk[1]
In [ ]:
Vnd = np.sum( Vnde, axis = 2, dtype = np.float ).reshape( VDn.shape[:2] + ( 1, ) )
Dn = np.sum( Dnk, axis = 1, dtype = np.float ).reshape( Dnk.shape[:1] + ( 1, ) )
EXnde = Vnde / Vnd
CDnk = Dnk / Dn
print CDnk
print EXnde
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# dat = np.load( './output/fbm_iqr_20150429-072004_18_0.50_200.npz' )
# dat = np.load( './output/fbm_iqr_20150429-072224_18_0.72_200.npz' )
# dat = np.load( './output/fbm_iqr_20150429-072516_18_0.95_200.npz' )
# dat = np.load( './output/fbm_iqr_20150429-083530_16_0.7250_200.npz' )
# dat = np.load( './output/fbm_iqr_20150429-083459_16_0.5000_200.npz' )
# dat = np.load( './output/fbm_iqr_20150429-090509_21_0.5000_1000.npz' )
dat = np.load( './output/fbm_iqr_20150429-150719_21_0.6500_1000.npz' )
Njn = dat['Njn']
Djnk = dat['Djnk']
Vjnde = dat['Vjnde']
Djn = np.sum( Djnk, axis = 2, dtype = np.float ).reshape( Djnk.shape[:2] + ( 1, ) )
Vjnd = np.sum( Vjnde, axis = 3, dtype = np.float ).reshape( Vjnde.shape[:3] + ( 1, ) )
Phat_jnk = Djnk / Djn
Qhat_jnde = Vjnde / Vjnd
Phat_nk = Phat_jnk.mean( axis = 0, dtype = np.float64 )
Qhat_nde = Qhat_jnde.mean( axis = 0, dtype = np.float64 )
sdPhat_nk = Phat_jnk.std( axis = 0, dtype = np.float64 )
sdQhat_nde = Qhat_jnde.std( axis = 0, dtype = np.float64 )
In [ ]:
# plt.figure( figsize = ( 16, 9 ) )
# for Phat_k in Phat_nk[:5] :
# plt.plot( np.log( Phat_k ) )
H = 0.65
theta = 2**( 1 - 1/H )
theoretical = [ theta*(1-theta)**k for k in range( 0, 10 ) ]
level = 5
plt.figure( figsize = ( 16, 9 ) )
plt.subplot(221)
plt.hist( Phat_jnk[:,level,0], color = 'gray' )
plt.title( "Simulated distribution: res. $\\delta 2^{%d}$" % level )
plt.xlabel( "$\\hat{\\mathbb{P}}(Z = 2)$" )
plt.axvline(theoretical[0], color='y', linewidth=2)
plt.axvline(Phat_nk[level,0], color='r', linewidth=2)
plt.axvline(Phat_nk[level,0] - sdPhat_nk[level,0], color='pink', linestyle='dashed', linewidth=2)
plt.axvline(Phat_nk[level,0] + sdPhat_nk[level,0], color='pink', linestyle='dashed', linewidth=2)
plt.subplot(222)
plt.hist( Phat_jnk[:,level,1], color = 'gray' )
plt.title( "Simulated distribution: res. $\\delta 2^{%d}$" % level )
plt.xlabel( "$\\hat{\\mathbb{P}}(Z = 4)$" )
plt.axvline(theoretical[1], color='y', linewidth=2)
plt.axvline(Phat_nk[level,1], color='r', linewidth=2)
plt.axvline(Phat_nk[level,1] - sdPhat_nk[level,1], color='pink', linestyle='dashed', linewidth=2)
plt.axvline(Phat_nk[level,1] + sdPhat_nk[level,1], color='pink', linestyle='dashed', linewidth=2)
plt.subplot(223)
plt.hist( Phat_jnk[:,level,2], color = 'gray' )
plt.title( "Simulated distribution: res. $\\delta 2^{%d}$" % level )
plt.xlabel( "$\\hat{\\mathbb{P}}(Z = 6)$" )
plt.axvline(theoretical[2], color='y', linewidth=2)
plt.axvline(Phat_nk[level,2], color='r', linewidth=2)
plt.axvline(Phat_nk[level,2] - sdPhat_nk[level,2], color='pink', linestyle='dashed', linewidth=2)
plt.axvline(Phat_nk[level,2] + sdPhat_nk[level,2], color='pink', linestyle='dashed', linewidth=2)
plt.subplot(224)
plt.hist( Phat_jnk[:,level,6], color = 'gray' )
plt.title( "Simulated distribution: res. $\\delta 2^{%d}$" % level )
plt.xlabel( "$\\hat{\\mathbb{P}}(Z = 14)$" )
plt.axvline(theoretical[6], color='y', linewidth=2)
plt.axvline(Phat_nk[level,6], color='r', linewidth=2)
plt.axvline(Phat_nk[level,6] - sdPhat_nk[level,6], color='pink', linestyle='dashed', linewidth=2)
plt.axvline(Phat_nk[level,6] + sdPhat_nk[level,6], color='pink', linestyle='dashed', linewidth=2)
In [ ]:
print Djnk.shape
In [ ]:
print Njn.mean( axis = 0, dtype = np.float64 )
print Njn.std( axis = 0, dtype = np.float64 )
print Phat_jnk.mean( axis = 0, dtype = np.float64 )
print Phat_jnk.std( axis = 0, dtype = np.float64 )
print Qhat_jnde.mean( axis = 0, dtype = np.float64 )
print Qhat_jnde.std( axis = 0, dtype = np.float64 )
In [ ]:
theta = 2**(1-1/0.50)
print theta
print 1 / np.sqrt( 2 / (1-theta) )
Here I confirm that my implementation is correct.
In [ ]:
Z = ( X - X[ 0 ] ) / delta
In [ ]:
%lprun -f xtree_integer_crossings for i in range( 10 ) : lht0, lhp0 = xtree_integer_crossings( T, Z )
# %timeit -n10 lht0, lhp0 = xtree_integer_crossings( T, Z )
In [ ]:
%lprun -f f_get_w_int for i in range( 1 ) : mht0, mhp0, _, _ = f_get_w_int( T, Z )
# %timeit -n10 mht0, mhp0, _, _ = f_get_w_int( T, Z )
In [ ]:
print np.allclose( mht0, lht0 )
print np.allclose( mhp0, lhp0 )
The following section concerns the correctness of the FGN generator.
In [ ]:
import numpy as np
from numpy.fft import fft, fft
rnd = np.random.RandomState( 123 )
N = 2**16 + 1
sigma = 1.0 ; H = (.99+1)/2
In [ ]:
R = np.arange( N, dtype = np.float )
R = sigma * sigma * .5 * (
np.abs( R - 1 ) ** ( 2.0 * H )
+ np.abs( R + 1 ) ** ( 2.0 * H )
- 2 * np.abs( R ) ** ( 2.0 * H ) )
In [ ]:
plt.plot( R )
In [ ]:
Z = np.real( fft( np.append( R, R[::-1][1:-1] ) ) )
In [ ]:
print len( R )
print len( Z )
print 2 * N - 2
In [ ]:
Z = np.sqrt( np.maximum( Z, 0.0 ) / ( 2 * N - 2 ) )
In [ ]:
W = rnd.randn( 2 * N - 2 ) + rnd.randn( 2 * N - 2 ) * 1j
W *= Z
C = fft( W )
u, v = ( np.real( C[ :N ] ), np.imag( C[ :N ] ) )
In [ ]:
u,v
In [ ]:
plt.figure( figsize = (16,9) )
plt.subplot(121)
plt.plot( u )
plt.subplot(122)
plt.plot( v )
In [ ]:
z = v - v.mean()
r = np.correlate( z,z, mode = 'full' )[-len(z):] / (len(z)*z.var())
In [ ]:
plt.plot( R, r )
In [ ]: