Togo Models for Temperatures and Rainfall with ENSO, SST, and TSI

Andy Long
Spring, 2018

Setting Up:

Importing Data

city="Togo";datahome="/home/andy/classes/2018Spring/MathModeling/projects/Final/Model/";
datahome="http://ceadserv1.nku.edu/longa/classes/mat375/finalproject/Data/v2/" ; (* web *)
(* I had to remove missing values from the V2 file *)
file=StringJoin[datahome,"All_Clean_v2_missing_removed.xlsx"];
list=Import[file];

list[[1,1,All]]

rainfall=Rest[Flatten[list[[1,All,4]]]];
location=Rest[Flatten[list[[1,All,1]]]];
longitude=Rest[Flatten[list[[1,All,6]]]];
latitude=Rest[Flatten[list[[1,All,7]]]];
elevation=Rest[Flatten[list[[1,All,8]]]];
decYear=Rest[Flatten[list[[1,All,10]]]];
(* Note that we're shifting time to 0 at the beginning of year 1961, the start of our study. This changes constant terms in models, essentially, but generally makes linear regression better conditioned: *)
decYear=decYear-1961;
minTemp=Rest[Flatten[list[[1,All,11]]]];
maxTemp=Rest[Flatten[list[[1,All,5]]]];
nvals=Length[rainfall];

Bring in the ENSO data, that I got from the Aussies:

ftp://ftp.bom.gov.au/anon/home/ncc/www/sco/soi/soiplaintext.html

ENSO=Import[StringJoin[datahome,"ENSO.txt"],"Data"];
moVals={0.04166666666992569,0.125,0.2083333333300743,0.2916666666699257,0.375,0.4583333333300743,0.5416666666699257,.625,0.7083333333300743,0.7916666666699257,0.875,0.9583333333300743};
ENSOdates=Flatten[Table[i+moVals[[j]],{i,0,54},{j,1,12}]];
ENSO=Table[First[ENSO[[j]]],{j,1,660}];
ENSO=ENSO-Mean[ENSO];
ENSOdata=Transpose[{ENSOdates,ENSO}];

Because we’d deleted some dates, and since some time series were shorter than others, I was obliged to do some ugly fenegling to get the dates to line up right. Don’t tell anyone about this... I’m ashamed:

ENSOvec={};
Do[
Do[
If[Abs[ENSOdata[[i,1]]-decYear[[j]]]<.01,
ENSOvec=AppendTo[ENSOvec,ENSOdata[[i,2]]]],{i,1,660}],
{j,1,nvals}];

ListPlot[ENSOdata]

https://climexp.knmi.nl/get_index.cgi

Bring in the SST data, that I got from the Dutch:

https://climexp.knmi.nl/select.cgi?id=someone@somewhere&field=hadisst1

I chose 
Get grid points, average area or generate subset
latitude 0 to 6 
(so 6 degrees down to the equator from the coast of Togo, and
longitude -4 to 6
4 degrees to the left of Togo, 4 degrees to the right.)

Then “make time series”. I grabbed the “raw data” from the result, and that’s what became the file SST.txt which I’ve got on my website....

SST=Import[StringJoin[datahome,"SST.txt"],"Data"];
SSTdates=Flatten[Table[i+moVals[[j]],{i,0,54},{j,1,12}]];
SST=Table[First[SST[[j]]],{j,1,660}];
SST=SST-Mean[SST];
SSTdata=Transpose[{SSTdates,SST}];

Because we’d deleted some dates, and since some time series were shorter than others, I was obliged to do some ugly fenegling to get the dates to line up right. Don’t tell anyone about this... I’m ashamed:

SSTvec={};
Do[
Do[
If[Abs[SSTdata[[i,1]]-decYear[[j]]]<.01,
SSTvec=AppendTo[SSTvec,SSTdata[[i,2]]]],{i,1,660}],
{j,1,nvals}];

ListPlot[SSTdata]

Bring in the TSI data, that I got from the Colorado:

http://lasp.colorado.edu/lisird/data/nrl2_tsi_P1M/

TSI=Import[StringJoin[datahome,"nrl2_tsi_P1M.csv"],"Data"];
Dimensions[TSI]
TSIdates=Flatten[Table[i+moVals[[j]],{i,0,54},{j,1,12}]];
TSIvec=Rest[TSI[[All,2]]];
TSIvec=TSIvec-Mean[TSIvec];
TSIdata=Transpose[{TSIdates,TSIvec}];

Because we’d deleted some dates, and since some time series were shorter than others, I was obliged to do some ugly fenegling to get the dates to line up right. Don’t tell anyone about this... I’m ashamed:

TSIvec={};
Do[
Do[
If[Abs[TSIdata[[i,1]]-decYear[[j]]]<.01,
TSIvec=AppendTo[TSIvec,TSIdata[[i,2]]]],{i,1,660}],
{j,1,nvals}];

ListPlot[TSIdata]

A couple of functions we need:

OurSquared replaces R-squared when we move to a model without an intercept

OurSquared[real_,pred_]:= 
Module[{mean=Mean[real],nvals=Length[real]},
1-((∑_(i=1)^nvals)(((real[[i]]-pred[[i]]))^2))/((∑_(i=1)^nvals)(((real[[i]]-mean))^2))
]

TestResiduals does just that -- for normality.

TestResiduals[resids_]:=
Module[
{s=StandardDeviation[resids],
m=Mean[resids],𝒹,lp,table,show},
𝒹=NormalDistribution[m,s];
lp=ListPlot[resids,PlotLabel->city];
DistributionFitTest[resids];
ℋ=DistributionFitTest[resids,Automatic,"HypothesisTestData"];
show=Show[
Histogram[resids,Automatic,"ProbabilityDensity"],Plot[PDF[ℋ["FittedDistribution"],x],{x,-5,5},PlotStyle->Thick]
];
table=DistributionFitTest[resids,𝒹,{"TestDataTable","AndersonDarling"}];
GraphicsGrid[{{lp,show},{table,
ProbabilityPlot[resids,𝒹]}}]
]

Global Togo Model for Minimum Temperature

We generate the Minimum Temperature Model:

finalmin=LinearModelFit[Transpose[{ENSOvec,minTemp}],{x},{x}]
finalmin["ParameterTable"]
finalmin["RSquared"]

mindata=Transpose[{latitude,longitude,elevation,ENSOvec,SSTvec,TSIvec,decYear,minTemp}];

finalmin=LinearModelFit[mindata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/5)],Cos[2*π*x/(1/5)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
lat,
long,
elev,
enso,
sst,
tsi
},{lat,long,elev,enso,sst,tsi,x}]
finalmin["ParameterTable"]
finalmin["AdjustedRSquared"]

