Sally and Connor’s Sotouboua analysis file, which I modified to extend some of their good ideas!

Andy Long
Spring, 2018

rainfalldata=Import["http://ceadserv1.nku.edu/longa/classes/mat375/data/BinaoData/Pluviometry/TOTAUXPLUVIOMETRIQUESMENSUELS1961-2015.xls"];
Rainfall=Flatten[Table[rainfalldata[[5,i]][[j]],{i,26,59},{j,3,14}]];
mintempdata=Import["http://ceadserv1.nku.edu/longa/classes/mat375/data/BinaoData/RawMonthly/TOGO-TEMPERATUREMOYENNEMINIMALE61-15.xls"];
maxtempdata=Import["http://ceadserv1.nku.edu/longa/classes/mat375/data/BinaoData/RawMonthly/TOGO-TEMPERATUREMOYENNEMAXIMALE61-15.xls"];

Length[Rainfall]

Maxtempyears=Table[maxtempdata[[5,i]][[2]],{i,8,41}];
Mintempyears=Table[mintempdata[[5,i]][[2]],{i,7,40}];

Mintempyears=Mintempyears-1982;
Maxtempyears=Maxtempyears-1982;

Mintemp=Flatten[Table[mintempdata[[5,i]][[j]],{i,7,40},{j,3,14}]];
Maxtemp=Flatten[Table[maxtempdata[[5,i]][[j]],{i,8,41},{j,3,14}]];

Mindecyears=Flatten[Table[Mintempyears[[i]]+(j/12),{i,1,34},{j,1,12}]];
Maxdecyears=Flatten[Table[Maxtempyears[[i]]+(j/12),{i,1,34},{j,1,12}]];

Functions preceeded by “Filter” had two points erased due to missing two months in the minimum and maximum temperature data, the months missing were January and February of 2013.

Filtermint=Delete[Mintemp,{{373},{374}}];
Filtermaxt=Delete[Maxtemp,{{373},{374}}];
Filterminy=Delete[Mindecyears,{{373},{374}}];
Filtermaxy=Delete[Maxdecyears,{{373},{374}}];
Filterrain=Delete[Rainfall,{{373},{374}}];

