Merge pull request #31 from LRydin/examples

Adding examples for reproducibility.

.github/workflows/CI.yml

@@ -12,14 +12,14 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: ['3.5', '3.6', '3.7', '3.8', '3.9', '3.10', pypy-3.6, pypy-3.7, pypy-3.8]
+        python-version: ['3.6', '3.7', '3.8', '3.9', '3.10', pypy-3.7, pypy-3.8, pypy-3.9]
         exclude:
-          - os: macos-latest
-            python-version: pypy-3.6
           - os: macos-latest
             python-version: pypy-3.7
           - os: macos-latest
             python-version: pypy-3.8
+          - os: macos-latest
+            python-version: pypy-3.9
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2019-2021 Leonardo Rydin Gorjão
+Copyright (c) 2019-2022 Leonardo Rydin Gorjão
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

docs/license.rst

@@ -3,7 +3,7 @@ License
 
 MIT License
 
-Copyright (c) 2019-2021 Leonardo Rydin Gorjão
+Copyright (c) 2019-2022 Leonardo Rydin Gorjão
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

examples/fig01.py

@@ -0,0 +1,146 @@
+# created by Leonardo Rydin Gorjão. Most python libraries are standard (e.g. via
+# Anaconda). If TeX is not present in the system comment out lines 12 to 15.
+
+import numpy as np
+
+# to install MFDFA just run pip install MFDFA
+from MFDFA import MFDFA
+from MFDFA import fgn
+
+import matplotlib
+import matplotlib.pyplot as plt
+matplotlib.rcParams['pgf.texsystem'] = 'pdflatex'
+matplotlib.rcParams.update({'font.family': 'serif', 'font.size': 18,
+    'axes.labelsize': 20,'axes.titlesize': 24, 'figure.titlesize' : 28})
+matplotlib.rcParams['text.usetex'] = True
+
+colours = ['#1b9e77','#d95f02','#7570b3']
+
+################################## fOU process #################################
+
+# %% Lets take a fractional Ornstein–Uhlenbeck process Xₜ with a time-dependent
+# diffusion or volatility σₜ, some drift of mean reverting term θ(Xₜ), and
+# fractional Gaussian noise with Hurst exponent H
+#
+#                           dXₜ = - θ(Xₜ)Xₜdt + σₜdBᴴₜ.
+
+# integration time and time sampling
+t_final = 100
+delta_t = .001
+
+# The parameters θ and σ
+theta = 1
+sigma = .5
+
+# The time array of the trajectory
+time = np.arange(0, t_final, delta_t)
+
+# Initialise the array y
+X1 = np.zeros(time.size)
+X2 = np.zeros(time.size)
+X3 = np.zeros(time.size)
+
+# Generate the fractional Gaussian noise
+dw1 = (t_final ** .3) * fgn(time.size, H=.3)
+dw2 = (t_final ** .5) * fgn(time.size, H=.5)
+dw3 = (t_final ** .7) * fgn(time.size, H=.7)
+
+# Integrate the process
+for i in range(1, time.size):
+    X1[i] = X1[i-1] - theta*X1[i-1] * delta_t + sigma * dw1[i]
+    X2[i] = X2[i-1] - theta*X2[i-1] * delta_t + sigma * dw2[i]
+    X3[i] = X3[i-1] - theta*X3[i-1] * delta_t + sigma * dw3[i]
+
+# %% MFDFA
+# Lag s from 3 to 1000 datapoints
+lag = np.unique(np.logspace(.6, 3, 118).astype(int))
+
+# q power variations, removing the 0 power
+q_list = np.linspace(-10, 10, 41)
+q_list = q_list[q_list!=.0]
+
+lag, dfa1 = MFDFA(X1, lag, q=q_list, order=1)
+lag, dfa2 = MFDFA(X2, lag, q=q_list, order=1)
+lag, dfa3 = MFDFA(X3, lag, q=q_list, order=1)
+
+# %% ################################## FIG01 ##################################
+fig, ax = plt.subplots(1,2, figsize=(16,4))
+axi1 = fig.add_axes([.31, .2, .18, .38])
+
+ax[0].loglog(lag, dfa1[...,23],'o', markersize=6, color = colours[0],
+    markerfacecolor='none', label=r'$H=0.3, q=2$')
+ax[0].loglog(lag, dfa2[...,23],'s', markersize=6, color = colours[1],
+    markerfacecolor='none', label=r'$H=0.5, q=2$')
+ax[0].loglog(lag, dfa3[...,23],'D', markersize=6, color = colours[2],
+    markerfacecolor='none', label=r'$H=0.7, q=2$')
+
+ax[0].set_ylim([5e-4,4e1])
+
+axi1.loglog(lag[10:], .02 * lag[10:]**1.3,'--', color = 'black', alpha= .5)
+axi1.loglog(lag[10:], .0002 * lag[10:]**1.3,'--', color = 'black', alpha= .5)
+axi1.loglog(lag[10::2], dfa1[10::2,[0]],  '^', markersize=6, color=colours[1],
+    markerfacecolor='none', label=r'$q=-\!\!10$')
+axi1.loglog(lag[10::2], dfa1[10::2,[17]], 'v', markersize=6, color=colours[1],
+    markerfacecolor='none', label=r'$q=-\!\!2$')
+axi1.loglog(lag[10::2], dfa1[10::2,[23]], '<', markersize=6, color=colours[1],
+    markerfacecolor='none', label=r'$q=2$')
+axi1.loglog(lag[10::2], dfa1[10::2,[39]], '>', markersize=6, color=colours[1],
+    markerfacecolor='none', label=r'$q=10$')
+
+axi1.set_xticks([])
+fig.text(.315, .49, r'$H=0.3$', fontsize=20)
+fig.text(.365, .225, r'$q=-10,-2,2,10$', fontsize=18)
+
+slopes1 = np.polynomial.polynomial.polyfit(np.log(lag)[5:45],
+    np.log(dfa1)[5:45],1)[1]-1
+slopes2 = np.polynomial.polynomial.polyfit(np.log(lag)[5:45],
+    np.log(dfa2)[5:45],1)[1]-1
+slopes3 = np.polynomial.polynomial.polyfit(np.log(lag)[5:45],
+    np.log(dfa3)[5:45],1)[1]-1
+
+ax[1].plot(q_list, q_list * slopes1 - 1, 'o', markersize=6, color=colours[0],
+    markerfacecolor='none', label=r'$H=0.3$')
+ax[1].plot(q_list, q_list * slopes2 - 1, 's', markersize=6, color=colours[1],
+    markerfacecolor='none', label=r'$H=0.5$')
+ax[1].plot(q_list, q_list * slopes3 - 1, 'D', markersize=6, color=colours[2],
+    markerfacecolor='none', label=r'$H=0.7$')
+
+ax[1].plot(q_list, q_list * .3 - 1, '-', color='black', alpha=.7)
+ax[1].plot(q_list, q_list * .5 - 1, '-', color='black', alpha=.7)
+ax[1].plot(q_list, q_list * .7 - 1, '-', color='black', alpha=.7)
+
+axi2 = fig.add_axes([.802, .2, .18, .38])
+axi2.plot(q_list, slopes1, 'o', markersize=6, color=colours[0],
+    markerfacecolor='none', label=r'$H=0.3$')
+axi2.plot(q_list, slopes2, 's', markersize=6, color=colours[1],
+    markerfacecolor='none', label=r'$H=0.5$')
+axi2.plot(q_list, slopes3, 'D', markersize=6, color=colours[2],
+    markerfacecolor='none', label=r'$H=0.7$')
+
+axi2.plot(q_list, np.ones_like(q_list) * .3, '-', color='black', alpha=.7)
+axi2.plot(q_list, np.ones_like(q_list) * .5, '-', color='black', alpha=.7)
+axi2.plot(q_list, np.ones_like(q_list) * .7, '-', color='black', alpha=.7)
+
+axi2.set_ylim([.1,.9])
+axi2.set_yticks([.3,.5,.7])
+axi2.set_xticks([])
+axi2.set_ylabel(r'$h(q)$',labelpad=3,fontsize=24)
+
+ax[0].set_ylabel(r'$F_q(s)$',labelpad=3,fontsize=24)
+ax[0].set_xlabel(r'segment size $s$',labelpad=3,fontsize=24)
+
+ax[1].set_ylabel(r'$\tau(q)$',labelpad=3,fontsize=24)
+ax[1].set_xlabel(r'$q$',labelpad=3,fontsize=24)
+
+ax[0].legend(fontsize=20, loc=2, handlelength=1.1, columnspacing=.6,
+    handletextpad=.2)
+ax[1].legend(fontsize=20, loc=2, handlelength=1.1, columnspacing=.6,
+    handletextpad=.2)
+
+fig.text(.005, .92, r'a)', fontsize=28)
+fig.text(.505, .92, r'b)', fontsize=28)
+
+fig.subplots_adjust(left=.07, bottom=.17, right=.99, top=.99, hspace=.06,
+    wspace=.15)
+
+# fig.savefig('fig01.pdf', trasparent=True)

