Custar Data (near BG), from 1982 to 2020 -- 38 years.

Andy Long
Spring, 2020 -- year of the Covid-19 virus.

Let’s read in the data, which was created with CusterData-Create.nb:

iris=Import["http://ceadserv1.nku.edu/longa//classes/mat375/days/docs/Fletcher/data/Custar.csv","CSV"];
header=iris[[1]]
Custar=iris[[2;;]];
Custar[[1]]
Dimensions[Custar]

Now we’ll split up that matrix Custar into splitData, and then we’ll start peeling off certain variables, e.g. maxes:

splitData=Transpose[Custar];
maxes=splitData[[3]];
RecordMaxTemp=splitData[[4]];
RecordMaxYear=splitData[[5]];
ListPlot[maxes,PlotLabel->header[[3]]];

Big Question: How do we handle time, and dates? 

There are actually 365.25 days in a year, and then we wrap. So if we do a regression against index, my feeling is that we should divide the index by 365.25 and call it years, with maybe a shift of half a day to center it on the day:

day=Table[i,{i,0,Length[maxes]-1}];
year=day/365.25;
dayOfYear=0.5*1/365.25 + Table[day[[i]]-Floor[day[[i]]/365.25]*365.25,{i,1,Length[day]}];
ListPlot[dayOfYear,PlotLabel->header[[3]]]

and every four years we’re back on track: the date, year, and dayOfYear agree, or cycle:

date=splitData[[1]];
Length[date]
date[[1]]
date[[365*4+1 +1]]
dayOfYear[[1]]
dayOfYear[[365*4+1 +1]]
year[[1]]
year[[365*4+1 +1]]

Let’s also pull mins, humidity, solar radiation, and wind speed:

mins=splitData[[6]];
RecordMinTemp=splitData[[7]];
RecordMinYear=splitData[[8]];
avgRelativeHumidity=splitData[[14]];
avgSolarRadiation=splitData[[19]];
avgWindSpeed=splitData[[21]];

ListPlot[Transpose[{dayOfYear,RecordMinTemp}],PlotLabel->header[[7]]]
ListPlot[Transpose[{dayOfYear,RecordMaxTemp}],PlotLabel->header[[4]]]

GraphicsGrid[
{
{ListPlot[mins,PlotLabel->header[[6]]],ListPlot[maxes,PlotLabel->header[[3]]]},
{ListPlot[Transpose[{dayOfYear,RecordMinTemp}],PlotLabel->header[[7]]],ListPlot[Transpose[{dayOfYear,RecordMaxTemp}],PlotLabel->header[[4]]]},
{ListPlot[avgRelativeHumidity,PlotLabel->header[[14]]],ListPlot[avgSolarRadiation,PlotLabel->header[[19]]]}
}
]

dtr=maxes-mins;
xs=Table[i,{i,0,Length[dtr]-1}];
xs=xs/365.25;

Clear[b,m];
data=Transpose[{xs,dtr}];
dtrplot=ListPlot[Transpose[{xs,dtr}],PlotLabel->"DTR - Diurnal Temperature Range"]
(* There are some missing values, hence this was suggested by
Thanks to: https://mathematica.stackexchange.com/questions/133355/linearmodelfit-and-missing-data
*)
dtrlm=LinearModelFit[Cases[data,{_?NumericQ..}],{x},x]
dtrfits=Table[dtrlm[x],{x,xs}];
dtrlm["ParameterTable"]
dtrlm["RSquared"]
dtrlm["AdjustedRSquared"]
dtrlm["ANOVATable"]
dtrlm["ParameterConfidenceIntervals"]
Show[dtrplot,Plot[dtrlm[x],{x,0,xs[[Length[xs]]]}]]

Clear[b,m];
data=Transpose[{xs,maxes}];
maxplot=ListPlot[data,PlotLabel->"Maximum -- Custar Data"]
maxlm=LinearModelFit[Cases[data,{_?NumericQ..}],{x, Cos[2 Pi x],Sin[2Pi x]},x]
maxlm["ParameterTable"]
maxlm["RSquared"]
maxlm["AdjustedRSquared"]
maxlm["ParameterConfidenceIntervals"]
data=Transpose[{xs,dtr}];
Show[maxplot,Plot[maxlm[x],{x,0,xs[[Length[xs]]]}]]
Histogram[maxlm["FitResiduals"]]

