HaarNoInterpolateUncertainty[timeraw_,xv_,calibcoord_,droppointslowpoints_,droppointshighpoints_]:=Module[{},{{dim}=Dimensions[xv],(*dimmax=dim; use this when fac=1*),dimmax=dim,(*a specific value for dimmax is only needed when fac=2, it is the maximum dimension to be encountered*),

(*calib=2.*),vecc=Table[10^-10 jjj,{jjj,1,dim}],(*this avoids having two identical time coordinates for different estimates of the same quantity at the same time*),time=timeraw+vecc,ϵmin=0.25,ϵmax=1-ϵmin,(*this is the range of reasonable values for estimating the flucutations, default is ϵmin=0.25, it means that the centre is asssumed to be between a fraction ϵmin and (1-ϵmin) between the extremes (rather than halfway*),calib=2.,

(*This is the "Haar calibration" i.e. the value that approxiamtely maps the Haar fluctuations onto difference and tendency fluctuations, the default value is 2., note that the minimum significant difference between successive time intervals is 10^-4, closer intervals will be put in smallest bin.  The largest is 10^6.*),

mindelta=0.1,δtimemax=time[[dim]]-time[[1]],

qmax=21,delq=0.1,qeps=0.0001`,Qtab=Table[q,{q,qeps,delq qmax+qeps-delq,delq}],

(*with these parameters we calucate the moments every delq=0.1 from 0 to q=2 (qmax values in all).  qeps is added to avoid exponenting to a zero power - a problem if the data contains zero values.  calib is an overall multiplier used only for the graphics output grRMS, grq1 that takes into account the calibration with respect to difference and anomaly fluctuations*)

points2=10^6,points1=IntegerPart[20. Log[10,(points2-0.001)/fac]],
δδt30={2,4,6,8,10,12,14,16,18,20,22},maxpoint=2 IntegerPart[dimmax/(2 fac)],δδt2=Flatten[Append[δδt30,Table[2  IntegerPart[((10.^(jj/20.)))/2],{jj,28,points1}]]],δδt=Select[δδt2,#<(maxpoint+0.001)&],points=Dimensions[δδt][[1]],δδt[[points]]=If[δδt[[points]]==dim,δδt[[points]]=δδt[[points]]-2,δδt[[points]]],

SS=Running[xv-Mean[xv]],δT=Table[0.,{kindex,1,200},{q,1,qmax}],
count=Table[10.^-15,{kindex,1,200}],

Do[{δt=δδt[[jj]],Sδt20=RotateLeft[SS,δt]-SS,Sδt10=RotateLeft[SS,δt/2]-SS,Sδt21=Sδt20-Sδt10,timeδt20=Abs[RotateLeft[time,δt]-time],timeδt10=Abs[RotateLeft[time,δt/2]-time],ϵ=timeδt10/timeδt20,(*we only select those estimates for which the "middle" point fraction ϵ is between a fraction ϵmin and (1-ϵmin) of the total separation (ideally it should be halfway; ϵ=1/2, here we use ϵmin=1/4 so that the centre point is within 1/4 of the true centre*),δRx= 2. Abs[Sδt10-Sδt21]/δt,(*note we divide by the number of contributions, not the time*),ϵS=Select[Transpose[Take[{timeδt20,δRx,ϵ},All,{1,dim-δt}]],ϵmax>#[[3]]>ϵmin&],goodpairs=Select[ϵS[[All,Range[1,2]]],#[[1]]>mindelta&],dimgood=Dimensions[goodpairs][[1]],Do[{index=80.+IntegerPart[20.  Log[10,goodpairs[[k,1]]]];index=Max[1,index];index=Min[200,index];δT[[index,All]]=δT[[index,All]]+goodpairs[[k,2]]^Qtab,count[[index]]=count[[index]]+1.},{k,1,dimgood}]},{jj,1,points}],

(*reversing just duplicates stats, it was not used*),

Print[" the number of flucutations used in the analysis is = ",Total[count]],

coordint=Table[(jj-80)/20.,{jj,0,199}],momsall=Log[10,10^-15+δT/count],
data=Select[Transpose[{coordint,momsall[[All,11]]}],#[[2]]>-5.&],calibnointer=0.,(*calibnointer=Log[10,δtimemax]-Max[data[[All,1]]]*),

moms=Table[Select[momsall[[All,kk]],#>-5.&],{kk,1,21}],dataq1=Select[Transpose[{coordint+calibnointer+Log[10,calibcoord],momsall[[All,11]]}],#[[2]]>-5.&],(*this lower cutoff is to avoid problems with spuriously low values.  Input data should not be too small.  If the typical series values are 10^-3 or lower, then it's bet to pre-multiply the series by a factor 10^3 or so so that typical values are of order unity*),

coordcalib=Transpose[dataq1][[1]],

dataRMS=Select[Transpose[{coordint+calibnointer+Log[10,calibcoord],momsall[[All,21]]/2.}],#[[2]]>-5.&],
grRMSsmoothnoint=ListPlot[Transpose[{Transpose[dataRMS][[1]],Log[10,calib]+Log[10,smooth1D[10.^Transpose[dataRMS][[2]]]]}],Joined->True,PlotStyle->Hue[0.9]],grq1smoothnoint=ListPlot[Transpose[{Transpose[dataq1][[1]],Log[10,calib]+Log[10,smooth1D[10.^Transpose[dataq1][[2]]]]}],Joined->True,PlotStyle->Hue[0.6]],

{points}=Dimensions[coordcalib],
Print[" The total number of points is ",points],
droppointslow=Max[1,1+droppointslowpoints],
droppointshigh=Max[droppointslow+2,points-droppointshighpoints],


(*droppointslow=1,droppointshigh=points, these indicate the range of data used in the fits  All the code below is purely for fitting, it is not strictly needed.  In this default mode all points between the smallest lag 1 and largest available lag (=points) are taken for regressions.  Indeed, the fitted parameters are only valid if the range of fitting is pretty linear on a log-log plot.  Often it is best to make the fits after viewing the plots*),


Print[" Fits are over the range   ", 10^coordcalib[[droppointslow]],"  to  ",10^coordcalib[[droppointshigh]], "  in the specified units"],
(*the calibrated coordinate*),

{inter,hqneg}=Transpose[Table[aa=Fit[Take[Transpose[{coordcalib,moms[[q]]}],{droppointslow,droppointshigh},{1,2}],{1,x1},x1];{aa[[1]],aa[[2]][[1]]},{q,1,qmax}]],Do[grref[i]=Plot[Log[10,calib]+inter[[i]]+hqneg[[i]] x1,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],{i,1,qmax}],Do[gr[i]=ListPlot[Transpose[{coordcalib,Log[10,calib]+moms[[i]]}],Joined->True,PlotRange->All,PlotStyle->{Hue[(i/qmax)],Thickness[0.005]},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)S(Δt)"},BaseStyle->{24,FontFamily->"Times"}],{i,1,qmax}],Do[grpts[i]=ListPlot[Transpose[{coordcalib,Log[10,calib]+moms[[i]]}],PlotRange->All,PlotStyle->{PointSize[0.008],Hue[(i/qmax)]}],{i,1,qmax}],grout=Show[Table[gr[i],{i,1,qmax}],Table[grpts[i],{i,1,qmax}],Table[grref[i],{i,1,qmax}]],psi=Table[{Qtab[[q]],hqneg[[q]]},{q,1,qmax}],grpsi=ListPlot[psi,Joined->True,PlotStyle->Hue[0.75`],AxesLabel->{"q","ξ(q)"},BaseStyle->{24,FontFamily->"Times"}],

ft0=Fit[Flatten[Append[{{{0,0.`}}},Take[psi,{1,2},{1,2}]],1],{1,xxy,(xxy^2)},xxy],psiprime0=ft0[[1]]+ft0[[2]][[1]],ft=Fit[Take[psi,{9,13},{1,2}],{1,xxy1,(xxy1^2)},xxy1],C1psi=ft[[1]]-ft[[3]][[1]],Hpsi=ft/.xxy1->1,Print["H=   ",Hpsi],Print["the RMS exponent= ξ(2)/2= ",psi[[21,2]]/2.,"  the difference with H is K(2)/2 = ",psi[[11,2]]-psi[[21,2]]/2.],Print["derivative C1=  ",C1psi],Kppq1=-2 ft[[3]][[1]],αpsi=Kppq1/C1psi,Print["α from psipp(1) =  ",αpsi],

(*The above are for regressions, parameter estimates over the range droppointslow,droppointshigh*),

C1data=Transpose[{coordcalib,Log[10,smooth1D[10.^(moms[[11]])]]-5.Log[10,smooth1D[10.^(moms[[12]]-moms[[10]])]]}],

grC1=ListPlot[C1data,Joined->True,PlotStyle->Hue[0.99],BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)((Δt)^(Subscript[C, 1]))"}],
grinfo=grC1,

aa=Fit[Take[C1data,{droppointslow,droppointshigh},All],{1,x1},x1],
{interC1,hqnegC1}={aa[[1]],aa[[2]][[1]]},grC1theory=Plot[interC1+hqnegC1 x1,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],
Print["C1 graph slope method =",hqnegC1],
grC1theorydata=Show[grC1,grC1theory,PlotRange->All],

grRMS=ListPlot[Transpose[{coordcalib,Log[10,calib]+moms[[21]]/2.}],Joined->True,PlotStyle->Hue[0.9],BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 2])(Δt)^(1/2)"}],grRMSsmooth=ListPlot[Transpose[{coordcalib,Log[10,calib]+Log[10,smooth1D[10^moms[[21]]]]/2.}],Joined->True,PlotStyle->Hue[0.9],BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 2])(Δt)^(1/2)"}],grcorsmooth=ListPlot[Transpose[{coordcalib,Log[10,calib]+Log[10,smooth1D[10^moms[[21]]]]}],Joined->True,PlotStyle->Hue[0.1],BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 2])(Δt)"}],
grboxsmooth=ListPlot[Transpose[{coordcalib,Log[10,calib]+Log[10,smooth1D[10^moms[[1]]]]}],Joined->True,PlotStyle->Hue[0.4],BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 0])(Δt)"}],
        grq1smooth=ListPlot[Transpose[{coordcalib,Log[10,calib]+Log[10,smooth1D[10^moms[[11]]]]}],Joined->True,PlotStyle->Hue[0.6],BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 1])(Δt)"}],
grq1=ListPlot[Transpose[{coordcalib,Log[10,calib]+moms[[11]]}],Joined->True,PlotStyle->Hue[0.6],BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 1])(Δt)"}],
grtheoryq1=Plot[Log[10,calib]+inter[[11]]+hqneg[[11]] x1,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],
grtheorydataq1=Show[grq1smooth,grtheoryq1],
grtheoryRMS=Plot[Log[10,calib]+inter[[21]]/2+hqneg[[21]] x1/2,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],
grtheorydataRMS=Show[grRMSsmooth,grtheoryRMS,PlotRange->All],

grallsmooth=Show[grq1smooth,grRMSsmooth,PlotRange->All],grall=Show[grq1,grRMS,PlotRange->All],gralldatatheory=Show[grtheorydataq1,grtheorydataRMS,PlotRange->All],

Kppdata=Transpose[{coordcalib,-25.Log[10,smooth1D[10.^(moms[[13]]+moms[[9]]-2 moms[[11]])]]}],
(*note -ve sign since second derivative of psi is estimated*),

grKppsmooth=ListPlot[Kppdata,Joined->True,PlotStyle->Hue[0.75],BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)((Δt)^(K'' (1)))"}],
aaKpp=Fit[Take[Kppdata,{droppointslow,droppointshigh},All],{1,x1},x1];
{interKpp,hqnegKpp}={aaKpp[[1]],aaKpp[[2]][[1]]},

Print["alpha from graph slope = ",hqnegKpp/hqnegC1],
grKpptheory=Plot[interKpp+hqnegKpp x1,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],
grKpptheorydata=Show[grKppsmooth,grKpptheory,PlotRange->All],
grallC1Kpp=Show[grC1theorydata,grKpptheorydata],

(* the above make various useful graphs.  The RMS (root mean square) grRMS, can be compared with the first order (q=1) graph, grq1 (together,grall) ; any divergence between them indicates intermittency.  This is more precisely estimated by the graph grC1 whose logarithmic slope is C1.  There is also a choice of smoothed or unsmoothed moments.  Smoothing is often justified for single series since the statistics at 20 lages per order of magnitude are often poor. grKppsmooth gives a curve whose logarithmic slope is an estimate of the second derivative of K(q) at q=1.  This is C1 alpha, hence it can be used to esitmate alpha.*) 

(*The main useful graphs are gralldatatheory which has RMS (top,magenta) and q=1 (bottom, cyan) data with regressions, and grallC1Kpp which has the C1 estimation graph (red)  with regression and the second derivative of K grpaph with regression (purple)*),

datacount=Select[Transpose[{coordint,count}],#[[2]]>10^-10.&], sigmadeltTsq=(10.^moms[[21]]-10.^(2. moms[[11]]))^0.5,
uncert=sigmadeltTsq/(datacount[[All,2]]-1.)^0.5,
numsd=2.,
(*the number of standard deviations to display*),
upper=10.^( moms[[11]])+numsd  uncert,
lower=10.^( moms[[11]])-numsd  uncert,mnsigma=Mean[Select[Transpose[{coordcalib,uncert,datacount[[All,2]]}],#[[3]]≥2.&][[All,2]]],
grsigma=ListPlot[Select[Transpose[{coordcalib,uncert,datacount[[All,2]]}],#[[3]]≥2.&][[All,Range[1,2]]],Joined->True,BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","S(Δt)/(n(Δt)^(1/2))"}],
grsigmanorm=ListPlot[Select[Transpose[{coordcalib,uncert/mnsigma,datacount[[All,2]]}],#[[3]]≥2.&][[All,Range[1,2]]],Joined->True,PlotStyle->Black,BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","σ(Δt)"}],
(*this is the sigma correctly calcluated on a linear scale*),sellupper=Select[Transpose[{coordcalib,Log[10,upper],datacount[[All,2]]}],#[[3]]≥2.&][[All,Range[1,2]]],selluppersmooth=Transpose[{sellupper[[All,1]],smooth1D[sellupper[[All,2]]]}],selllower=Select[Transpose[{coordcalib,Log[10,lower],datacount[[All,2]]}],#[[3]]≥2.&&#[[2]]>-5.&][[All,Range[1,2]]],
selllowersmooth=Transpose[{selllower[[All,1]],smooth1D[selllower[[All,2]]]}],(*this is the smoothed version*),grq1uncertaintysmooth=ListPlot[{selluppersmooth,Transpose[{coordcalib,Log[10,smooth1D[10.^( moms[[11]])]]}],selllowersmooth},PlotStyle->{{Red,Dashing[0.05]},Black,{Red,Dashing[0.05]}},Joined->True,BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","(Subscript[Log, 10])(Subscript[S, 1])(Δt)"}],

grq1uncertainty=ListPlot[{testt1=Select[Transpose[{coordcalib,Log[10,upper],datacount[[All,2]]}],#[[3]]≥2.&][[All,Range[1,2]]];Transpose[{testt1[[All,1]],smooth1D[testt1[[All,2]]]}],Transpose[{coordcalib,Log[10,smooth1D[10.^( moms[[11]])]]}],testt2=Select[Transpose[{coordcalib,Log[10,lower],datacount[[All,2]]}],#[[3]]≥2.&][[All,Range[1,2]]];Transpose[{testt2[[All,1]],smooth1D[testt2[[All,2]]]}]},PlotStyle->{{Red,Dashing[0.05]},Black,{Red,Dashing[0.05]}},Joined->True,BaseStyle->{FontFamily->"Times",FontSize->32},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","(Subscript[Log, 10])(Subscript[S, 1])(Δt)"}],

(*this  uses the second moment to estimate the standard deviation of the first moment.  numsd is the number of standard deviations that are used for the uncertainty bounds). Note we  divided by (N-1)^0.5 the principale curve is for the mean (q=1) the RMS is used for the unceratinty in the q=1 moment (combined with information from the number of fluctuations avialble in the analysis)*),

grschmitt=ListPlot[Log[10,Table[{psi[[q]][[1]],psi[[q]][[1]] psiprime0-psi[[q]][[2]]},{q,1,qmax}]],PlotStyle->Hue[0.9`]],alph=Fit[Log[10,Table[{psi[[q]][[1]],psi[[q]][[1]] psiprime0-psi[[q]][[2]]},{q,1,qmax}]],{1,xx},xx],Print["Schmitt α=  ",alph[[2]][[1]]],grschmitttheory=Plot[alph/.xx->xxx,{xxx,-1.`,0.5`}],grschmitttall=Show[grschmitttheory,grschmitt],
(*this uese the second derivative of psi(q) at q=0 to estimate alpha.  Only valid for alpha>1*)
};{Show[gralldatatheory],Show[grallC1Kpp]}]

The above does a Haar analysis without interpolating the data.   The uncertainties are estimated.  See grq1uncertainty.
NOTE THAT small values get truncated:  multiply by a large number first!  This also calculates the uncertainties using the number of data points, see: grq1uncertainty

ConfidenceCorNoInter[rrcor_,numint_,sdconfidence_]:=Module[{},{dimm=Dimensions[numint][[1]],rrrcor=Table[Max[Min[rrcor[[jkk]],1-10^-3.],-1+10^-3.],{jkk,1,dimm}],zr= Log[(1+rrrcor)/(1-rrrcor)]/2.,
neffectivesamplesize=Table[Max[numint[[jk]]-3,1],{jk,1,dimm}],minneffectivesamplesize=Table[If[neffectivesamplesize[[jk]]==1,1,0],{jk,1,dimm}],coninterval=sdconfidence ((1+0.5 rrrcor^2)/neffectivesamplesize)^0.5,zu=zr+coninterval,zl=zr-coninterval,ru=(E^(2 zu)-1)/(E^(2 zu)+1),rl=(E^(2 zl)-1)/(E^(2 zl)+1)};{minneffectivesamplesize,rl,rrrcor,ru}]

This is for the correlation with HaarNoInterpolate.  minneffectivesamplesize = 1 when insufficient data, 0 otherwise ...

ConfidenceCor[rrcor_,numint_,sdconfidence_]:=Module[{},{dimm=Dimensions[numint][[1]],rrrcor=Table[Max[Min[rrcor[[jkk]],1-10^-3.],-1+10^-3.],{jkk,1,dimm}],zr=Log[(1+rrrcor)/(1-rrrcor)]/2.,
neffectivesamplesize=Table[Max[If[numint[[jk]]<5,2 numint[[jk]],numint[[jk]]]-3,1],{jk,1,dimm}],coninterval=sdconfidence ((1+0.5 rrrcor^2)/neffectivesamplesize)^0.5,zu=zr+coninterval,zl=zr-coninterval,ru=(E^(2 zu)-1)/(E^(2 zu)+1),rl=(E^(2 zl)-1)/(E^(2 zl)+1),minneffectivesamplesize=Table[If[neffectivesamplesize[[jk]]==1,1,0],{jk,1,dimm}]};{minneffectivesamplesize,rl,rrrcor,ru}]

this gives the confidence interval (width sdconfidence number of standard deviations) around the measured correlation coefficeint using the theory in Bonnet and Wright 2000.  numint is from the numinterval value in StructureMomentDisAllq2D21newHaarErrors.  this takes into account the possiblity of more than one realization (= num).  There is also a factor of 2 to take into account the series reversal .  The minneffectivesamplesize parameter is 1 when the smaple size is too small, otherwise, =0.

HaarCorrelationNoInterpolate[comtime_,data1_,data2_,calibcoord_]:=Module[{},{sdconf=1.,HaarNoInterpolate[comtime,data1,calibcoord,1,1],moms1=moms[[21]],moms1q1=moms[[11]],dimmen=Dimensions[moms1][[1]],grRMSsmooth1=grRMSsmooth,HaarNoInterpolate[comtime,data2,calibcoord,1,1],moms2=moms[[21]],moms2q1=moms[[11]],grRMSsmooth2=grRMSsmooth,HaarNoInterpolate[comtime,data1-data2,calibcoord,1,1],momsdifference=moms[[21]],correlation=0.5 (10.^moms1+10.^moms2-10.^momsdifference),correlationcoefLog=Transpose[{coordcalib,Log[10,Abs[correlation]/(10.^((moms1+moms2)/2.))]}],correlationcoefLinear=Transpose[{coordcalib,correlation/(10.^((moms1+moms2)/2.))}],datacount=Select[Transpose[{coordint,count}],#[[2]]>10^-10.&],confintervals=ConfidenceCorNoInter[correlationcoefLinear[[All,2]],datacount[[All,2]],sdconf],grconlow=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[2]]}],Joined->True,PlotStyle->{Red,Dashing[0.05]}],grconhigh=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[4]]}],Joined->True,PlotStyle->{Red,Dashing[0.05]}],grconmean=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[3]]}],Joined->True,PlotStyle->Black],
grscatterq1=Show[ListPlot[10^Transpose[{moms1q1,moms2q1}],Joined->True,BaseStyle->{32,FontFamily->"Arial"},AxesLabel->{"(Subscript[q1, 1])","(Subscript[q1, 2])"}],Graphics[Table[{Hue[jj/dimmen],PointSize[0.05 (jj/dimmen)^0.5],Point[10.^Transpose[{moms1q1,moms2q1}][[jj]]]},{jj,1,dimmen}]]],
correlationq1=Correlation[10.^Transpose[{moms1q1,moms2q1}]][[2,1]],
correlationRMS=Correlation[10.^Transpose[{moms1,moms2}/2.]][[2,1]],
Print[" The scale by scale correlation of the q=1 fluctuations is",correlationq1],
Print[" The scale by scale correlation of the RMS fluctuations is",correlationRMS],grscatterRMS=Show[ListPlot[10^Transpose[{moms1,moms2}/2.],Joined->True,BaseStyle->{32,FontFamily->"Arial"},AxesLabel->{"(Subscript[RMS, 1])","(Subscript[RMS, 2])"}],Graphics[Table[{Hue[jj/dimmen],PointSize[0.05 (jj/dimmen)^0.5],Point[10.^Transpose[{moms1,moms2}/2.][[jj]]]},{jj,1,dimmen}]]],(*the RMS correlation is for the points at all scales in scatter plot, in the latter the area of points is proportional to the log of the time scale*),grsuff=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[1]]}],Joined->True,PlotStyle->{Brown,Thickness[0.008]}],grconfall=Show[grsuff,grconmean,grconhigh,grconlow,PlotRange->All],corrlowpickzero=Pick[Transpose[{correlationcoefLinear[[All,1]],confintervals[[2]]}],confintervals[[1]],0],grconlowzero=ListPlot[corrlowpickzero,Joined->True,PlotStyle->{Red,Dashing[0.05]}],corrhighpickzero=Pick[Transpose[{correlationcoefLinear[[All,1]],confintervals[[4]]}],confintervals[[1]],0],grconhighzero=ListPlot[corrhighpickzero,Joined->True,PlotStyle->{Red,Dashing[0.05]}],
corrlinearpickzero=Pick[Transpose[{correlationcoefLinear[[All,1]],confintervals[[3]]}],confintervals[[1]],0],grconmeanzero=ListPlot[corrlinearpickzero,Joined->True,PlotStyle->Black],grconfall=Show[grsuff,grconmean,grconhigh,grconlow,PlotRange->All],grconfallzero=Show[grconmeanzero,grconhighzero,grconlowzero,PlotRange->All,BaseStyle->{32,FontFamily->"Arial"},AxesLabel->{"(Subscript[Log, 10])Δt","ρ"}],(*This shows the confidence limits with sdconf number of standard devations, high, low (red), mean (black) and insufficient (brown) when the latter is 1 (otherwise it is 0), the graphs with "zero" at the end are those that only show the correlations that have enough points (more than 3 data points)*),grcorr=ListPlot[correlationcoefLog,Joined->True,PlotStyle->Blue,BaseStyle->{24,FontFamily->"Arial"}],correlationlinear=Transpose[{Transpose[correlationcoefLinear][[1]],smooth1D[Transpose[correlationcoefLinear][[2]]]}],grcorrLinear=ListPlot[correlationlinear,PlotStyle->Red,Joined->True,PlotRange->All,BaseStyle->{24,FontFamily->"Arial"}],grcorrall=Show[grRMSsmooth1/.Hue[_,_,_]->Black,grRMSsmooth2/.Hue[_,_,_]->Brown,grcorr,PlotRange->All],grcorralllin=Show[grRMSsmooth1/.Hue[_,_,_]->Black,grRMSsmooth2/.Hue[_,_,_]->Brown,grcorrLinear,AxesOrigin->{Log[10,calibcoord],0.},PlotRange->All],Print[" black is the first data set fluctuations, brown is the second dat a set (both on log-log scale)  red is the correlation on a log-linear scale"]};Show[grcorralllin]]

