The Calculus of Climate Change
Andy Long, MAT128
Spring, 2023

Case Study: 

Washington’s Treasured Cherry Blossoms Prompt Reflection on Local Climate Change

Over the years, the trees have been blooming earlier and earlier. Does it matter? Yes and no, experts say.

By Emma Ricketts (Inside Climate News)
March 18, 2023

“Each year, the National Park Service calculates when “peak’’ bloom of the cherry trees—when 70 percent of the blossoms are open—is expected to occur. This year, it is anticipated between March 22 and March 25. While directly in line with last year’s peak, that’s a full two weeks earlier than the area’s historic average.”

March 22 is day 81 (82 on leap years -- that in itself is an interesting problem). So they are claiming that the historic average is sometime between day 95 and 98.

Three questions of interest to the mathematician:

	1. Is it true that the trees have been blooming earlier and earlier? 
	3. What is the “historic average”?
	3. If trees are blooming earlier, then what is a good model for the change in 
	    blooming schedule of the trees?

Import the data, from https://www.epa.gov/climate-indicators/cherry-blossoms:

(* data=Import["/Users/longa/Desktop/NKU/classes/2023Spring/mat128/climate/cherry-blossoms.csv"] *)

Now we arrange and plot the data:

(* It's such a small data set that it is easy enough to explicitly add it to this file: *)
data={{"Year","YoshinoPeakBloomDate"},{1921,79},{1922,97},{1923,99},{1924,104},{1925,86},{1926,101},{1927,79},{1928,99},{1929,90},{1930,91},{1931,101},{1932,106},{1933,99},{1934,105},{1935,91},{1936,98},{1937,104},{1938,84},{1939,89},{1940,104},{1941,102},{1942,95},{1943,94},{1944,100},{1945,79},{1946,82},{1947,102},{1948,88},{1949,88},{1950,99},{1951,96},{1952,100},{1953,86},{1954,96},{1955,92},{1956,97},{1957,98},{1958,108},{1959,96},{1960,105},{1961,92},{1962,97},{1963,93},{1964,102},{1965,105},{1966,95},{1967,96},{1968,90},{1969,99},{1970,106},{1971,98},{1972,102},{1973,101},{1974,93},{1975,93},{1976,83},{1977,85},{1978,102},{1979,92},{1980,97},{1981,93},{1982,97},{1983,97},{1984,94},{1985,97},{1986,92},{1987,87},{1988,91},{1989,88},{1990,74},{1991,88},{1992,98},{1993,101},{1994,95},{1995,92},{1996,95},{1997,85},{1998,86},{1999,95},{2000,77},{2001,96},{2002,92},{2003,92},{2004,91},{2005,99},{2006,89},{2007,91},{2008,86},{2009,91},{2010,90},{2011,88},{2012,80},{2013,99},{2014,100},{2015,100},{2016,85},{2017,84},{2018,95},{2019,91},{2020,80},{2021,87},{2022,80}};year=Rest[data[[All,1]]]-1970; (* 1970 was chosen as the reference year for reasons to be seen later *)
peak=Rest[data[[All,2]]];
data = Transpose[{year,peak}];
pl=ListPlot[data,PlotLabel->"Peak Day (from January 1)",AxesLabel->{"Year from 1970","Day"},ImageSize->600];
Show[pl,Plot[N[Mean[peak]],{x,-50,55},PlotStyle->Red]]
N[Mean[peak]]

The mean value of Peak Day over the data set is Day 93.5 -- generally around April 4rd/5th. That’s called the “null model” (no change); the question is whether we can reject the null model? Has there been change?

Analysis

Confidence is 95%, by default. We’ll go with the default. So statements made  reflect the notion that  we will reject a null in favor of an alternative if it would happen less than 5% of the time just by chance.

Let’s try a linear regression model, which would argue that the form of the peak day model is

			Peak(year) = Day1970 + slope (year - 1970)

lm=LinearModelFit[data,t,t]
Show[pl,Plot[lm[t],{t,-50,55}]]
lm["RSquared"]
lm["ParameterConfidenceIntervalTable"]
Histogram[lm["FitResiduals"]]

{bands90[t_],bands95[t_],bands99[t_],bands999[t_]}=Table[lm["MeanPredictionBands",ConfidenceLevel->cl],{cl,{.90,.95,.99,.999}}];
pun=Show[pl,Plot[{lm[x],bands90[x],bands95[x],bands99[x],bands999[x]},{x,-50,55},Filling->{2->{1},3->{2},4->{3},5->{4}},PlotLegends->{"model","90% confidence bands","95% confidence bands","99% confidence bands","99.9% confidence bands"}]];
cis=lm["ParameterConfidenceIntervals",ConfidenceLevel->.95]
Manipulate[
Show[pun,
Plot[m*x+b,{x,-50,55},PlotStyle->{Thick,Red}]],
{{b,Mean[cis[[1,All]]]},cis[[1,1]],cis[[1,2]]},
{{m,Mean[cis[[2,All]]]},cis[[2,1]],cis[[2,2]]}
]

We can reject the null hypothesis of no change, in favor of the hypothesis that there has been change over  time. The slope parameter is negative across the spectrum of the 95% confidence interval, effectively rejecting a possibility of slope 0 -- constant Peak Day. 

The histogram of residuals is quite skew -- whereas it should be normal. 

The 1970 “historical” value is estimated to be day 94.... but that would not be an historical average, since the slope is negative, and so it would have been higher in the past. Let’s think about how we might model that historical average.

Let’s try a quadratic model:

qm=LinearModelFit[data,{t,t^2},t]
Show[pl,Plot[qm[t],{t,-50,55}]]
qm["RSquared"]
qm["ParameterConfidenceIntervalTable"]
Histogram[qm["FitResiduals"]]
peakYear=t/.Solve[-0.06453265960958882` -2* 0.001951425428533258` t==0,t];
peakYear=peakYear[[1]];
{peakYear+1970, qm[peakYear]}

{bands90[t_],bands95[t_],bands99[t_],bands999[t_]}=Table[qm["MeanPredictionBands",ConfidenceLevel->cl],{cl,{.90,.95,.99,.999}}];
Show[pl,Plot[{qm[x],bands90[x],bands95[x],bands99[x],bands999[x]},{x,-50,55},Filling->{2->{1},3->{2},4->{3},5->{4}},PlotLegends->{"model","90% confidence bands","95% confidence bands","99% confidence bands","99.9% confidence bands"}]]

Again we can reject the null hypothesis of no change, in favor of the hypothesis that there has been change over  time. In particular the quadratic term among the parameters is negative across the spectrum of the 95% confidence interval, suggesting an “umbrella” -- and that there is an accelerating decrease in the Peak Day. 

The histogram of residuals is again rather skew -- whereas it should be normal. 

The peak of this model occurs around 1953, with a maximum of Day 96.

Let’s try a cutpoint model:

To get at the “historical average”, we might assume that things started off with a historical average (constant value), and then, at some point, a change began to kick in, with a likely decreasing Peak Day.

cm=NonlinearModelFit[data,If[t<cut,constant,constant+slope(t-cut)],{{cut,-10},constant,slope},t]
Show[pl,Plot[cm[t],{t,-50,55}]]
cm["ParameterConfidenceIntervalTable"]
Histogram[cm["FitResiduals"]]

{bands90[t_],bands95[t_],bands99[t_],bands999[t_]}=Table[cm["MeanPredictionBands",ConfidenceLevel->cl],{cl,{.90,.95,.99,.999}}];
Show[pl,Plot[{cm[x],bands90[x],bands95[x],bands99[x],bands999[x]},{x,-50,55},Filling->{2->{1},3->{2},4->{3},5->{4}},PlotLegends->{"model","90% confidence bands","95% confidence bands","99% confidence bands","99.9% confidence bands"}]]

This model suggests that the historical average was within the band of Days 93.5 and 97.5; then, as in the previous cases, there is a significantly negative change, starting around 1970 (actually 1971). 

The histogram of residuals is better. 

One thing that we don’t like about this model is that abrupt (slope discontinuity).

Let’s try a cutpoint model with a switch to a quadratic:

cqm=NonlinearModelFit[data,If[t<cut,constant,constant+quad(t-cut)^2],{{cut,0},constant,{quad,-.2}},t]
Show[pl,Plot[cqm[t],{t,-50,55}]]
cqm["ParameterConfidenceIntervalTable"]
Histogram[cqm["FitResiduals"]]

{bands90[t_],bands95[t_],bands99[t_],bands999[t_]}=Table[cqm["MeanPredictionBands",ConfidenceLevel->cl],{cl,{.90,.95,.99,.999}}];
Show[pl,Plot[{cqm[x],bands90[x],bands95[x],bands99[x],bands999[x]},{x,-50,55},Filling->{2->{1},3->{2},4->{3},5->{4}},PlotLegends->{"model","90% confidence bands","95% confidence bands","99% confidence bands","99.9% confidence bands"}]]

This model is interesting: it shows a historical  average of about Day 95.5. But it also doesn’t promise a significant quadratic term, meaning that its value might be 0 -- meaning a constant model.... This model would say that the historical average was in play until about 1955, at which point the downturn began; but the uncertainty in that cut-point parameter is very high -- the 95% confidence interval basically includes the whole time span of data collection.

A final word:

“The cherry trees also face the reality of sea level rise. Water levels in the Tidal Basin are approximately four feet higher than they were when it was built 80 years ago, Litterst said. “We’ve had to remove more than a dozen cherry trees because their roots simply can’t take the constant inundation of water.” 

The AP ran a similar story on March 1, 2023: DC's cherry blossoms coming early due to confusing weather. (By Ashraf Khalil)