17 Time Series Analysis

Some Working Data

A random walk

randata = NestList[(# + Random[] - 0.5)&, 0, 1999];

rstream = OpenWrite["rseries.dat"];
Do[
Write[rstream, randata[[i]]],{i,1,2000}];
Close[rstream];

ListPlot[randata]

A noisy oscillatory Series

osccidata = Table[N[Sin[x] + Cos[Pi x] + 0.1*Random[]],
{x, 0, 10, 0.01}];

Needs["Utilities`BinaryFiles`"]

ostream = OpenWriteBinary["oseries.dat"];
Do[
WriteBinary[ostream, osccidata[[j]]],
{j, 1, 1001}];
Close[ostream]

ListPlot[osccidata]

A chaotic Series

NDSolve[{x'[t] == -3 (x[t] - y[t]),
	y'[t] == -x[t] z[t] + 26.5 x[t] - y[t],
	z'[t] ==  x[t] y[t] - z[t],
	x[0] == 0.53755, 
	y[0] ==-0.0071619,
	z[0] == 0.025124},
		 {x, y, z}, {t, 0, 10}, MaxSteps -> 1500]

cdata = Table[Evaluate[x[(j-1)/100] + y[(j-1)/100] -  
3*z[(j-1)/100] /. %],
{j, 1, 1001}];

cdata[[1]]

chaosstream = OpenWrite["cseries.dat"];
Do[
Write[chaosstream, cdata[[j,1]]],
{j, 1, 1001}];
Close[chaosstream]

ListPlot[Flatten[cdata]]

Moving Averages

Needs["Statistics`MovingAverage`"]

madata = ReadList["rseries.dat", Number];

averaged = MovingAverage[madata, 100];

modmadata = Transpose[{Range[2000], madata}];
modaveraged = Transpose[{Range[101, 2000], averaged}];

ListPlot[modaveraged, 
PlotJoined -> True, PlotStyle -> {Thickness[0.002]},
Epilog -> Map[{PointSize[0.001], Point[#]}&, modmadata]]

CompareAverage[data_, n_] := Module[{moddata, 
modaveraged, inter},
inter = MovingAverage[data, n];
moddata = Transpose[{Range[Length[data]], data}];
modaveraged = Transpose[{Range[n + 1, Length[data]], 
inter}];
ListPlot[modaveraged, 
PlotJoined -> True, PlotStyle -> {Thickness[0.002]},
Epilog -> Map[{PointSize[0.001], Point[#]}&, moddata]];
modaveraged
]

mynewave = CompareAverage[madata, 200];

Fourier Analysis

osccidata = ReadList["oseries.dat", Number];

ListPlot[Abs[Fourier[osccidata]]]

tran = Chop[Fourier[osccidata], 0.1];

itran = Chop[InverseFourier[tran], 0.001];

ListPlot[itran, PlotJoined -> True]

Spectral Analysis Package

osccidata = ReadList["oseries.dat", Number];
randata = ReadList["rseries.dat", Number];
chaosdata = ReadList["cseries.dat", Number];

osccidata[[1]]

Length[osccidata]

spectralmatrix[dat_, windowsize_, range_]:=
Table[Sum[ dat[[n+ k - 1]] dat[[n + m - 1]],
{n, 1, range - windowsize + 1}],
{k, 1, windowsize}, {m, 1, windowsize}];

Timing[ mymatrix = spectralmatrix[osccidata, 8, 100];]

TableForm[ mymatrix, TableSpacing -> {1,1}]

specentry[i_,j_,data_,range_,windowsize_,correls_] :=
  specentry[j,i,data,range,windowsize,correls] /; j > i

specentry[i_,j_,data_,range_,windowsize_, correls_] :=
   specentry[i-1,j-1, data, range, windowsize, correls] +
   data[[range-windowsize+i]] data[[range-windowsize+j]]-
   data[[i-1]] data[[j-1]];

specentry[i_,1,data_,range_,windowsize_, correls_] :=
   correls[[i]];

betterspec[data_,windowsize_,range_] := Module[{correls},
    data1 = Take[data,range];
    data2 = Join[Drop[data1,1-windowsize], 
                   Table[0,{windowsize-1}]];
    correls = Table[ Apply[Plus, data1*RotateRight[data2,n-1]],
                    {n,1,windowsize}];
    Table[specentry[i,j,data1,range,windowsize,correls],
          {i,1,windowsize},{j,1,windowsize}]]

This new function runs much faster:

Timing[ mymatrix = betterspec[osccidata,8,100]; ]

TableForm[ mymatrix, TableSpacing -> {1,1}]

{vals, vecs} = Eigensystem[mymatrix]

sorter[{evals_, evecs_}] := 
Module[{bigset},
bigset = MapThread[Join[{#1}, #2]&, {evals, evecs}];
bigset = Reverse[Sort[bigset]];
bigset = Transpose[Map[{First[#], Rest[#]}&, bigset]]
]

{svals, svecs} = sorter[{vals, vecs}]

displayEval[u_List] := 
ListPlot[Log[10, Reverse[Sort[u]]], 
AxesOrigin -> {0,0}, PlotStyle -> PointSize[0.02]]

displayEval[svals]

displayEvec[svecs[[1]]]

displayEvec[svecs[[2]]]

dynamAve[vecs_, data_, range_, n_] := 
Table[
	Sum[
		vecs[[k, 1]]*Sum[vecs[[k, u]] data[[i + u - 1]], 
		{u, 1, Length[vecs]}],
	{k, 1, n}], 
{i, 1, range}]

projecttwo = ListPlot[dynamAve[svecs, osccidata, 900, 2],
PlotJoined -> True]

origplot = ListPlot[Take[osccidata, 900]]

Needs["Spectral`"]

$ContextPath

?Spectral`*

?Spectralmatrix

newtest = Spectralmatrix[randata, 8, 1000];

Spectralmatrix[randata, 6, 5]

Spectralmatrix[randata, 6, 3000]

esystem = Sorter[Eigensystem[newtest]];

DisplayEval[esystem[[1]]]

DisplayEvec[esystem[[2, 1]]]

BeginPackage["Spectral`"]

Spectral::usage = "Spectral.m is a package for time series analysis."

Spectralmatrix::usage = "Spectralmatrix[data, windowsize,
range] computes XT.X for X the matrix of windowed data"

Sorter::usage = "Sorter[Eigensystem] sorts an 
Eigensystem into decreasing order of eigenvalues"

DisplayEval::usage = "DisplayEval[u_List] does a standard log
plot of eigenvalues in decresing order"

DisplayEvec::usage = "DisplayEval[u_List] does a standard plot to
show components of an eigenvector"

DynamAve::usage = "DynamAve[vecs, data, range, n] computes the 
projection of the data onto the eigenvectors corresponding to
the n largest eigenvectors."

Begin["`Private`"]

(* Error messages for the exported objects *)

Spectral::badrangeone = "Range is less than Window Size!";
Spectral::badrangetwo = "Range is bigger than Length of data!";

(* Function Definitions *)

specentry[i_,j_,data_,range_,windowsize_,correls_] :=
  specentry[j,i,data,range,windowsize,correls] /; j > i

specentry[i_,j_,data_,range_,windowsize_, correls_] :=
   specentry[i-1,j-1, data, range, windowsize, correls] +
   data[[range-windowsize+i]] data[[range-windowsize+j]]-
   data[[i-1]] data[[j-1]];

specentry[i_,1,data_,range_,windowsize_, correls_] :=
   correls[[i]];

Spectralmatrix[data_,windowsize_,range_] := Module[{correls},
    data1 = Take[data,range];
    data2 = Join[Drop[data1,1-windowsize], 
                   Table[0,{windowsize-1}]];
    correls = Table[ Apply[Plus, data1*RotateRight[data2,n-1]],
                    {n,1,windowsize}];
    Table[specentry[i,j,data1,range,windowsize,correls],
          {i,1,windowsize},{j,1,windowsize}]] /; Length[dat] >= range && range >= windowsize;
Spectralmatrix[dat_, windowsize_, range_]:=  
Message[Spectral::badrangeone, range] /; range < windowsize;
Spectralmatrix[dat_, windowsize_, range_]:=  
Message[Spectral::badrangetwo, range] /; Length[dat] < range

Sorter[{evals_, evecs_}] := 
Module[{bigset},
bigset = MapThread[Join[{#1}, #2]&, {evals, evecs}];
bigset = Reverse[Sort[bigset]];
bigset = Transpose[Map[{First[#], Rest[#]}&, bigset]] ]

DisplayEval[u_List] := 
ListPlot[Log[10, Reverse[Sort[u]]], 
AxesOrigin -> {0,0}, PlotStyle -> PointSize[0.02]]

DisplayEvec[u_List] := 
ListPlot[u, PlotRange -> {-1, 1}, 
AxesOrigin -> {0, 0}, PlotStyle -> PointSize[0.02]]

DynamAve[vecs_, data_, range_, n_] := 
Table[
	Sum[
		vecs[[k, 1]]*Sum[vecs[[k, u]] data[[i + u - 1]], {u, 1, Length[vecs]}],
		{k, 1, n}], 
{i, 1, range}]

End[]

EndPackage[]