examples/fig02.py

@@ -0,0 +1,128 @@
+# created by Leonardo Rydin Gorjão. Most python libraries are standard (e.g. via
+# Anaconda). If TeX is not present in the system comment out lines 13 to 16.
+
+import numpy as np
+from scipy.stats import levy_stable
+
+# to install MFDFA just run pip install MFDFA
+from MFDFA import MFDFA
+from MFDFA import fgn
+
+import matplotlib
+import matplotlib.pyplot as plt
+matplotlib.rcParams['pgf.texsystem'] = 'pdflatex'
+matplotlib.rcParams.update({'font.family': 'serif', 'font.size': 18,
+    'axes.labelsize': 20,'axes.titlesize': 24, 'figure.titlesize' : 28})
+matplotlib.rcParams['text.usetex'] = True
+
+colours = ['#1b9e77','#d95f02','#7570b3']
+
+# %% ######################### Lévy distributions ##############################
+
+# integration time and time sampling
+t_final = 1000
+delta_t = .001
+
+# The time array of the trajectory
+time = np.arange(0, t_final, delta_t)
+
+alpha = 1.75
+X4 = levy_stable.rvs(alpha=alpha, beta=0, size=time.size)
+
+alpha = 1.25
+X5 = levy_stable.rvs(alpha=alpha, beta=0, size=time.size)
+
+alpha = .75
+X6 = levy_stable.rvs(alpha=alpha, beta=0, size=time.size)
+
+# %%
+lag = np.unique(np.logspace(0.6, 3, 118).astype(int))
+
+q_list = np.linspace(-10,10,41)
+q_list = q_list[q_list!=0.0]
+lag, dfa4 =  MFDFA(X4, lag, q=q_list, order=3)
+lag, dfa5 = MFDFA(X5, lag, q=q_list, order=3)
+lag, dfa6 = MFDFA(X6, lag, q=q_list, order=3)
+
+# %% ################################## FIG02 ##################################
+fig, ax = plt.subplots(1,2, figsize=(16,4))
+
+ax[0].loglog(lag[9::2], dfa5[9::2,[0]], '^', markersize=6,
+    markerfacecolor='none', color=colours[1], label=r'$q=-10$')
+ax[0].loglog(lag[9::2], dfa5[9::2,[23]],'<', markersize=6,
+    markerfacecolor='none', color=colours[1], label=r'$q=2$')
+ax[0].loglog(lag[9::2], dfa5[9::2,[10]],'D', markersize=6,
+    markerfacecolor='none', color=colours[1], label=r'$q=-5$')
+ax[0].loglog(lag[9::2], dfa5[9::2,[29]],'H', markersize=6,
+    markerfacecolor='none', color=colours[1], label=r'$q=5$')
+ax[0].loglog(lag[9::2], dfa5[9::2,[17]],'v', markersize=6,
+    markerfacecolor='none', color=colours[1], label=r'$q=-2$')
+ax[0].loglog(lag[9::2], dfa5[9::2,[39]],'>', markersize=6,
+    markerfacecolor='none', color=colours[1], label=r'$q=10$')
+
+ax[0].loglog(lag[9::2], 0.03 * lag[9::2]**(1/1.25),'--', color='black')
+ax[0].loglog(lag[9::2], 4e4 * lag[9::2]**(1/10),'--', color='black')
+
+
+ax[0].set_ylabel(r'$F_q(s)$', labelpad=7, fontsize=24)
+ax[0].set_xlabel(r'segment size $s$', labelpad=3, fontsize=24)
+ax[0].set_ylim([9e-2, 2e5])
+fig.text(.18, .22, r'$\alpha = 1.25$', fontsize=22)
+
+slopes4 = np.polynomial.polynomial.polyfit(np.log(lag)[25:75],
+    np.log(dfa4)[25:75],1)[1]
+slopes5 = np.polynomial.polynomial.polyfit(np.log(lag)[25:75],
+    np.log(dfa5)[25:75],1)[1]
+slopes6 = np.polynomial.polynomial.polyfit(np.log(lag)[25:75],
+    np.log(dfa6)[25:75],1)[1]
+
+ax[1].plot(q_list, slopes4, 'o', markersize=9, markerfacecolor='none',
+    color=colours[0], label=r'$\alpha=1.75$')
+ax[1].plot(q_list, slopes5, 's', markersize=9, markerfacecolor='none',
+    color=colours[1], label=r'$\alpha=1.25$')
+ax[1].plot(q_list, slopes6, 'D', markersize=9, markerfacecolor='none',
+    color=colours[2], label=r'$\alpha=0.75$')
+
+ax[1].plot(q_list[:20], np.ones(q_list[:20].size)*(1/1.75),'-', color='black',
+    alpha=.7)
+ax[1].plot(q_list[:20], np.ones(q_list[:20].size)*(1/1.25),'-', color='black',
+    alpha=.7)
+ax[1].plot(q_list[:20], np.ones(q_list[:20].size)*(1/0.75),'-', color='black',
+    alpha=.7)
+
+ax[1].set_ylim([.001, 1.49])
+ax[1].set_ylabel(r'$h(q)$',labelpad=3, fontsize=24)
+ax[1].set_xlabel(r'$q$',labelpad=3, fontsize=24)
+
+axi2 = fig.add_axes([.85, .6, .135, .37])
+axi2.plot(q_list, q_list * slopes4 - 1,'o', markersize=6, color=colours[0],
+    markerfacecolor='none')
+axi2.plot(q_list, q_list * slopes5 - 1,'s', markersize=6, color=colours[1],
+    markerfacecolor='none')
+axi2.plot(q_list, q_list * slopes6 - 1,'D', markersize=6, color=colours[2],
+    markerfacecolor='none')
+
+axi2.set_ylim([-16, 1.2])
+axi2.set_yticks([-12, -6, 0])
+axi2.set_xlabel(r'$q$', labelpad=3, fontsize=24)
+axi2.set_ylabel(r'$\tau(q)$', labelpad=-5, fontsize=24)
+
+ax[0].legend(loc=4, handletextpad=.3, handlelength=.5, ncol=3,
+    columnspacing=.65)
+ax[1].legend(loc=3, handletextpad=.3, handlelength=.5, ncol=3,
+    columnspacing=.65)
+
+locmaj = matplotlib.ticker.LogLocator(base=10.0, subs=(1.0, ), numticks=100)
+ax[0].yaxis.set_major_locator(locmaj)
+
+locmin = matplotlib.ticker.LogLocator(base=10.0, subs=np.arange(2, 10) * .1,
+                                      numticks=100)
+ax[0].yaxis.set_minor_locator(locmin)
+ax[0].yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
+
+fig.text(.005, .92, r'a)',fontsize=28)
+fig.text(.505, .92, r'b)',fontsize=28)
+
+fig.subplots_adjust(left=.07, bottom=.17, right=.99, top=.99, hspace=.06,
+wspace=.15)
+# fig.savefig('fig02.pdf', trasparent=True)