Keeling Data and Oscillatory Behavior

Now let’s add in the sinusoidal variation of the Keeling data:

keeling=Import["http://ceadserv1.nku.edu/longa/classes/2018spring/mat375/mathematica/Keeling.csv"];
y=Rest[Transpose[keeling][[5]]];
t=Rest[Transpose[keeling][[3]]];
tstart=1958; (* let's establish the start date, then model from there: *)
t=t-tstart;
lm=LinearModelFit[Transpose[{t,y}],{x,x^2,x^3,Sin[4 Pi x],Cos[4 Pi x],Sin[2 Pi x],Cos[2 Pi x],Sin[2 Pi (x /8)],Cos[2 Pi (x /8)],Sin[2 Pi (x /16)],Cos[2 Pi (x /16)],Sin[2 Pi (x /32)],Cos[2 Pi (x /32)],Sin[2 Pi (x /64)],Cos[2 Pi (x /64)]},x]

lm["ParameterTable"]
lm["RSquared"]
beta=lm["BestFitParameters"]

ListPlot[lm["FitResiduals"]]
(* Does this model look better? *)
model[x_]:=beta.{1,x,x^2,x^3,Sin[4 Pi x],Cos[4 Pi x],Sin[2 Pi x],Cos[2 Pi x],Sin[2 Pi (x /8)],Cos[2 Pi (x /8)],Sin[2 Pi (x /16)],Cos[2 Pi (x /16)],Sin[2 Pi (x /32)],Cos[2 Pi (x /32)],Sin[2 Pi (x /64)],Cos[2 Pi (x /64)]}
model[t[[1]]]
n=Length[y];
p1=ListPlot[Table[{t[[i]],y[[i]]},{i,1,n}]];
p2=Plot[model[x],{x,Min[t],Max[t]}]
Show[p1,p2]
(* Let’s have a look at the residuals: *)
ListPlot[Table[{t[[i]],y[[i]]-model[t[[i]]]},{i,1,n}]]

Something’s going on in those residuals, but who knows what? Those are trends over 20 year spans of some sort....

p1=ListPlot[Transpose[{t,y}]];
Length[t]
Length[y]
p2=Plot[model[x],{x,Min[t],Max[t]}]
Show[p1,p2]
(* Let’s have a look at the residuals: *)
ListPlot[Table[{t[[i]],y[[i]]-model[t[[i]]]},{i,1,n}]]

We might like to scan across the years, to see how well this model fits at various points:

(* Does this model look better? *)
p1=ListPlot[Table[{t[[i]]+tstart,y[[i]]},{i,1,n}]];
Manipulate[Show[Plot[model[x-tstart],{x,a,a+5}],p1],{a,1959,2014}]

Now we’re going to find the day of the peak of the oscillation, and the size of the amplitude of the oscillation.

psine=beta[[5]]
pcosine=beta[[6]]
amplitude=Sqrt[psine^2+pcosine^2]
theta=ArcTan[pcosine/psine]
Plot[{psine*Sin[2Pi t]+pcosine*Cos[2 Pi t],amplitude*Sin[theta+2Pi t]},{t,0,2}]
soln=Solve[theta+2Pi x==Pi/2,x]
peak=tstart+x/.soln[[1]]
(peak-Floor[peak])*365.25 (* About April 23rd *)

rsquareValues={};
sValues={};
For[s=0,s<5,s=s+.1,
lm=LinearModelFit[Transpose[{t,y}],{x,x^2,x^3,Sin[ Pi x/(1+s x)],Cos[ Pi x/(1+s x)]},x];
rsquareValues=Append[rsquareValues,lm["AdjustedRSquared"]];
sValues=Append[sValues,s];
]
ListPlot[Transpose[{sValues,rsquareValues}],PlotRange->All]

marker1=Graphics[{Red,Disk[]}];
Manipulate[
Module[
{maxrsquared=0,
rsquared,
rsquareValues={},
sValues={}},
For[s=0.01,s<10,s=s+.1,
lm=LinearModelFit[Transpose[{t,y}],{x,x^2,x^3,Sin[2 Pi x],Cos[2 Pi x],Sin[2 Pi x/3.75],Cos[2 Pi x/3.75],Sin[2 Pi x/30],Cos[2 Pi x/30],Sin[2 Pi x/16],Cos[2 Pi x/16],Sin[  c Pi x/(1+s)],Cos[ c Pi x/(1+s)]},x];
lm=LinearModelFit[Transpose[{t,y}],{x,x^2,x^3,Sin[2 Pi x],Cos[2 Pi x],Sin[2 Pi x/4],Cos[2 Pi x/4],Sin[2 Pi x/32],Cos[2 Pi x/32],Sin[2 Pi x/16],Cos[2 Pi x/16],Sin[  c Pi x/(1+s)],Cos[ c Pi x/(1+s)]},x];
rsquared=lm["AdjustedRSquared"];
rsquareValues=Append[rsquareValues,rsquared];
If[rsquared>maxrsquared,maxrsquared=rsquared];
sValues=Append[sValues,s];
];
ListPlot[Transpose[{sValues,rsquareValues}],PlotRange->All,PlotLabel->maxrsquared]
]
,{{c,1},0.1,20,.1} 
]

So here’s my best model:

lm=LinearModelFit[Transpose[{t,y}],{x,x^2,Sin[2 Pi x],Cos[2 Pi x],Sin[2 Pi x/32],Cos[2 Pi x/32],Sin[2 Pi x/16],Cos[2 Pi x/16]},x];
lm["AdjustedRSquared"]
lm["ParameterTable"]
ListPlot[lm["FitResiduals"]]
bestbeta=lm["BestFitParameters"];
model[x_]:=bestbeta.{1,x,x^2,Sin[2 Pi x],Cos[2 Pi x],Sin[2 Pi x/32],Cos[2 Pi x/32],Sin[2 Pi x/16],Cos[2 Pi x/16]}
model[t[[1]]]
n=Length[y];
p1=ListPlot[Table[{t[[i]],y[[i]]},{i,1,n}]];
p2=Plot[model[x],{x,Min[t],Max[t]}]
Show[p1,p2]

rsquareValues={};
sValues={};
For[s=0.001,s<20,s=s+.05,
lm=LinearModelFit[Transpose[{t,y}],{x,x^2,x^3,Sin[Pi x],Cos[Pi x],Sin[Pi x/(1+s )],Cos[Pi x/(1+s )]},x];
rsquareValues=Append[rsquareValues,lm["AdjustedRSquared"]];
sValues=Append[sValues,s];
]
plotrsq=ListPlot[Transpose[{sValues,rsquareValues}],PlotRange->All]

marker1=Graphics[{Red,Disk[]}];
Manipulate[
Module[
{lm,rsquared},
lm=LinearModelFit[Transpose[{t,y}],{x,x^2,x^3,Sin[Pi x],Cos[Pi x],Sin[Pi x/(1+s )],Cos[Pi x/(1+s )]},x];
rsquared=lm["AdjustedRSquared"];
Show[plotrsq,ListPlot[{{s,rsquared}},PlotMarkers->{marker1,.045},PlotLabel->rsquared]]
]
,{{s,14.7},0.001,20,.005} (* 6.2, 14.7 *)
]