Keeling Analysis: 
non-linear regression

Andy Long
Spring, 2020

Here’s one way you might Import the CO2 data: straight off the web.

dat=Import["ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt","Table"];
DeleteCases[dat,{__String},{-1}];
dat[[96,1]]
dat[[71]]
Dimensions[dat]
(* First data line in the file: *)
Table[dat[[73,i]],{i,1,7}] 
dim=Dimensions[dat][[1]]
DecimalDate=Table[dat[[i,3]]-1958,{i,74,dim}] ;
CO2data=Table[dat[[i,4]],{i,74,dim}] ;
KeelingData=Transpose[{DecimalDate,CO2data}];
p1=ListPlot[KeelingData]

We’ll do some linear regression to get starting parameters:

Clear[a,b,c]
lm=LinearModelFit[
KeelingData,{x,x^2,Sin[2Pi x],Cos[2Pi x]},x]

{astart,bstart,cstart,d,e}=lm["BestFitParameters"]
ampstart=Sqrt[d^2+e^2]
x0start=-ArcTan[e/d]/(2 Pi)
lm["ParameterTable"]
lm["ANOVATable"]
Show[p1,Plot[lm[x],{x,0,Max[DecimalDate]}]]

Examining the RSquared as calculated by the non-linear regression function NonlinearModelFit, and that produced by LinearModelFit:

SSTuncorrected=(CO2data).(CO2data)
fits=Table[lm[x],{x,DecimalDate}];
SSModel=(fits).(fits)
SSModel/SSTuncorrected
lm["RSquared"]

Now we’ll do non-linear regression to get the parameters:

Clear[a,b,c,amp,x0,period]
{astart,bstart,cstart,ampstart,x0start}
nlm=NonlinearModelFit[
KeelingData,
a+b x+c x^2 +amp Sin[2 Pi (x-x0)/(1-period)],
{{a ,astart},{b,bstart},{c,cstart},
{amp,ampstart},{x0,x0start},{period,0}},
x]
{astart,bstart,cstart,ampstart,x0start}=nlm["BestFitParameters"]
nlm["ParameterTable"]
nlm["ANOVATable"]

Examining the RSquared as calculated by the non-linear regression function NonlinearModelFit:

fits=Table[nlm[x],{x,DecimalDate}];
SSModel=(fits).(fits)
SSModel/SSTuncorrected
nlm["RSquared"]

Show[p1,Plot[nlm[x],{x,0,Max[DecimalDate]}]]

Not much change. Now let’s allow the amplitude of the sinusoidal variation be a linear function, and see if it’s significantly changing over time:

Clear[a,b,c,amp,x0,period,m]
{astart,bstart,cstart,ampstart,x0start}
nlmFinal=NonlinearModelFit[
KeelingData,
a+b x+c x^2 +(amp+m x) Sin[2 Pi (x-x0)/(1-period)],
{{a ,astart},{b,bstart},{c,cstart},
{amp,ampstart},{x0,x0start},{period,0},{m,0}},
x]
{aend,bend,cend,ampend,x0end,periodend, mend}=nlmFinal["BestFitParameters"]
nlmFinal["ParameterTable"]

Examining the RSquared as calculated by the non-linear regression function NonlinearModelFit:

fits=Table[nlmFinal[x],{x,DecimalDate}];
SSModel=(fits).(fits)
SSModel/SSTuncorrected
nlmFinal["RSquared"]

Show[p1,Plot[nlmFinal[x],{x,0,Max[DecimalDate]}]]
ListPlot[Transpose[{DecimalDate,nlmFinal["FitResiduals"]}]]
phase=x0/.x0end
phaseInDaysFromNewYear=phase*365.25
periodFinal=1+period/.periodend
periodFinal= periodFinal*365.25

Let’s compare the models:

Plot[nlmFinal[x]-lm[x],{x,0,Max[DecimalDate]}]
Plot[nlmFinal[x]-nlm[x],{x,0,Max[DecimalDate]}]
Plot[nlm[x]-lm[x],{x,0,Max[DecimalDate]}]

Plot[nlmFinal[x]-lm[x],{x,0,10}]
Plot[nlmFinal[x]-nlm[x],{x,0,10}]
Plot[nlm[x]-lm[x],{x,0,10}]

A manipulate command let’s us see how the models are doing in their fit to the data:

Manipulate[
Module[{},
Show[
Plot[
{nlmFinal[x],nlm[x],lm[x]},
{x,start,start+int}]
,p1]
],
{start,0,Max[DecimalDate]},
{int,1,10}
]

I didn’t ask you for the exponential model, but Flerlage, et al. did a nice job in their report, and concluded that this was the best model:

“Mathematica found that the best power model fit for the data is fit by the nonlinear model: CO2=257.184 +56.2922 e^ (0.0163503 * (years since 1958)) + 2.83537 * sin(2π ((years since 1958)-1.05748)] )”

Clear[a,b,c,amp,x0,period]
{astart,bstart,cstart,ampstart,x0start,periodstart}={257.184,56.2922,0.0163503,2.83537,.05748,0}
nlmExpo=NonlinearModelFit[
KeelingData,
a+b Exp[c x] +amp Sin[2 Pi (x-x0)/(1-period)],
{{a ,astart},{b,bstart},{c,cstart},
{amp,ampstart},{x0,x0start},{period,0}},
x]
{aend,bend,cend,ampend,x0end,periodend, mend}=nlmExpo["ParameterTable"]
nlmExpo["ANOVATable"]

Show[p1,Plot[nlmExpo[x],{x,0,Max[DecimalDate]}]]
phase=x0/.x0end
phaseInDaysFromNewYear=phase*365.25
periodFinal=1+period/.periodend
periodFinal= periodFinal*365.25

Examining the RSquared as calculated by the non-linear regression function NonlinearModelFit:

fits=Table[nlmExpo[x],{x,DecimalDate}];
SSModel=(fits).(fits)
SSModel/SSTuncorrected
nlmExpo["RSquared"]

Plot[nlmFinal[x]-nlmExpo[x],{x,0,Max[DecimalDate]}]
NSolve[nlmExpo[x]==450,x]
NSolve[nlm[x]==450,x]
NSolve[lm[x]==450,x]

NSolve[254.99001926392936` +58.251194601701116` (E^(0.01603937562541139` x))+2.846689094213021` Sin[6.28608558285868` (-0.07262448861204018`+x)]==450,x]

ListPlot[Transpose[{DecimalDate,nlmExpo["FitResiduals"]}]]

1958+x/.FindRoot[nlmExpo[x]-450,{x,80}]
1958+x/.FindRoot[nlm[x]-450,{x,80}]
1958+x/.FindRoot[lm[x]-450,{x,80}]