finalmin=LinearModelFit[mindata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
enso,
sst,
tsi
},{lat,long,elev,enso,sst,tsi,x}]
finalmin["ParameterTable"]
finalmin["AdjustedRSquared"]

The minimum residuals are tested against a normal, and are not normal:

minResiduals=finalmin["FitResiduals"];
TestResiduals[minResiduals]

Minmax=Max[minTemp];
Minmin=Min[minTemp];
Minavg=Mean[Minmax,Minmin];
Manipulate[
Module[
{
end=t+n*12,
yearmin,
yearmax,
p1, p2
},
yearmin=Min[Part[decYear,t;;end]];
yearmax=Max[Part[decYear,t;;end]];
p1=ListPlot[Table[{decYear[[i]],finalmin[latitude[[i]],longitude[[i]],elevation[[i]],ENSOvec[[i]],SSTvec[[i]],TSIvec[[i]],decYear[[i]]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{Minmin,Minmax}},Joined->True,PlotLabel->{location[[t]],decYear[[t]]},AxesLabel->{"Year","Min"}];
p2=ListPlot[Table[{decYear[[i]],minTemp[[i]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{Minmin,Minmax}},PlotStyle->PointSize[Medium],PlotLabel->{location[[t]],decYear[[t]]},AxesLabel->{"Year","Min"}];
Show[p1,p2]
],
{{t,1},1,nvals-n*12,12},
{{n,1},1,5,1}
]

Global Togo Model for Maximum Temperature

We generate the Maximum Temperature Model:

finalmax=LinearModelFit[Transpose[{ENSOvec,SSTvec,maxTemp}],{x,y},{x,y}]
finalmax["ParameterTable"]
finalmax["RSquared"]

maxdata=Transpose[{latitude,longitude,elevation,ENSOvec,SSTvec,TSIvec,decYear,maxTemp}];

finalmax=LinearModelFit[maxdata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
enso,
sst,
tsi
},{lat,long,elev,enso,sst,tsi,x}]
finalmax["ParameterTable"]
finalmax["AdjustedRSquared"]

The maximum residuals are tested against a normal, and are not bad!

maxResiduals=finalmax["FitResiduals"];
TestResiduals[maxResiduals]

Maxmax=Max[maxTemp];
Maxmin=Min[maxTemp];
Maxavg=Mean[Maxmax,Maxmin];
Manipulate[
Module[
{
end=t+n*12,
yearmin,
yearmax,
p1, p2
},
yearmin=Min[Part[decYear,t;;end]];
yearmax=Max[Part[decYear,t;;end]];
p1=ListPlot[Table[{decYear[[i]],finalmax[latitude[[i]],longitude[[i]],elevation[[i]],ENSOvec[[i]],SSTvec[[i]],TSIvec[[i]],decYear[[i]]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{Maxmin,Maxmax}},Joined->True,PlotLabel->{location[[t]],decYear[[t]]},AxesLabel->{"Year","Max"}];
p2=ListPlot[Table[{decYear[[i]],maxTemp[[i]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{Maxmin,Maxmax}},PlotStyle->PointSize[Medium],PlotLabel->{location[[t]],decYear[[t]]},AxesLabel->{"Year","Max"}];
Show[p1,p2]
],
{{t,1},1,nvals-n*12,12},
{{n,1},1,5,1}
]

Global Togo Model for Rainfall Data

Time now for the rainfall: let’s see how it depends on max and min temperatures:

finalrain=LinearModelFit[Transpose[{ENSOvec,SSTvec,rainfall}],{x,y},{x,y}]
finalrain["ParameterTable"]
finalrain["RSquared"]

raindata=Transpose[{decYear,latitude,longitude,elevation,minTemp,maxTemp,ENSOvec,SSTvec,TSIvec,rainfall}];

finalrain=LinearModelFit[raindata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
mint,
maxt,
enso,
sst,
tsi
},{x,lat,long,elev,mint,maxt,enso,sst,tsi}]
finalrain["ParameterTable"]

Since the constant term is not significant, it seems wise to drop it. I haven’t done so just yet, however. There are only modest changes to the other terms:

finalrainNoConstant=LinearModelFit[raindata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
mint,
maxt,
enso,
sst,
tsi
},{x,lat,long,elev,mint,maxt,enso,sst,tsi},IncludeConstantBasis->False]
finalrainNoConstant["ParameterTable"]

Note that when we drop the constant term, the R^2 rockets up. We don’t really believe that, however, so we consult “Our^2”, and it shows that they’re really about the same:

finalrain["AdjustedRSquared"]
finalrainNoConstant["AdjustedRSquared"]
(* Note the ugly -- but easy -- way to compute the estimates *)
OurSquared[rainfall,rainfall-finalrain["FitResiduals"]]
OurSquared[rainfall,rainfall-finalrainNoConstant["FitResiduals"]]

rainestimates=rainfall-finalrainNoConstant["FitResiduals"];
newrain=HeavisideTheta[rainestimates]*rainestimates;
rainresiduals=rainfall-newrain;
OurSquared[rainfall,newrain]
ListPlot[Transpose[{rainfall, newrain}]]
ListPlot[Transpose[{rainfall, rainresiduals}]]
ListPlot[Transpose[{rainfall, rainresiduals}]]
TestResiduals[rainResiduals]
ListPlot[Transpose[{decYear,newrain}]]

Rainmin=Min[rainfall];
Rainmax=Max[rainfall];
Rainavg=Mean[{Rainmin,Rainmax}];
Manipulate[
Module[
{
end=t+n*12,
yearmin,
yearmax,
p1, p2
},
yearmin=Min[Part[decYear,t;;end]];
yearmax=Max[Part[decYear,t;;end]];
p1=ListPlot[Table[{decYear[[i]],finalrainNoConstant[decYear[[i]],latitude[[i]],longitude[[i]],elevation[[i]],minTemp[[i]],maxTemp[[i]],ENSOvec[[i]],SSTvec[[i]],TSIvec[[i]]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{Rainmin,Rainmax}},Joined->True,PlotLabel->{location[[t]],decYear[[t]]},AxesLabel->{"Year","Rain"}];
p2=ListPlot[Table[{decYear[[i]],rainfall[[i]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{Rainmin,Rainmax}},PlotStyle->PointSize[Medium],PlotLabel->{location[[t]],decYear[[t]]},AxesLabel->{"Year","Rain"}];
Show[p1,p2]
],
{{t,1},1,nvals-n*12,12},
{{n,1},1,5,1}
]

The rainfall residuals are not normal :(                :

rainResiduals=finalrainNoConstant["FitResiduals"];
TestResiduals[rainResiduals]

sstdata=Transpose[{decYear,SSTvec}];
sstReg=LinearModelFit[sstdata,{
x,
x^2,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)]
},{x}]
sstReg["ParameterTable"]
sstReg["AdjustedRSquared"]

sstdata=Transpose[{decYear,SSTvec}];
sstReg=LinearModelFit[sstdata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)]
},{x}]
sstReg["ParameterTable"]
sstReg["AdjustedRSquared"]

sstdata=Transpose[{decYear,TSIvec,SSTvec}];
sstReg=LinearModelFit[sstdata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)],
tsi
},{x,tsi}]
sstReg["ParameterTable"]
sstReg["AdjustedRSquared"]

SSTparams=sstReg["BestFitParameters"];
SSTFn[x_]:=SSTparams.{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)]
}