this determines the correlation and correlation coefficient between the fluctuations of the series 1 and 2 but without interpolation
Note: “correlation” is no longer defined with a smoothing....

smooth1D[field_]:=Module[{},{{dim}=Dimensions[field],smoothfield=field/2.+(RotateRight[field,1]+RotateLeft[field,1])/4.,smoothfield[[1]]=(field[[1]]+field[[2]])/2.,smoothfield[[dim]]=(field[[dim]]+field[[dim-1]])/2.};smoothfield]

1-D smoother: 1-2-1 weighting

StructureMomentDisAllq2D21newHaarErrorsPrintSwitch[xv_,calibcoord_,droppointslowpoints_,droppointshighpoints_,pswitch_]:=Module[{xx,dim},{{num,dim}=Dimensions[xv],(*if pswitch=1, then print.  Otherwise no print*),
(*xv is a two dimensional array.  num is the nubmer of independent series, dim is the legth of each.  The parameters droppointslow and droppointshighpoints are the number of points at low and high lags respectively that should be dropped before esitmating exponents.  If there is a single series, in the argument, put curly brackets around it, e.g. {xv}.  Alternatively, we could think of this as analysing a 2-D field (num X dim).  In this case we could also analyse the transpose of xv (so as to study the anisotropy. calibcoord is for calibration of the time/space axis it is the resolution in some basic units for example if the data has monthly resolution, then using calibcoord=1/12 give a result in units of years.  The range of scales over which fluctuations are calculated is 2 units to dim i.e. a ratio of dim/2 (we use the convention that the largest lag is equal to the length of the series.  fac=1 is the default.  If fac=2, then periodic data are assumed oand the largest lag is half the full length.)*),

(*calib=2.*),

(*This is the "Haar calibration" i.e. the value that approxiamtely maps the Haar flucutaitons onto difference and tendency fluctuations, the default value is 2.*),

qmax=21,delq=0.1,qeps=0.0001`,Qtab=Table[q,{q,qeps,delq qmax+qeps-delq,delq}],

(*with these parameters we calucate the moments every delq=0.1 from 0 to q=2 (qmax values in all).  qeps is added to avoid exponenting to a zero power - a problem if the data contains zero values.  calib is an overall multiplier used only for the graphics output grRMS, grq1 that takes into account the calibration with respect to difference and tendency fluctuations*)

points2=10^6,points1=IntegerPart[20. Log[10,(points2-0.001)/fac]],
δδt30={2,4,6,8,10,12,14,16,18,20,22},maxpoint=2 IntegerPart[dim/(2 fac)],δδt2=Flatten[Append[δδt30,Table[2  IntegerPart[((10.^(jj/20.)))/2],{jj,28,points1}]]],δδt=Select[δδt2,#<(maxpoint+0.001)&],points=Dimensions[δδt][[1]],δδt[[points]]=If[δδt[[points]]==dim,δδt[[points]]=δδt[[points]]-2,δδt[[points]]],coord= Log[10.,δδt],

If[pswitch==1,Print[" The total number of points is ",points],Continue],

(*this establishes the (integer) lags that we will use; they must be even and are spaced at roughly 20 per order of magnitude (although this is not possible at the smallest lags, hence the complications).  The variable coord has the log10 values of the lags and is quite useful, δδt has the actual coordinates without logs, points is the total number of lags the possibility that the lag is exactly equal to the size of the system (and hence imply a subscript zero) is avoided by reducing the lag by 2 units in those rare special cases*)

moms=Table[0,{qmax},{points}],(*initializing the table with the moments: moms*),
numinterval= num  Table[IntegerPart[(dim-1)/δδt[[jjj]]],{jjj,1,points}],
(* this is the number of data points in a lag of length δδt[[jjj]], used for confidence level statistics, note that the reversal of the series is not taken into account since the result is not very statistically different!*)

Do[{R=Running[xv[[i]]-Mean[xv[[i]]]],RR=Reverse[R],


(*the reveresed series to get better large lag statistics*)

moms=moms+Transpose[Table[δt=δδt[[jj]];Rδt=RotateRight[R,δt];Rδt12=RotateRight[R,δt/2];δRxsign=2 (R+Rδt-2 Rδt12)/δt;
δRx= Abs[δRxsign];RRδt=RotateRight[RR,δt];RRδt12=RotateRight[RR,δt/2];
δRRxsign=2  (RR+RRδt-2 RRδt12)/δt;δRRx=Abs[δRRxsign];Mean[(((#1)^Qtab)&)/@Flatten[{Take[δRx,{δt+1,dim,δt}],Take[δRRx,{δt+1,dim,δt}]}]],{jj,1,points}]]},{i,1,num}],moms=Log[10,(moms/num)],

(*moms has the logs of all the moments and all the lags *)

droppointslow=Max[1,1+droppointslowpoints],
droppointshigh=Max[droppointslow+2,points-droppointshighpoints],


(*droppointslow=1,droppointshigh=points, these indicate the range of data used in the fits  All the code below is purely for fitting, it is not strictly needed.  In this default mode all points between the smallest lag 1 and largest available lag (=points) are taken for regressions.  Indeed, the fitted parameters are only valid if the range of fitting is pretty linear on a log-log plot.  Often it is best to make the fits after viewing the plots*),

coordcalib=coord+Log[10,calibcoord],
If[pswitch==1,Print[" Fits are over the range   ", 10^coordcalib[[droppointslow]],"  to  ",10^coordcalib[[droppointshigh]], "  in the specified units"],Continue],
(*the calibrated coordinate*),

(*{inter,hqneg}=Transpose[Table[aa=Fit[Take[Transpose[{coordcalib,moms[[q]]}],{droppointslow,droppointshigh},{1,2}],{1,x1},x1];{aa[[1]],aa[[2]][[1]]},{q,1,qmax}]],*)

{inter,hqneg}=Transpose[Table[lm=LinearModelFit[Take[Transpose[{coordcalib,moms[[q]]}],{droppointslow,droppointshigh},{1,2}],x,x];aa=Normal[lm];pmm=Normal[lm["ParameterConfidenceIntervals"]];sdd[q]=(pmm[[2,2]]-pmm[[2,1]])/4.;{aa[[1]],aa[[2]][[1]]},{q,1,qmax}]],sdHpsi=Mean[Table[sdd[q],{q,9,13}]],

Do[grref[i]=Plot[Log[10,calib]+inter[[i]]+hqneg[[i]] x1,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],{i,1,qmax}],Do[gr[i]=ListPlot[Transpose[{coordcalib,Log[10,calib]+moms[[i]]}],Joined->True,PlotRange->All,PlotStyle->{Hue[(i/qmax)],Thickness[0.005]},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)S(Δt)"},BaseStyle->{24,FontFamily->"Times"}],{i,1,qmax}],Do[grpts[i]=ListPlot[Transpose[{coordcalib,Log[10,calib]+moms[[i]]}],PlotRange->All,PlotStyle->{PointSize[0.008],Hue[(i/qmax)]}],{i,1,qmax}],grout=Show[Table[gr[i],{i,1,qmax}],Table[grpts[i],{i,1,qmax}],Table[grref[i],{i,1,qmax}]],psi=Table[{Qtab[[q]],hqneg[[q]]},{q,1,qmax}],grpsi=ListPlot[psi,Joined->True,PlotStyle->Hue[0.75`],AxesLabel->{"q","ξ(q)"},BaseStyle->{24,FontFamily->"Times"}],

ft0=Fit[Flatten[Append[{{{0,0.}}},Take[psi,{1,2},{1,2}]],1],{1,xx,(xx^2)},xx],psiprime0=ft0[[1]]+ft0[[2]][[1]],ft=Fit[Take[psi,{9,13},{1,2}],{1,xx,(xx^2)},xx],C1psi=ft[[1]]-ft[[3]][[1]],Hpsi=ft/.xx->1,If[pswitch==1,Print["H=   ",Hpsi," ±",sdHpsi],Continue],If[pswitch==1,Print["the RMS exponent= ξ(2)/2= ",psi[[21,2]]/2.,"±",sdd[21]/2.,"  the difference with H is K(2)/2 = ",psi[[11,2]]-psi[[21,2]]/2.],Continue],If[pswitch==1,Print["derivative C1=  ",C1psi],Continue],Kppq1=-2 ft[[3]][[1]],αpsi=Kppq1/C1psi,If[pswitch==1,Print["α from psipp(1) =  ",αpsi],Continue],

(*The above are for regressions, parameter estimates over the range droppointslow,droppointshigh*),

C1data=Transpose[{coordcalib,Log[10,smooth1D[smooth1D[10.^(moms[[11]]-(5.(moms[[12]]-moms[[10]])))]]]}],

grC1=ListPlot[C1data,Joined->True,PlotStyle->Hue[0.99],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)((Δt)^(Subscript[C, 1]))"}],
grinfo=grC1,

lm=LinearModelFit[Take[C1data,{droppointslow,droppointshigh},All],x1,x1],aa=Normal[lm];pmm=Normal[lm["ParameterConfidenceIntervals"]];sdddC1=(pmm[[2,2]]-pmm[[2,1]])/4.;
{interC1,hqnegC1}={aa[[1]],aa[[2]][[1]]},grC1theory=Plot[interC1+hqnegC1 x1,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],
If[pswitch==1,Print["C1 graph slope method =",hqnegC1," ±",sdddC1],Continue],
grC1theorydata=Show[grC1,grC1theory,PlotRange->All],

grRMS=ListPlot[Transpose[{coordcalib,Log[10,calib]+moms[[21]]/2.}],Joined->True,PlotRange->All,PlotStyle->Hue[0.9],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 2])(Δt)^(1/2)"}],RMSdata=Transpose[{coordcalib,Log[10,calib]+Log[10,smooth1D[10^moms[[21]]]]/2.}],grRMSsmooth=ListPlot[RMSdata,Joined->True,PlotRange->All,PlotStyle->Hue[0.9],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 2])(Δt)^(1/2)"}],grcorsmooth=ListPlot[Transpose[{coordcalib,Log[10,calib]+Log[10,smooth1D[10^moms[[21]]]]}],Joined->True,PlotStyle->Hue[0.1],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 2])(Δt)"}],
grboxsmooth=ListPlot[Transpose[{coordcalib,Log[10,calib]+Log[10,smooth1D[10^moms[[1]]]]}],Joined->True,PlotStyle->Hue[0.4],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 0])(Δt)"}],
q1data=Transpose[{coordcalib,Log[10,calib]+Log[10,smooth1D[10^moms[[11]]]]}],
        grq1smooth=ListPlot[q1data,PlotRange->All,Joined->True,PlotStyle->Hue[0.6],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 1])(Δt)"}],