Clear[b,m];
data=Transpose[{xs,mins}];
minplot=ListPlot[data,PlotLabel->"Minimum -- Custar Data"]
minlm=LinearModelFit[Cases[data,{_?NumericQ..}],{x, Cos[2 Pi x],Sin[2Pi x]},x]
minfits=Table[minlm[x],{x,xs}];
minlm["ParameterTable"]
minlm["RSquared"]
minlm["AdjustedRSquared"]
minlm["ParameterConfidenceIntervals"]
data=Transpose[{xs,dtr}];
Show[minplot,Plot[minlm[x],{x,0,xs[[Length[xs]]]}]]
Histogram[minlm["FitResiduals"]]

Manipulate[
Show[Plot[minlm[x],{x,start,start+nyrs}],minplot],
{{start,0},0,37},
{{nyrs,1},0,5}
]

m=Length[mins]
bins={-30,30,3}
minDiffs=Table[mins[[i]]-mins[[i+1]],{i,1,m-1}];
Histogram[minDiffs,bins]
n=Length[maxes]
maxDiffs=Table[maxes[[i]]-maxes[[i+1]],{i,1,n-1}];
Histogram[maxDiffs,bins]

data=Transpose[{mins[[2;;m]],minDiffs}];
mindiffplot=ListPlot[data]
data=Cases[data,{_?NumericQ..}];
mindifflm=LinearModelFit[data,{x},x]
mindifflm["ParameterTable"]
mindifflm["RSquared"]
mindifflm["AdjustedRSquared"]
mindifflm["ParameterConfidenceIntervals"]
Show[mindiffplot,Plot[mindifflm[x],{x,-20,80}]]
Histogram[mindifflm["FitResiduals"]]

data=Transpose[{maxes[[2;;n]],maxDiffs}];
maxdiffplot=ListPlot[data]
data=Cases[data,{_?NumericQ..}];
maxdifflm=LinearModelFit[data,{x},x]
maxdifflm["ParameterTable"]
maxdifflm["RSquared"]
maxdifflm["AdjustedRSquared"]
maxdifflm["ParameterConfidenceIntervals"]
Show[maxdiffplot,Plot[maxdifflm[x],{x,-10,110}]]
Histogram[maxdifflm["FitResiduals"]]

I want to estimate the distribution of “next day’s temperature extreme” as a function of the current temperature. 
The idea is that, if you’re near an extreme, you’re much more likely to “regress to the mean”. So at each temp T, you’ll as P(tomorrow temp=t|]