Manipulate[
Module[
{
end=t+n*12,
SSTmin,
SSTmax,
SSTavg,
yearmin,
yearmax,
p1, p2
},
yearmin=Min[Part[decYear,t;;end]];
yearmax=Max[Part[decYear,t;;end]];
SSTmin=Min[Part[SSTvec,t;;end]];
SSTmax=Max[Part[SSTvec,t;;end]];
SSTavg=Mean[{SSTmin,SSTmax}];
SSTmin=SSTavg-1.5(SSTavg-SSTmin);
SSTmax=SSTavg+1.5(SSTavg-SSTmin);
p1=ListPlot[Table[{decYear[[i]],sstReg[decYear[[i]],TSIvec[[i]]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{SSTmin,SSTmax}},Joined->True,PlotLabel->decYear[[t]],AxesLabel->{"Year","SST"}];
p2=ListPlot[Table[{decYear[[i]],SSTvec[[i]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{SSTmin,SSTmax}},PlotStyle->PointSize[Medium],PlotLabel->decYear[[t]],AxesLabel->{"Year","SST"}];
Show[p1,p2,PlotRange->All]
],
{{t,1},1,nvals-n*12,12},
{{n,1},1,5,1}
]

sstReg[0,TSIvec[[1]]]

I just want to take a look at all 147 years of SST data, to see if those periods show up over all, of if they are perhaps artifactual....

The first plot is recent times -- 1961 -- 2015

the second is for the whole time, from 1870 to present.