grq1=ListPlot[Transpose[{coordcalib,Log[10,calib]+moms[[11]]}],Joined->True,PlotRange->All,PlotStyle->Hue[0.6],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 1])(Δt)"}],
grtheoryq1=Plot[Log[10,calib]+inter[[11]]+hqneg[[11]] x1,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],
grtheorydataq1=Show[grq1smooth,grtheoryq1],
grtheoryRMS=Plot[Log[10,calib]+inter[[21]]/2+hqneg[[21]] x1/2,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],
grtheorydataRMS=Show[grRMSsmooth,grtheoryRMS,PlotRange->All],

grallsmooth=Show[grq1smooth,grRMSsmooth,PlotRange->All],grall=Show[grq1,grRMS,PlotRange->All],gralldatatheory=Show[grtheorydataq1,grtheorydataRMS,PlotRange->All],

Kppdata=Transpose[{coordcalib,-25.Log[10,smooth1D[10.^(moms[[13]]+moms[[9]]-2 moms[[11]])]]}],
(*note -ve sign since second derivative of psi is estimated*),

grKppsmooth=ListPlot[Kppdata,Joined->True,PlotStyle->Hue[0.75],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)((Δt)^(K'' (1)))"}],
aaKpp=Fit[Take[Kppdata,{droppointslow,droppointshigh},All],{1,x1},x1];
{interKpp,hqnegKpp}={aaKpp[[1]],aaKpp[[2]][[1]]},
If[pswitch==1,Print["alpha from graph slope = ",hqnegKpp/hqnegC1],Continue],
grKpptheory=Plot[interKpp+hqnegKpp x1,{x1,coordcalib[[droppointshigh]],coordcalib[[droppointslow]]}],
grKpptheorydata=Show[grKppsmooth,grKpptheory,PlotRange->All],
grallC1Kpp=Show[grC1theorydata,grKpptheorydata],


