diff --git a/triangle_autocorrelations.m b/triangle_autocorrelations.m
new file mode 100644
index 0000000000000000000000000000000000000000..5786766af6002554b0c31352a5edbe5f3f20491b
--- /dev/null
+++ b/triangle_autocorrelations.m
@@ -0,0 +1,59 @@
+(* ::Package:: *)
+
+(* ::Section:: *)
+(*Data import*)
+
+
+gsraw=Import[NotebookDirectory[]<>"data/graphdata_timeevol.m"];
+(* gsraw=SortBy[gsraw,{#[[1,1]]&,#[[1,2]]&}]; (* Sort by n and then by tau. The {} forces a *stable* sort because otherwise Mathematica sorts also on triangle count and other things. *) *)
+
+
+gdata=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}];
+(* Data format: *)
+(* gdata[[ tau index, n index, run index , datatype index ]] *)
+(* datatype index:
+1: {n,tau}
+2: #triangles time sequence
+3: degree sequence
+4: GCM1 starting triangle counts
+5: GCM2 starting triangle counts
+6: GCM1 time sequence
+7: GCM2 time sequence
+*)
+nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]];
+taulabels=Map["tau = "<>ToString[#]&,gdata[[All,1,1,1,2]]];
+
+
+(* ::Section:: *)
+(*Triangle counts autocorrelation*)
+
+
+doAutoCorrelate[dataset_,maxlag_]:=Module[{avg,data2,result},
+avg=Mean[dataset];
+data2=dataset-avg;
+result=ParallelTable[Mean[data2[[1+lag;;]]*data2[[;;-1-lag]]],{lag,0,maxlag}];
+result=result/result[[1]]
+]
+
+
+doPeriodicAutocorrelate[dataset_,maxlag_]:=Module[{f,result},
+f=Fourier[dataset-Mean[dataset]];
+result=InverseFourier[f*Conjugate[f]];
+result=result[[1;;Min[Floor[Length[result]/2],maxlag]]];
+result=Chop[result];
+result=result/result[[1]]
+]
+
+
+dataX=gdata[[1,1]][[3,2]];
+ListPlot[dataX,Joined->True,PlotRange->All]
+dataX=dataX[[15000;;50000]];
+
+
+cor1=doAutoCorrelate[dataX,2000];
+
+
+cor2=doPeriodicAutocorrelate[dataX,20000];
+
+
+ListPlot[{cor1,cor2},Joined->True,PlotRange->{-0.5,1}]