𝒟=SmoothKernelDistribution[Transpose[{Filtermint,Filtermaxt}]];
contour=ContourPlot[PDF[𝒟,{x,y}]//Evaluate,{x,Min[Filtermint],Max[Filtermint]},{y,Min[Filtermaxt],Max[Filtermaxt]},PlotRange->All,AxesLabel->{"Min","Max"},PlotPoints->50,Contours->15,ColorFunction->"TemperatureMap"]

dens=Plot3D[PDF[𝒟,{x,y}]//Evaluate,{x,Min[Filtermint],Max[Filtermint]},{y,Min[Filtermaxt],Max[Filtermaxt]},PlotRange->All,PlotPoints->50,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","PDF"}]
ListPointPlot3D[Transpose[{Filtermint,Filtermaxt,Filterrain}],PlotRange->All,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","Rainfall"}]

raintable=Table[{i+1981,j,Rainfall[[12*(i-1)+j]]},{i,1,34,1},{j,1,12,1}];
ListPointPlot3D[raintable,PlotRange->All,ColorFunction->"TemperatureMap",AxesLabel->{"Year","Month of the Year","Rainfall"}]

days=Import["http://www.nku.edu/~longa/classes/2018spring/mat375/mini-projects/mini-project-no3/Mango/Mango.dates.txt","List"];
days=Rest[days];

Filterdays=Drop[days,252];

Filterdays=Delete[Filterdays,{{373},{374}}];

minmodel=LinearModelFit[Transpose[{Filterdays,Filtermint}],{1,x,
Sin[2Pi (1 x/365.25)],Cos[2Pi (1 x/365.25)],
Sin[2Pi (2 x/365.25)],Cos[2Pi (2 x/365.25)]
},x];
maxmodel=LinearModelFit[Transpose[{Filterdays,Filtermaxt}],{1,x,
Sin[2Pi (1 x/365.25)],Cos[2Pi (1 x/365.25)],
Sin[2Pi (2 x/365.25)],Cos[2Pi (2 x/365.25)],
Sin[2Pi (4 x/365.25)],Cos[2Pi (4x/365.25)]},x];rainmodel=LinearModelFit[Transpose[{Filterdays,Filterrain}],{1,x,
Sin[2Pi (1 x/365.25)],Cos[2Pi (1x/365.25)],
Sin[2Pi (2 x/365.25)],Cos[2Pi (2 x/365.25)]},x];
minmodel["ParameterTable"]
maxmodel["ParameterTable"]
rainmodel["ParameterTable"]

minmodel["RSquared"]
maxmodel["RSquared"]
rainmodel["RSquared"]
minmodel["ParameterConfidenceIntervals"]
maxmodel["ParameterConfidenceIntervals"]
rainmodel["ParameterConfidenceIntervals"]

raindataplot=ListPlot[Transpose[{Filterdays,Filterrain}]];
Manipulate[
rainmodelplot=Plot[rainmodel[x],{x,t*365.25,(t+n)*365.25},PlotRange->{{t*365.25,(t+n)*365.25},{Min[Filterrain],Max[Filterrain]}},PlotLabel->"Model Versus Rain Data",AxesLabel->{"Days","Rainfall"}];
Show[rainmodelplot,raindataplot],
{{t,21},21,50,1},
{{n,1},1,5,1}]

maxtempdataplot=ListPlot[Transpose[{Filterdays,Filtermaxt}]];
Manipulate[
maxmodelplot=Plot[maxmodel[x],{x,t*365.25,(t+n)*365.25},PlotRange->{{t*365.25,(t+n)*365.25},{Min[Filtermaxt],Max[Filtermaxt]}},PlotLabel->"Model Versus Max Temp Data",AxesLabel->{"Days","Max Temp"}];
Show[maxmodelplot,maxtempdataplot],
{{t,21},21,50,1},
{{n,1},1,5,1}]

mintempdataplot=ListPlot[Transpose[{Filterdays,Filtermint}]];
Manipulate[
minmodelplot=Plot[minmodel[x],{x,t*365.25,(t+n)*365.25},PlotRange->{{t*365.25,(t+n)*365.25},{Min[Filtermint],Max[Filtermint]}},PlotLabel->"Model Versus Min Temp Data",AxesLabel->{"Days","Min Temp"}];
Show[minmodelplot,mintempdataplot],
{{t,21},21,50,1},
{{n,1},1,5,1}]

All functions proceeded by “no” are data with No Outliers.

Rainno=Delete[Filterrain,{{21},{37},{56},{59},{60},{61},{73},{84},{106},{108},{131},{139},{149},{160},{168},{193},{196},{250},{256},{258},{259},{260},{261},{262},{264},{265},{267},{268},{269},{270},{271},{275},{276},{289},{300},{312},{372},{383},{388}}];
Length[Rainno]

Mintempno=Delete[Filtermint,{{21},{37},{56},{59},{60},{61},{73},{84},{106},{108},{131},{139},{149},{160},{168},{193},{196},{250},{256},{258},{259},{260},{261},{262},{264},{265},{267},{268},{269},{270},{271},{275},{276},{289},{300},{312},{372},{383},{388}}];
Length[Mintempno]

Maxtempno=Delete[Filtermaxt,{{21},{37},{56},{59},{60},{61},{73},{84},{106},{108},{131},{139},{149},{160},{168},{193},{196},{250},{256},{258},{259},{260},{261},{262},{264},{265},{267},{268},{269},{270},{271},{275},{276},{289},{300},{312},{372},{383},{388}}];
Length[Maxtempno]

Daysno=Delete[Filterdays,{{21},{37},{56},{59},{60},{61},{73},{84},{106},{108},{131},{139},{149},{160},{168},{193},{196},{250},{256},{258},{259},{260},{261},{262},{264},{265},{267},{268},{269},{270},{271},{275},{276},{289},{300},{312},{372},{383},{388}}];
Length[Daysno]

minmodel=LinearModelFit[Transpose[{Daysno,Mintempno}],{1,x,
Sin[2Pi (1 x/365.25)],Cos[2Pi (1 x/365.25)],
Sin[2Pi (2 x/365.25)],Cos[2Pi (2 x/365.25)]
},x];
maxmodel=LinearModelFit[Transpose[{Daysno,Maxtempno}],{1,x,
Sin[2Pi (1 x/365.25)],Cos[2Pi (1 x/365.25)],
Sin[2Pi (2 x/365.25)],Cos[2Pi (2 x/365.25)],
Sin[2Pi (4 x/365.25)],Cos[2Pi (4x/365.25)]},x];
rainmodel=LinearModelFit[Transpose[{Daysno,Rainno}],{1,
Sin[2Pi (1 x/365.25)],Cos[2Pi (1x/365.25)],
Sin[2Pi (2 x/365.25)],Cos[2Pi (2 x/365.25)]},x]
minmodel["ParameterTable"]
maxmodel["ParameterTable"]
rainmodel["ParameterTable"]

minmodel["RSquared"]
maxmodel["RSquared"]
rainmodel["RSquared"]
minmodel["ParameterConfidenceIntervals"]
maxmodel["ParameterConfidenceIntervals"]
rainmodel["ParameterConfidenceIntervals"]

𝒟=SmoothKernelDistribution[Transpose[{Mintempno,Maxtempno}]];
contour=ContourPlot[PDF[𝒟,{x,y}]//Evaluate,{x,Min[Mintempno],Max[Mintempno]},{y,Min[Maxtempno],Max[Maxtempno]},PlotRange->All,PlotLabel->"Min. v. Max. as a function of Rainfall",PlotPoints->50,Contours->15,ColorFunction->"TemperatureMap"]

dens=Plot3D[PDF[𝒟,{x,y}]//Evaluate,{x,Min[Mintempno],Max[Mintempno]},{y,Min[Maxtempno],Max[Maxtempno]},PlotRange->All,PlotPoints->50,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","PDF"}]
ListPointPlot3D[Transpose[{Mintempno,Maxtempno,Rainno}],PlotRange->All,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","Rainfall"}]

ael: tempDiffModel is a model that I did, following the lead of Connor and Sally: only I used temperature difference and sinusoidal variation in time (which is not incorporated into the model below). I also did mine on the unfiltered data -- they did a lot of processing, to remove suspected outliers. 

I just want to show the impact of the quadratic temperature difference term (the time dependence is just going to place the modeled point at the proper place on the following quadratic surface):

tempDiffModel[min_,max_]:=355.441+(max-min)(-29.5096 +0.716411*(max-min))
tempDiffModelPlot=Plot3D[tempDiffModel[x,y],{x,Min[Mintempno],Max[Mintempno]},{y,Min[Maxtempno],Max[Maxtempno]},PlotRange->{{Min[Mintempno],Max[Mintempno]},{Min[Maxtempno],Max[Maxtempno]}},PlotPoints->50,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","Rain Model"}];
dataPlot=ListPointPlot3D[Transpose[{Mintempno,Maxtempno,Rainno}],PlotRange->All,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","Rainfall"},PlotStyle->PointSize[0.03]];
Show[dataPlot,tempDiffModelPlot]

ael: Now I’ll do the same using their data (and we note that the linear term in time is not significant):

u=ConstantArray[1,Length[Maxtempno]];
v1=Maxtempno-Mintempno;
v2=v1^2;
v3=Daysno;
v4=Sin[2Pi Daysno/365.25];
v5=Cos[2Pi Daysno/365.25];
aelRainModel=LinearModelFit[{Transpose[{u,v1,v2,v3,v4,v5}],Rainno}]
aelRainModel["ParameterTable"]
aelRainModel["AdjustedRSquared"]

(* ael: these are the estimates: *)
modelPoints=Table[aelRainModel[1,v1[[i]],v2[[i]],v3[[i]],v4[[i]],v5[[i]]],{i,1,Length[v1]}];
Length[modelPoints]
(* ael: this function is the trend, ignoring time *)
modelPointFunction[min_,max_]:=aelRainModel[1,max-min,(max-min)^2,0,0,0];

ael Now we’re going to plot everything, and see how it looks. The colored points are data; drab points are modeled. And the surface is the trend:

tempDiffModelPlot=Plot3D[modelPointFunction[x,y],{x,Min[Mintempno],Max[Mintempno]},{y,Min[Maxtempno],Max[Maxtempno]},PlotRange->{{Min[Mintempno],Max[Mintempno]},{Min[Maxtempno],Max[Maxtempno]}},PlotPoints->50,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","Rain Model"}];
modelPointPlot=ListPointPlot3D[Transpose[{Mintempno,Maxtempno,modelPoints}],PlotRange->All,AxesLabel->{"Min","Max","Rainfall"},PlotStyle->PointSize[0.03]];

Show[dataPlot,tempDiffModelPlot,modelPointPlot]

Now let’s do it again, only using a general quadratic of the two temperatures, rather than a quadratic in their difference:

v=Maxtempno;
w=Mintempno;
vs=v^2;
ws=w^2;
vw=v*w;
v3=Daysno/365.25;
v4=Sin[2Pi v3];
v5=Cos[2Pi v3];
aelRainModelGeneral=LinearModelFit[{Transpose[{u,v,w,vs,ws,vw,v3,v4,v5}],Rainno}]
aelRainModelGeneral["ParameterTable"]
aelRainModelGeneral["AdjustedRSquared"]

(* ael: these are the estimates: *)
modelPointGeneral=Table[aelRainModelGeneral[1,v[[i]],w[[i]],vs[[i]],ws[[i]],vw[[i]],v3[[i]],v4[[i]],v5[[i]]],{i,1,Length[v1]}];
(* ael: this function is the trend, ignoring time *)
modelPointGeneralFunction[min_,max_]:=aelRainModelGeneral[1,max,min,max^2,min^2,min*max,0,0,0];

tempDiffModelPlot=Plot3D[modelPointGeneralFunction[x,y],{x,Min[Mintempno],Max[Mintempno]},{y,Min[Maxtempno],Max[Maxtempno]},PlotRange->{{Min[Mintempno],Max[Mintempno]},{Min[Maxtempno],Max[Maxtempno]}},PlotPoints->50,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","Rain Model"}];
modelPointGeneralPlot=ListPointPlot3D[Transpose[{Mintempno,Maxtempno,modelPointGeneral}],PlotRange->All,AxesLabel->{"Min","Max","Rainfall"},PlotStyle->PointSize[0.03]];

Show[dataPlot,tempDiffModelPlot,modelPointGeneralPlot]

So it looks like the model with simply Max-Min is better: simpler, and almost as high an RSquared. 

Why is this model useful? 

	1. It’s simple.
	2. It incorporates temperature in a sensible way (as a function of the variation in temperature, measured by the difference between maxes and mins).
	3.  It incorporates time in a sensible way (as sinusoidal -- yearly -- variation).
	4. Once we derive models for temperatures Max(t) and Min(t), we have a model for rainfall (and we can estimate the impact increasing temperatures will have on rainfall). 
	
This is a very good step. Thanks much to Sally and Connor for this approach.

One suggestion, however: while the choice of Max-Min is sensible, as a measure of variability, why not consider
Max - slope Min
That is, some other slope than 1:

Manipulate[
Module[
{v1,v2,mod},
v1=Maxtempno-slope*Mintempno;
v2=v1^2;
mod=LinearModelFit[{Transpose[{u,v1,v2,v4,v5}],Rainno}];
{mod["ParameterTable"],
mod["AdjustedRSquared"]}
],
{{slope,1},.1,2,.01}
]

slopes=Table[slope,{slope,0,2,.01}];
rsquared=Table[
Module[{v1=Maxtempno-slope*Mintempno,v2,model},
v2=v1^2;
model=LinearModelFit[{Transpose[{u,v1,v2,v4,v5}],Rainno}];
{slope,model["AdjustedRSquared"]}
],{slope,0,2,.01}];
ListPlot[rsquared]

ael: In this case, it seems that a choice of Max-.5 Min gives rise to a better R-squared. It corresponds graphically to a slightly different parabolic trough. 

Now I’ll create the final model, which I’ll call the JaCS rainfall model (after Jacob E., who first indicated that differences in Max and Min temperatures are interesting; Connor, and Sally):

u=ConstantArray[1,Length[Maxtempno]];
linearMinMaxChoice[min_,max_]:=max-0.5 min
v1=linearMinMaxChoice[Mintempno,Maxtempno];
v2=v1^2;
aelRainModelFinal=LinearModelFit[{Transpose[{u,v1,v2,v4,v5}],Rainno}]
aelRainModelFinal["ParameterTable"]
aelRainModelFinal["AdjustedRSquared"]

(* ael: these are the estimates: *)
modelPointFinal=Table[aelRainModelFinal[1,v1[[i]],v2[[i]],v4[[i]],v5[[i]]],{i,1,Length[v1]}];
(* ael: this function is the trend, ignoring time *)
modelPointFunctionFinal[min_,max_]:=aelRainModelFinal[1,linearMinMaxChoice[min,max],(linearMinMaxChoice[min,max])^2,0,0];

ael Now we’re going to plot everything, and see how it looks. The colored points are data; drab points are modeled. And the surface is the trend:

finalModelPlot=Plot3D[modelPointFunctionFinal[x,y],{x,Min[Mintempno],Max[Mintempno]},{y,Min[Maxtempno],Max[Maxtempno]},PlotRange->{{Min[Mintempno],Max[Mintempno]},{Min[Maxtempno],Max[Maxtempno]}},PlotPoints->50,ColorFunction->"TemperatureMap",AxesLabel->{"Min","Max","Rain Model"}];
modelPointFinalPlot=ListPointPlot3D[Transpose[{Mintempno,Maxtempno,modelPointFinal}],PlotRange->All,AxesLabel->{"Min","Max","Rainfall"},PlotStyle->PointSize[0.03]];

Show[dataPlot,finalModelPlot,modelPointFinalPlot]

Now let’s see how well our model does versus the actual data as we rove over the data set:

maxRain=Max[Rainno];
scalerain=Rainno*.03/maxRain;
finalRainModel[min_,max_,day_]:=aelRainModelFinal[1,linearMinMaxChoice[min,max],(linearMinMaxChoice[min,max])^2,Sin[2Pi day/365.25],Cos[2Pi day/365.25]];
man=Manipulate[
p1=ListPlot[{{Mintempno[[t]],Maxtempno[[t]]}},PlotRange->All,PlotStyle->{Black,PointSize[0.02+scalerain[[t]]]}];
p2=ListPlot[{{minmodel[Daysno[[t]]],maxmodel[Daysno[[t]]]}},PlotRange->All,PlotStyle->{Blue,PointSize[0.02+Max[{0,finalRainModel[minmodel[Daysno[[t]]],maxmodel[Daysno[[t]]],Daysno[[t]]]}]*.03/maxRain]}];
Show[contour,p1,p2],
{{t,1},1,Length[Rainno],1}
]

ManToGif[man_,name_String,step_Integer: 1]:=
Module[
{moviename},
moviename=name<>If[$OperatingSystem==="MacOSX",".mov",".avi"];
Export[name<>".gif",
Import[Export[moviename,man],"ImageList"][[1;;-1;;step]]];
DeleteFile[moviename];
]

(* ManToGif[man,"SotoubouaJaCSModel",1] *)

And now we can predict the future....

Manipulate[
Module[
{time=(year-1960)*365.25,
min,max
},
min=minmodel[time];
max=maxmodel[time];
p2=ListPlot[{{min,max}},PlotRange->All,PlotStyle->{Blue,PointSize[0.02+Max[{0,finalRainModel[min,max,time]}]*.03/maxRain]}];
Show[contour,p2]
],
{{year,2015},1960,2100,1.0/12}
]

man=Manipulate[
Module[
{time=(year-1960)*365.25,
min,max
},
min=minmodel[time];
max=maxmodel[time];
p2=ListPlot[{{min,max}},PlotRange->All,PlotStyle->{Blue,PointSize[0.02+Max[{0,finalRainModel[min,max,time]}]*.03/maxRain]}];
Show[contour,p2]
],
{{year,2014},2014,2016,1.0/1200}
]

(* ManToGif[man,"SotoubouaJaCSModel2",1] *)

According to our model, there appears to be little change in modelled rainfall at the end of the century:

endOfCenturyRainfall=Table[
Module[
{time=(year-1960)*365.25,
min,max
},
min=minmodel[time];
max=maxmodel[time];
{year-2099,finalRainModel[min,max,time]}
],{year,2099,2100,.05}
];
endOfStudyRainfall=Table[
Module[
{time=(year-1960)*365.25,
min,max
},
min=minmodel[time];
max=maxmodel[time];
{year-2014,finalRainModel[min,max,time]}
],{year,2014,2015,.05}
];
pnow=ListPlot[endOfStudyRainfall]
pthen=ListPlot[endOfCenturyRainfall]
Show[pnow,pthen]

nn=7.5
time1=(2015-1960)*365.25+30*nn;
min1=minmodel[(2015-1960)*365.25];
max1=maxmodel[(2015-1960)*365.25];
finalRainModel[min1,max1,time1]
time2=(2099-1960)*365.25+30*nn;
min2=minmodel[(2099-1960)*365.25];
max2=maxmodel[(2099-1960)*365.25];
finalRainModel[min2,max2,time2]