(* the above make various useful graphs.  The RMS (root mean square) grRMS, can be compared with the first order (q=1) graph, grq1 (together,grall) ; any divergence between them indicates intermittency.  This is more precisely estimated by the graph grC1 whose logarithmic slope is C1.  There is also a choice of smoothed or unsmoothed moments.  Smoothing is often justified for single series since the statistics at 20 lages per order of magnitude are often poor. grKppsmooth gives a curve whose logarithmic slope is an estimate of the second derivative of K(q) at q=1.  This is C1 alpha, hence it can be used to esitmate alpha.*) 

(*The main useful graphs are gralldatatheory which has RMS (top,magenta) and q=1 (bottom, cyan) data with regressions, and grallC1Kpp which has the C1 estimation graph (red)  with regression and the second derivative of K grpaph with regression (purple)*),

Ratioq1toRMSsmoothdata=Transpose[{coordcalib,Log[10,smooth1D[smooth1D[10^(moms[[11]]-moms[[21]]/2.)]]]}],
grRatioq1toRMSsmooth=ListPlot[Ratioq1toRMSsmoothdata,Joined->True,PlotRange->All,PlotStyle->Hue[0.9],BaseStyle->{FontFamily->"Times",FontSize->24},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)Δt","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(Subscript[S, 1])(Δt)/(Subscript[S, 2])(Δt)^(1/2)"}],
gralphaall=If[pswitch==1,alphaestimate[Ratioq1toRMSsmoothdata,C1data],Continue],

