Interpolated Keeling Data

Joey Stevens

(* keeling=Import["\\\\fileserv2.hh.nku.edu\\stevensj8$\\Keeling.csv"]; *)
keeling=Import["http://ceadserv1.nku.edu/longa/classes/mat375/Mathematica/Keeling.csv"];
years=Delete[keeling[[All, 3]],1];
interpolated= Delete[keeling[[All, 5]], 1];

Basic Linear Regression

lp=ListPlot[interpolated]

linearModel= LinearModelFit[interpolated,x,x]

linearModelPlot=Plot[linearModel[x], {x, 0, 720}, PlotStyle->{Thick,Red}];
dataPlot=ListPlot[interpolated];
Show[linearModelPlot, dataPlot, Frame-> True, PlotRange-> All]

Diagnostics

linearModel["ParameterTable"]

linearModel["ParameterConfidenceIntervals"]

linearModel["RSquared"]

The linear model shows that a linear term (x) is necessary in the model and has an (R^2) value of 0.977756, but there is obviously some curvature to the data. It is more than likely a quadratic, or higher term, is necessary.

Linear Regression with quadratic term

lp

linearModelQuad= LinearModelFit[interpolated,{x,x^2},x]

linearModelQuadPlot= Plot[linearModelQuad[x], {x, 0, 720}, PlotStyle->{Thick,Red}];
dataPlot=ListPlot[interpolated];
Show[linearModelQuadPlot, dataPlot, Frame-> True, PlotRange-> All]

Diagnostics

linearModelQuad["ParameterTable"]

linearModelQuad["ParameterConfidenceIntervals"]

linearModelQuad["RSquared"]

Now, we see with the confidence intervals that a quadratic term is needed and subsequently the (R^2) value increased by roughly 1.5%. This model is getting closer but we still need to account for the oscillations in the data which means trigonometric functions are needed.

Linear Regression with Cos[x] and Sin[x]

linearModelTrig=LinearModelFit[interpolated, {x,x^2,Sin[x/365], Cos[x/365]},{x}]

linearModelTrigPlot= Plot[linearModelTrig[x], {x, 0, 720}, PlotStyle->{Thick,Red}];
dataPlot=ListPlot[interpolated];
Show[linearModelTrigPlot, dataPlot, Frame-> True, PlotRange-> All]

Diagnostics

linearModelTrig["ParameterTable"]
ListPlot[linearModelTrig["FitResiduals"], Filling->Axis,PlotLabel->"Residuals Plot for Detailed Keeling Data",Axes->True,AxesLabel->{"\!\(\*
StyleBox[\"Month\",\nFontWeight->\"Plain\"]\)","CO2 Levels\n(Parts per Million)"}]

linearModelTrig["ParameterConfidenceIntervals"]

linearModelTrig["RSquared"]

Both the Cosine and Sine terms are significant terms in the model and the (R^2) value has increased with the added terms.

Cubic Model

linearModelCubic=LinearModelFit[interpolated, {x,x^2,x^3,Sin[x/365], Cos[x/365]},{x}]

linearModelCubicPlot= Plot[linearModelCubic[x], {x, 0, 720}, PlotStyle->{Thick,Red}];
dataPlot=ListPlot[interpolated];
Show[linearModelCubicPlot, dataPlot, Frame-> True, PlotRange-> All]

Diagnostics

linearModelCubic["ParameterTable"]
ListPlot[linearModelCubic["FitResiduals"], Filling->Axis,PlotLabel->"Residuals Plot for Detailed Keeling Data",Axes->True,AxesLabel->{"\!\(\*
StyleBox[\"Month\",\nFontWeight->\"Plain\"]\)","CO2 Levels\n(Parts per Million)"}]

linearModelCubic["ParameterConfidenceIntervals"]
linearModelCubic["RSquared"]

The cubic term is not significant in the model, as well as the trigonometric functions, so this is probably not the best model to use.

Cubic Model -- alternate period (yearly)

linearModelCubic=LinearModelFit[interpolated, {x,x^2,x^3,Sin[2 Pi x / 12], Cos[2 Pi x / 12]},{x}]

linearModelCubicPlot= Plot[linearModelCubic[x], {x, 0, 720}, PlotStyle->{Thick,Red}];
dataPlot=ListPlot[interpolated];
Show[linearModelCubicPlot, dataPlot, Frame-> True, PlotRange-> All]

linearModelCubic["ParameterTable"]
ListPlot[linearModelCubic["FitResiduals"], Filling->Axis,PlotLabel->"Residuals Plot for Detailed Keeling Data\nCubic Model",Axes->True,AxesLabel->{"\!\(\*
StyleBox[\"Month\",\nFontWeight->\"Plain\"]\)","CO2 Levels\n(Parts per Million)"}]

linearModelCubic["ParameterConfidenceIntervals"]
linearModelCubic["RSquared"]

So what have we learned? Let’s go back to that quadratic trig model, and find out what it was “thinking”...

monthperiod=2.0 Pi 365
period=monthperiod/12 (* in years *)

linearModelTrig["ParameterTable"]

params=linearModelTrig["BestFitParameters"]
a=params[[4]]
b=params[[5]]

amplitude=Sqrt[a^2+b^2]
theta=ArcTan[b/a]

Let’s just check that we’ve produced an equivalent model:

p1=Plot[a Sin[2 Pi x /monthperiod]+b Cos[2 Pi x /monthperiod],{x,0,700}];
p2=Plot[1*amplitude Sin[2 Pi x /monthperiod+ theta],{x,0,700}];
Show[p1,p2]

Now let’s discover about when the carbon absorption peaks (and when it pits!):

yearpeak=(Pi/2-theta)/(2 Pi)*period
yearpit=(3 Pi/2-theta)/(2 Pi)*period