Togo Models for Temperatures and Rainfall with ENSO

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];

Now I 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}];
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}];

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,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)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)],
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
enso
},{lat,long,elev,enso,x}]
finalmin["ParameterTable"]
finalmin["AdjustedRSquared"]

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

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

Global Togo Model for Maximum Temperature

We generate the Maximum Temperature Model:

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

maxdata=Transpose[{latitude,longitude,elevation,ENSOvec,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/(1/5)],Cos[2*π*x/(1/5)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
enso
},{lat,long,elev,enso,x}]
finalmax["ParameterTable"]
finalmax["AdjustedRSquared"]

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

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

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,rainfall}],{x},{x}]
finalrain["ParameterTable"]
finalrain["RSquared"]

raindata=Transpose[{decYear,latitude,longitude,elevation,minTemp,maxTemp,ENSOvec,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)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
mint,
maxt,
enso
},{x,lat,long,elev,mint,maxt,enso}]

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)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
mint,
maxt,
enso
},{x,lat,long,elev,mint,maxt,enso},IncludeConstantBasis->False]
finalrain["ParameterTable"]
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"]]

The rainfall residuals are not normal :(                :

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