data=Transpose[{maxes[[1;;n-1]],maxDiffs}];
data=Cases[data,{_?NumericQ..}];
low=Min[Transpose[data][[1]]];
high=Max[Transpose[data][[1]]];
mean=Mean[Transpose[data][[1]]];
Manipulate[
Histogram[
Transpose[
Select[data,temp-degrees≤ #[[1]]≤temp+degrees &]
][[2]],
{-30,30,5},
PlotLabel->temp
],
{{temp,mean},low,high},
{{degrees,5},1,10}
]

data=Transpose[{mins[[1;;n-1]],minDiffs}];
data=Cases[data,{_?NumericQ..}];
low=Min[Transpose[data][[1]]];
high=Max[Transpose[data][[1]]];
mean=Mean[Transpose[data][[1]]];
Manipulate[
Histogram[
Transpose[
Select[data,temp-degrees≤ #[[1]]≤temp+degrees &]
][[2]],
{-40,40,5},
PlotLabel->temp
],
{{temp,mean},low,high},
{{degrees,5},1,10}
]

data=Transpose[{mins[[1;;n-1]],minDiffs}];
data=Cases[data,{_?NumericQ..}];
low=Min[Transpose[data][[1]]];
high=Max[Transpose[data][[1]]];
mean=Mean[Transpose[data][[1]]];
Manipulate[
vals=Transpose[
Select[data,temp-degrees≤ #[[1]]≤temp+degrees &]
][[2]];
edist=EstimatedDistribution[vals,If[temp>mean,GumbelDistribution[α,β],
ExtremeValueDistribution[α,β]]];
h=DistributionFitTest[vals,edist,"HypothesisTestData"];
Show[
Histogram[vals,{-40,40,5},"ProbabilityDensity",PlotLabel->temp],
Plot[PDF[h["FittedDistribution"],x],{x,-40,40 },PlotStyle->Thick]
],
{{temp,mean},low,high},
{{degrees,5},1,10}
]

ExtremeValueDistribution
The extreme value distribution gives the asymptotic distribution of the maximum value in a sample from a distribution such as the normal distribution.

The asymptotic distribution of the minimum value, also sometimes called an extreme value distribution, is implemented in the Wolfram Language as (RefLink[GumbelDistribution, paclet:ref/GumbelDistribution]). »

GumbelDistribution

data=Transpose[{maxes[[1;;n-1]],maxDiffs}];
data=Cases[data,{_?NumericQ..}];
low=Min[Transpose[data][[1]]];
high=Max[Transpose[data][[1]]];
mean=Mean[Transpose[data][[1]]];
mean=Mean[{high,low}];
Manipulate[
vals=Transpose[
Select[data,temp-degrees≤ #[[1]]≤temp+degrees &]
][[2]];
(*
ℋ=DistributionFitTest[vals,Automatic,"HypothesisTestData"];
Print[ℋ["TestDataTable",All]];
*)
eGumbel=EstimatedDistribution[vals,GumbelDistribution[α,β]];
eExtreme=EstimatedDistribution[vals,ExtremeValueDistribution[α,β]];
hGumbel=DistributionFitTest[vals,eGumbel,"HypothesisTestData"];
hExtreme=DistributionFitTest[vals,eExtreme,"HypothesisTestData"];
h=If[hGumbel["PValue"] > hExtreme["PValue"] ,hGumbel,hExtreme];
e=If[hGumbel["PValue"] > hExtreme["PValue"] ,eGumbel,eExtreme];
Print[e];
Show[
Histogram[vals,{-40,40,5},"ProbabilityDensity",PlotLabel->h["PValue"] ],
Plot[PDF[h["FittedDistribution"],x],{x,-40,40 },PlotStyle->Thick]
],
{{temp,mean},low,high},
{{degrees,5},1,10}
]

data=Transpose[{mins[[1;;n-1]],minDiffs,mins[[1;;n-1]]-minfits[[1;;n-1]]}];
Dimensions[data]
data=Cases[data,{_?NumericQ..}];
Dimensions[data]
splitData=Transpose[data];
low=Min[splitData[[3]]];
high=Max[splitData[[3]]];
mean=Mean[splitData[[3]]];
Manipulate[
Histogram[
Transpose[
Select[data,temp-degrees≤ #[[3]]≤temp+degrees &]
][[2]],
{-40,40,5},
PlotLabel->temp
],
{{temp,mean},low,high},
{{degrees,2},0.5,4}
]

data=Transpose[{maxes[[1;;n-1]],maxDiffs,maxes[[1;;n-1]]-maxfits[[1;;n-1]]}];
Dimensions[data]
data=Cases[data,{_?NumericQ..}];
Dimensions[data]
splitData=Transpose[data];
low=Min[splitData[[3]]];
high=Max[splitData[[3]]];
mean=Mean[splitData[[3]]];
Manipulate[
Histogram[
Transpose[
Select[data,temp-degrees≤ #[[3]]≤temp+degrees &]
][[2]],
{-40,40,5},
PlotLabel->temp
],
{{temp,mean},low,high},
{{degrees,2},0.5,4}
]

Clear[b,m];
data=Transpose[{xs,avgSolarRadiation}];
humidplot=ListPlot[Transpose[{xs,avgSolarRadiation}],PlotLabel->"avgSolarRadiation"]

(* There are some missing values, hence this was suggested by
Thanks to: https://mathematica.stackexchange.com/questions/133355/linearmodelfit-and-missing-data
*)
humlm=LinearModelFit[Cases[data,{_?NumericQ..}],{x,x^2},x]
humfits=Table[humlm[x],{x,xs}];
humlm["ParameterTable"]
humlm["RSquared"]
humlm["AdjustedRSquared"]
humlm["ParameterConfidenceIntervals"]
estimates=humlm["BestFitParameters"]
1982+estimates[[2]]/(-2 * estimates[[3]])

Show[humidplot,Plot[humlm[x],{x,0,xs[[Length[xs]]]}]]