grschmitt=ListPlot[Log[10,Table[{psi[[q]][[1]],psi[[q]][[1]] psiprime0-psi[[q]][[2]]},{q,1,qmax}]],PlotStyle->Hue[0.9`]],alph=Fit[Log[10,Table[{psi[[q]][[1]],psi[[q]][[1]] psiprime0-psi[[q]][[2]]},{q,1,qmax}]],{1,xx},xx],If[pswitch==1,Print["Schmitt α=  ",alph[[2]][[1]]],Continue],grschmitttheory=Plot[alph/.xx->xxx,{xxx,-1.`,0.5`}],grschmitttall=Show[grschmitttheory,grschmitt],
(*this uese the second derivative of psi(q) at q=0 to estimate alpha.  Only valid for alpha>1*)};{If[pswitch==1,{Show[gralldatatheory,ImageSize->600],Show[grallC1Kpp,ImageSize->600]},Continue]}]

This has a switch parameter that allows the printing to be shut down : if pswitch = 1, then print, otherwise, no print.  Otherwise this is the same as StructureMomentDisAllq2D21newHaarErrors.

HaarCorrelationPrintSwitch[data1_,data2_,calibcoord_,pswitch_]:=Module[{},{sdconf=1.,StructureMomentDisAllq2D21newHaarErrorsPrintSwitch[{data1},calibcoord,1,1,pswitch],moms1=moms[[21]],grRMSsmooth1=grRMSsmooth,StructureMomentDisAllq2D21newHaarErrorsPrintSwitch[{data2},calibcoord,1,1,pswitch],moms2=moms[[21]],grRMSsmooth2=grRMSsmooth,StructureMomentDisAllq2D21newHaarErrorsPrintSwitch[{data1-data2},calibcoord,1,1,pswitch],momsdifference=moms[[21]],correlation=0.5 (10.^moms1+10.^moms2-10.^momsdifference),correlationcoefLog=Transpose[{coordcalib,Log[10,Abs[correlation]/(10.^((moms1+moms2)/2.))]}],correlationcoefLinear=Transpose[{coordcalib,correlation/(10.^((moms1+moms2)/2.))}],
confintervals=ConfidenceCor[correlationcoefLinear[[All,2]],numinterval,sdconf],

grconlow=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[2]]}],Joined->True,PlotStyle->{Red,Dashing[0.05]}],grconhigh=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[4]]}],Joined->True,PlotStyle->{Red,Dashing[0.05]}],grconmean=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[3]]}],Joined->True,PlotStyle->Black],grsuff=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[1]]}],Joined->True,PlotStyle->{Brown,Thickness[0.008]}],grconfall=Show[grsuff,grconmean,grconhigh,grconlow,PlotRange->All],(*This shows the confidence limits with sdconf number of standard devations, high, low (red), mean (black) and insufficient (brown) when the latter is 1 (otherwise it is 0)*),
grnoconf=ListPlot[Transpose[{correlationcoefLinear[[All,1]],confintervals[[1]]}],Joined->True,PlotStyle->{Brown,Thickness[0.008]}],grconfall=Show[grconmean,grconhigh,grconlow,PlotRange->All],(*the same but without the brown indicating insufficient data*),

grcorr=ListPlot[correlationcoefLog,Joined->True,PlotStyle->Blue,BaseStyle->{24,FontFamily->"Arial"}],correlationlinear=Transpose[{Transpose[correlationcoefLinear][[1]],smooth1D[Transpose[correlationcoefLinear][[2]]]}],grcorrLinear=ListPlot[correlationlinear,PlotStyle->Red,Joined->True,PlotRange->All,BaseStyle->{24,FontFamily->"Arial"}],grcorrall=Show[grRMSsmooth1/.Hue[_,_,_]->Black,grRMSsmooth2/.Hue[_,_,_]->Brown,grcorr,PlotRange->All]};grcorralllin=Show[grRMSsmooth1/.Hue[_,_,_]->Black,grRMSsmooth2/.Hue[_,_,_]->Brown,grcorrLinear,AxesOrigin->{Log[10,calibcoord],0.},PlotRange->All]]

This has a switch parameter that allows the printing to be shut down : if pswitch = 1, then print, otherwise, no print.

Below is an example of the use of some of  the above software.  The data and results are described in:

Spiridonov, A., L. Balakauskas, S. Lovejoy, 2022: Longitudinal expansion fitness of brachiopod genera controlled by the Wilson cycle, Global and Planetary Change, https://doi.org/10.1016/j.gloplacha.2022.103926.

http://www.physics.mcgill.ca/~gang/eprints/eprintLovejoy/neweprint/Brachiopods.proofs.1-s2.0-S092181812200193X-main.pdf

ageSong={0.`,1.`,2.`,3.`,4.`,5.`,6.`,7.`,8.`,9.`,10.`,11.`,12.`,13.`,14.`,15.`,16.`,17.`,18.`,20.`,22.`,23.`,24.`,26.`,27.`,28.`,29.`,30.`,32.`,33.`,34.`,35.`,37.`,39.`,41.`,42.`,43.`,44.`,45.`,46.`,47.`,48.`,49.`,50.`,51.`,52.`,53.`,54.`,55.`,56.`,57.`,58.`,59.`,60.`,61.`,62.`,63.`,66.`,67.`,68.`,69.`,70.`,71.`,72.`,73.`,74.`,75.`,76.`,77.`,78.`,79.`,80.`,82.`,84.`,86.`,87.`,88.`,89.`,90.`,91.`,92.`,93.`,94.`,95.`,96.`,97.`,98.`,99.`,100.`,101.`,102.`,103.`,104.`,105.`,106.`,107.`,108.`,109.`,112.`,116.`,117.`,118.`,119.`,120.`,121.`,122.`,123.`,124.`,125.`,126.`,127.`,128.`,129.`,130.`,131.`,132.`,133.`,134.`,135.`,136.`,137.`,138.`,139.`,140.`,141.`,142.`,145.`,150.`,153.`,154.`,156.`,157.`,158.`,159.`,160.`,161.`,162.`,164.`,165.`,166.`,167.`,168.`,169.`,170.`,171.`,172.`,173.`,174.`,175.`,176.`,177.`,178.`,179.`,180.`,181.`,182.`,183.`,184.`,185.`,186.`,187.`,188.`,189.`,190.`,191.`,192.`,193.`,195.`,197.`,198.`,199.`,200.`,201.`,202.`,203.`,204.`,205.`,206.`,207.`,209.`,210.`,211.`,212.`,214.`,217.`,222.`,226.`,230.`,231.`,232.`,234.`,235.`,236.`,237.`,238.`,240.`,242.`,244.`,245.`,246.`,247.`,248.`,249.`,250.`,251.`,252.`,253.`,254.`,255.`,256.`,257.`,258.`,259.`,260.`,261.`,262.`,264.`,265.`,266.`,268.`,270.`,271.`,272.`,274.`,275.`,276.`,277.`,279.`,282.`,285.`,288.`,290.`,291.`,292.`,293.`,294.`,295.`,296.`,297.`,298.`,299.`,300.`,301.`,302.`,303.`,304.`,305.`,306.`,307.`,308.`,310.`,311.`,312.`,313.`,314.`,315.`,316.`,317.`,318.`,319.`,321.`,322.`,323.`,324.`,325.`,326.`,327.`,328.`,329.`,330.`,331.`,332.`,333.`,334.`,335.`,336.`,339.`,345.`,347.`,348.`,349.`,350.`,351.`,352.`,353.`,355.`,356.`,357.`,358.`,359.`,360.`,361.`,362.`,363.`,364.`,365.`,366.`,367.`,368.`,369.`,370.`,371.`,372.`,373.`,374.`,375.`,376.`,379.`,381.`,382.`,383.`,384.`,385.`,386.`,387.`,388.`,389.`,390.`,391.`,392.`,393.`,394.`,395.`,396.`,398.`,400.`,403.`,405.`,406.`,407.`,409.`,410.`,412.`,414.`,415.`,416.`,421.`,422.`,423.`,424.`,425.`,426.`,427.`,428.`,429.`,430.`,431.`,432.`,433.`,434.`,435.`,436.`,437.`,438.`,439.`,441.`,442.`,443.`,444.`,445.`,446.`,447.`,448.`,449.`,450.`,451.`,452.`,453.`,454.`,459.`,460.`,462.`,466.`,468.`,472.`,474.`,477.`,480.`,482.`,484.`,485.`,493.`,494.`,495.`,496.`,497.`,498.`}

Song age in Myrs

TSong

TSong={18.22`,15.38`,16.19`,17.16`,17.47`,17.13`,15.19`,15.67`,16.11`,16.99`,18.94`,17.78`,16.57`,15.25`,16.74`,17.09`,17.64`,16.83`,15.43`,16.64`,19.03`,18.61`,18.49`,14.09`,15.24`,14.47`,15.22`,15.67`,21.8`,21.25`,21.79`,22.17`,16.21`,16.8`,17.82`,18.07`,18.43`,15.74`,21.54`,19.8`,17.41`,21.81`,23.75`,22.98`,20.88`,22.07`,22.38`,20.43`,26.45`,20.06`,18.33`,18.34`,20.06`,17.15`,17.49`,15.81`,15.62`,18.81`,17.09`,17.92`,15.82`,15.96`,15.89`,16.25`,17.3`,16.95`,16.45`,20.26`,19.08`,19.47`,18.06`,16.7`,27.63`,22.07`,29.91`,28.72`,29.79`,29.71`,30.45`,31.97`,31.93`,31.75`,31.66`,30.03`,24.98`,28.37`,27.86`,24.77`,25.44`,24.66`,24.1`,24.01`,23.54`,22.89`,22.08`,20.9`,19.91`,20.6`,19.45`,22.86`,24.28`,24.33`,22.68`,21.66`,20.22`,21.49`,21.42`,20.87`,20.74`,20.84`,20.5`,17.87`,25.24`,23.4`,21.25`,21.3`,22.`,20.69`,18.56`,18.17`,16.29`,16.94`,19.39`,17.8`,18.06`,17.25`,17.85`,19.3`,19.29`,20.44`,24.39`,19.12`,20.51`,20.07`,21.8`,19.4`,16.76`,19.43`,23.52`,22.66`,25.16`,24.88`,25.34`,20.05`,23.02`,25.71`,21.07`,23.05`,26.2`,29.07`,24.46`,26.51`,25.41`,26.58`,24.08`,22.23`,18.88`,18.38`,20.02`,18.85`,21.91`,22.54`,21.14`,21.42`,20.42`,20.44`,21.67`,21.16`,22.88`,16.44`,19.66`,20.81`,15.24`,12.08`,15.61`,18.35`,19.65`,20.04`,20.98`,21.46`,22.81`,23.53`,24.41`,26.29`,24.21`,19.93`,21.44`,23.46`,23.52`,24.17`,22.64`,19.57`,24.72`,22.27`,25.25`,27.97`,26.11`,31.39`,28.53`,22.98`,23.87`,25.98`,29.28`,29.8`,29.94`,23.84`,18.57`,20.21`,20.54`,21.29`,22.86`,23.11`,22.93`,20.56`,23.25`,18.79`,19.3`,18.5`,17.`,20.59`,18.42`,17.3`,17.63`,17.61`,16.54`,14.38`,15.06`,17.41`,14.54`,15.45`,17.42`,15.66`,14.13`,16.4`,14.15`,12.32`,16.91`,15.71`,14.23`,12.85`,17.47`,17.2`,20.`,22.45`,20.51`,18.1`,20.49`,22.13`,21.25`,22.18`,21.73`,20.51`,15.43`,16.29`,13.5`,11.28`,10.94`,11.77`,9.42`,9.37`,10.54`,11.18`,12.37`,12.97`,13.5`,12.33`,15.36`,15.84`,15.13`,13.76`,14.04`,13.`,13.47`,13.04`,13.02`,14.87`,18.54`,20.6`,20.12`,20.33`,21.14`,23.04`,23.92`,23.67`,27.34`,29.36`,28.04`,27.94`,27.17`,27.03`,28.43`,27.57`,31.35`,29.4`,29.85`,30.03`,30.77`,31.74`,33.1`,34.73`,32.95`,32.65`,31.22`,33.52`,31.64`,30.89`,31.11`,24.59`,25.77`,24.35`,25.59`,31.7`,34.44`,25.15`,24.89`,24.31`,28.26`,27.53`,26.08`,25.03`,25.03`,26.26`,24.92`,24.42`,23.73`,24.54`,33.18`,30.55`,27.`,26.7`,28.9`,29.86`,34.42`,35.55`,33.54`,32.46`,31.8`,29.39`,30.03`,27.84`,32.54`,30.69`,28.36`,28.02`,25.73`,26.92`,28.04`,26.02`,30.22`,28.52`,31.28`,30.14`,31.26`,30.6`,29.17`,29.21`,26.87`,26.14`,24.71`,26.14`,25.28`,29.78`,28.25`,29.04`,28.37`,27.79`,27.93`,27.86`,27.13`,31.14`,30.56`,34.71`,32.74`,33.79`,35.55`,38.49`,37.41`,40.61`,39.55`,42.7`,43.88`,46.87`,48.88`,50.19`,47.37`,42.97`,45.82`}

Song temps in oC

Dimensions[TSong]

Dimensions[TSong]

ListPlot[Transpose[{ageSong,TSong}],Joined->True]

ageextbio={517.5`,511.5`,506.75`,502.5`,498.75`,495.5`,491.75`,487.45`,481.55`,473.85`,468.65`,462.85`,455.7`,449.1`,444.3`,442.1`,439.65`,435.95`,431.95`,428.95`,426.5`,424.3`,421.1`,415.`,409.2`,400.45`,390.5`,385.2`,377.45`,365.55`,352.8`,338.8`,327.05`,319.2`,311.1`,305.35`,301.3`,297.2`,292.8`,284.7`,275.8`,270.55`,266.95`,262.5`,257.035`,253.17`,251.685`,249.2`,244.6`,239.5`,232.5`,218.25`,204.9`,200.3`,195.05`,186.75`,178.4`,172.2`,169.3`,167.2`,164.8`,160.4`,154.7`,148.55`,142.4`,136.35`,131.15`,127.2`,119.`,106.75`,97.2`,91.85`,88.05`,84.95`,77.85`,69.05`,63.8`,58.8`,51.9`,44.55`,39.65`,35.95`,31.`,25.565`,19.5`,13.795`,8.4765`,3.9605`}

these are the ages of the biodata base that we used.  The first 4 points are too old (>500Myrs) and were eliminated

TSonginter=Interpolation[Transpose[{ageSong,TSong}]]

extdatabSong[0]=Table[TSonginter[ageextbio[[j]]],{j,5,88}]

this is the Song temperatures interpolated to the biodata base times only the first 500 Myrs were used

extdatab[2]={0.5837130112576242`,0.6097358610287622`,0.7961794722984478`,1.2967095708817937`,1.1501409580221502`,0.17115546764562467`,0.7692660692846008`,0.6584276082222609`,1.0615333360790011`,0.30807093507981503`,0.06565387045962608`,0.5654702674634876`,0.2325636224045869`,0.69150506718677`,0.4019580812639008`,0.09425269631646545`,0.09905549384448015`,0.26594288790738646`,0.1713508403880515`,0.13338301394044666`,0.12072806514572935`,0.07979911020630251`,0.23066159834633612`,0.07864021679036312`,0.0998037016090045`,0.3514105441208212`,0.4320481590406668`,0.41783833206638343`,0.3004229441834304`,0.6462736324763594`,0.20724956502250813`,0.3243854732658298`,0.17853581075204938`,0.1729626383379231`,0.3409411608388917`,0.1466111840465934`,0.03453318300821983`,0.16272053268135792`,-0.008207455689638807`,0.11663776848362792`,0.19174062733659525`,0.14848854292968197`,0.11440382191037805`,0.3524188549019554`,0.09445772432034216`,1.4447871577340559`,0.48655153400359497`,0.2659111463050304`,0.05521728152348486`,0.3392442857656595`,0.3553943942647876`,0.46981032768673187`,0.8351259704093075`,0.15690932120396617`,0.16638457577443255`,0.25811947222794884`,0.1529768429255419`,0.19108601289484803`,0.25674068555964774`,0.0733887587640002`,0.18008261808540357`,0.1085565698726875`,0.22023109651340686`,0.5644022724985689`,0.48111962110982326`,0.08647620288802271`,0.1626928200079793`,0.13505540519606163`,0.34704747727791885`,0.30128227173297634`,0.4500994545342102`,0.28341899221213196`,-0.012648422920668272`,0.08527331955047664`,0.4216985515363849`,0.5137097080295049`,0.10627637181182406`,0.08211946640178583`,0.09273091732599967`,0.013829482055036138`,0.0002887345119478402`,0.2632798820215299`,-0.004660671275221842`,-0.00029348274720398805`,0.04439471049187963`,0.10500228520153045`,0.0720708402754312`,0.2546820297331648`}

extinction rates from the biopaleo data base.

Dimensions[extdatab[2]]

HaarNoInterpolateUncertainty[ageextbio[[Range[5,88]]],extdatabSong[0],1.,3,8]

Song temperature analysis over last 500 Myrs resampled at the biodatabase times.

grq1uncertainty

the dashed lines show them 2 standard deviation level uncertainties (the parameter numsd=2.).

datacount

for each lag, this give the number of fluctuations that were used.

ListPlot[datacount,Joined->True]

number of data points as a function of the log10 of the lag Δt

grpsi

The scaling exponent of the moments every 0.1 up to q = 2.  The fact that it is nearly linear means that C1 ≈ 0., very small intermittency.

HaarCorrelationNoInterpolate[ageextbio[[Range[5,88]]],extdatabSong[0],extdatab[2][[Range[5,88]]],1.]

The analysis of the temperature (black, extinction rate (brown) and correlations (red) over the paleobiodatabase.

grconfallzero

This gives the temeprature- extinctino rate correlations with 1 standard deviation confidence limits (the parameter sdconf=1)