Climate change can be defined as the variation in the joint probability distribution describing the state of the atmosphere, oceans, and fresh water including ice (Hsiang and Kopp, 2018). These are complex, multidimensional physical systems, and the various features of climate change have been described using a diverse set of summary statistics. One of the most important aspects of climate change is the evolving distribution of temperature, which has important implications for growth, urbanization, agriculture, health, productivity, and poverty (Kahn et al., 2019;Colacito et al., 2019;Carleton et al., 2020). A diverse set of indicators have been used to measure such distributional temperature variation, including, for example, mean temperature, temperature range, hot and cold spell duration, frost days, growing season length, ice days, heating and cooling degree days, and start of spring dates (Masson-Delmotte et al., 2018;Reidmiller et al., 2018). Of course, the level of temperature -the central tendency of the distribution -has attracted the most attention, in particular regarding the upward trend in the average daily temperature (AVG) (Rivas and Gonzalo, 2020). In contrast, less attention has been given to temperature volatility. 1Temperature volatility can be measured by the diurnal temperature range (DTR), which is the difference between the daily maximum temperature (MAX) and minimum temperature (MIN). Similar to changes in temperature averages, changes in temperature ranges and variability can also have important effects on environmental and human health (Davy et al., 2017). For example, the incidence of temperature extremes such as heat waves depends critically on how the whole distribution of temperature is shifting -including both central tendency and variability. Of course, such temperature extremes can have notable adverse effects on society and the economy. Temperature variability can stress workers and lower labor productivity, but it can also have direct effects on output. A salient example is agriculture, whose output is a function of capital, labor, and weather inputs.2 Indeed, the very viability of certain agricultural sub-industries, notably wine and coffee production, is crucially dependent on temperature ranges. For example, Robinson (2006) notes that Diurnal temperature variation is of particular importance in viticulture. Wine regions situated in areas of high altitude experience the most dramatic swing in temperature variation during the course of a day. In grapes, this variation has the effect of producing high acid and high sugar content as the grapes' exposure to sunlight increases the ripening qualities while the sudden drop in temperature at night preserves the balance of natural acids in the grape. (p. 691) To better understand the full nature of the changing distribution of temperature, we examine DTR in select cities in the United States over the past half-century. We allow for time-varying coefficients, which provide a straightforward yet flexible approximation to more general nonlinear effects. Although our focus is on DTR, we also provide a parallel analysis for AVG, which allows valuable interpretive context and contrast. Our work reveals an evolving DTR conditional mean seasonal pattern, in contrast to the fixed AVG conditional mean seasonal pattern. The previous research literature that examined DTR struggled for some time to develop firm conclusions about the dynamics of temperature variability. Even the direction of the trend in DTR has been somewhat contentious (Alexander and Perkins, 2013). Recent work has established that the downward trend in DTR in many locations (Sun et al., 2019) reflects a more rapid warming of MIN than MAX -generally the result of nighttime lows rising faster than daytime highs (Davy et al., 2017). However, this differential trending of MIN and MAX, or "diurnal asymmetry," is not geographically uniform because of variation in vegetation, cloud cover, and other factors (Jackson and Forster, 2010;Sun and Pinker, 2014). Along with this trend in temperature variability, seasonal variation in DTR has also been considered by a few authors who describe a lower temperature range in winter than at other times (Ruschy et al., 1991;Durre and Wallace, 2001). There is also some evidence that the seasonality of DTR in the United States may be changing over time (Qu et al., 2014). To capture as much variation as possible in the distribution of DTR -including trend and seasonal -we use linear time series models with time-varying coefficients to provide simple yet powerful representations. We proceed as follows. In section 2, we provide an introductory analysis for a representative city, Philadelphia. Then, in section 3, we broaden the analysis to include fifteen geographically dispersed U.S. cities. We conclude in section 4. We introduce and illustrate our approach by studying temperature data measured at the Philadelphia airport (PHL) in a step-by-step fashion. We present most results graphically, while regression results on which these graphs are based appear in Appendix A. 3 The underlying data are the daily MAX and MIN measured in degrees Fahrenheit, obtained from the U.S. National Ocean and Atmospheric Administration's Global Historical Climate Network database (GHCN-daily).4 Our sample period is from 01/01/1960 to 12/31/2017, which covers the period of almost all recent climate change. The daily MAX and MIN are informative of both the central tendency and variability of the daily continuous-time temperature record. In particular, the daily average temperature, AVG=(MAX+MIN)/2, is a natural measure of central tendency, and the daily temperature range, DTR=MAX-MIN, is a natural measure of volatility or variability.5 DTR is a natural and intuitive estimator of daily volatility, and it is also highly efficient statistically. The "daily range" has a long and distinguished tradition of use in econometrics due to its good properties in estimating underlying quadratic variation from discretely-sampled data (Alizadeh et al., 2002). Interestingly, although AVG has been studied extensively (Raftery et al., 2017), DTR has been studied much less. In Figure 1, we show kernel estimates of the unconditional densities of AVG and DTR. The bimodal shape of the AVG density reflects the strong seasonality in AVG. The "winter mode" is around 40 • F, and the "summer mode" is around 75 • F. The AVG density contrasts sharply with the unimodal approximately-symmetric density of DTR, which is centered around 19 • F and much less dispersed.
In Figure 2, we display time series plots of the entire data sample of AVG and DTR with fitted linear trends superimposed. The regression is

where Y is AVG or DTR, c is a constant, and T IM E is a time trend (that is, T IM E t = t and t = 1, ..., T ). Here and throughout, we use heteroskedasticity and autocorrelation consistent (HAC) standard errors to assess statistical significance (Newey and West, 1987). The AVG trend slopes upward and is statistically significant, which is consistent with the overall global warming during this period. The steepness of this trend is surprising, as the AVG trend grows by nearly five degrees Fahrenheit over the course of the 57-year 1960-2017 sample. This increment is a bit more than twice as much as the average global increase over the same period (Rudebusch, 2019). The faster upward trend in the Philadelphia airport average temperature likely reflects two key factors: (1) average temperatures in growing cities tend to rise more quickly due to an increasing urban heat island effect and (2) average land temperatures generally grow more quickly than the global average, which includes ocean areas that are slow to warm. As for Philadelphia temperature variability, DTR also has a significant trend, and it slopes downward, dropping by more than two degrees over the course of the sample. The downward DTR trend arises from different trends in the underlying MAX and MIN -a diurnal asymmetry. Both trend upward, but MIN is on a steeper incline as evening temperatures warm more quickly. Hence, the spread between MAX and MIN tends to shrink, and DTR decreases over time. The relatively muted upward trend in MAX is generally ascribed to increased cloud cover, soil moisture, and precipitation, which lead to diminished incoming solar radiation and increased daytime surface evaporative cooling; however, with local variation in these meterological elements, a downward DTR trend is not found at all locations (Dai et al., 1999;Davy et al., 2017;Vinnarasi et al., 2017). The overall picture, then, involves not only an upward trend in AVG, but also a gradual tightening of daily fluctuations around that trend. Warming is not only happening, but also happening more reliably. In Figure 3, we show the actual and fitted values from regressions of de-trended AVG and DTR on 12 monthly seasonal dummies,

where Y is de-trended AVG or DTR -the residuals from regression (1) -and D it = 1 if day t is in month i, and 0 otherwise.6 This model is effectively an intercept regression for deviations from trend, allowing for a different intercept each month. As shown in the top panel of Figure 3, AVG displays pronounced seasonality. The seasonality is highly significant and is responsible for a large amount AVG variation. The R 2 of the seasonal AVG regression (2) is .81. As with the upward trend in AVG, strong seasonality in deviations of AVG from its trend is hardly surprising -it's cold in the winter and hot in the summer. There is also significant seasonality in DTR, as shown in the bottom panel of Figure 3. The DTR seasonality was hard to detect visually in the time series plot of Figure 2, because it is buried in much more noise than that of AVG. The R 2 of the seasonal DTR regression (2) is only .07. In Figure 4, we show the estimated monthly seasonal factors for AVG (left panel) and DTR (right panel). They are simply the 12 estimated coefficients on the 12 monthly dummies in the seasonal regression (2). The seasonal pattern for AVG is as expected -smooth and unimodal, high in the summer and low in the winter, achieving its maximum in July and its minimum in January. In contrast, the seasonal pattern for DTR is clearly bi-modal, with one mode in April-May and one in October. DTR's two annual peaks (spring and fall) and two annual troughs (winter and summer) contrast sharply with AVG's single annual peak (summer) and single annual trough (winter). This "twin-peaks" or "M-shaped" DTR pattern is common across many U.S. cites. Moreover, as we shall show, in many locations, the DTR seasonal pattern has evolved noticeably over time with climate change, whereas the AVG seasonal pattern has remained stable. The AVG and DTR trends documented thus far are trends in level. More subtle are trends in seasonality -that is, trends in the tent-shaped AVG seasonal pattern and the M-shaped DTR seasonal pattern. We now explore the possibility of such evolving seasonality by allowing for linear trends in the seasonal factors. 7 Mechanically, this involves regressing de-trended AVG or DTR not only on 12 monthly dummies, but also those same 12 dummies interacted with time,

where Y is de-trended AVG or DTR, D it = 1 if day t is in month i and 0 otherwise, and T IM E t = t. Regression (3) can capture linearly-trending seasonal deviations from a linear trend. Effectively, it allows for a different intercept each month, with those intercepts themselves potentially trending at different rates. In the special case where all interaction coefficients are zero, it collapses to fixed seasonal deviations from linear trend, as explored in section 2.3. For AVG, there are no gains from estimating the more flexible seasonal specification (3). The interaction terms are universally insignificantly different from zero, clearly indicating no change over time in the AVG seasonal pattern. In the left panel of Figure 5, we show the estimated seasonal factors for AVG for the first year (1960) and last year (2017) of our sample. This range provides the maximum contrast, but the two seasonal patterns are nevertheless essentially identical. The results for DTR, however, are very different. Unlike the AVG seasonal, which does not evolve, the DTR seasonal changes significantly over time. The January-through-March DTR interaction coefficients are significantly positive, indicating that the winter DTR low is increasing. In addition, all May-through-October interaction coefficients are negative, and the October coefficient is large and highly significantly negative.
This corresponds to progressively lower DTR highs in Octobers, so that the fall DTR peak is gradually vanishing. Both effects (higher winter DTR lows, and lower fall DTR highs) are visually apparent in the right panel of Figure 5, in which we contrast the estimated DTR M-shaped seasonal pattern in the first year (1960) and last year (2017) of our sample. We now expand our analysis to include data from the airports of the fifteen U.S. cities shown in Figure 6. 8 As with the Philadelphia case study in section 2, we obtain the underlying daily MAX and MIN data, from which we construct daily AVG and DTR, from the U.S. National Ocean and Atmospheric Administration's GHCN-daily, https://www.ncdc.noaa. gov/ghcn-daily-description. Our sample period is 01/01/1960-12/31/2017.9

We choose these city weather reporting stations because all of them have had temperature derivatives traded on the Chicago Merchantile Exchange (CME) (Campbell and Diebold, 2005). Consideration of such CME cities is of interest for several reasons. First, these locations cover a diverse set of climates, so they can provide a check of the robustness of our Philadelphia results. Second, they are urban locations that represent large numbers of people and a sizable share of economic activity -one reason that their CME contracts are traded. Finally, the valuations of weather derivatives traded in financial markets depend on the evolution of the stochastic structure of temperature dynamics, which is precisely the focus of our modeling efforts and so naturally paired with the CME cities. The full set of historically-traded cities includes: Atlanta, ATL; Boston, BOS; Baltimore Washington, BWI; Chicago, ORD; Cincinnati, CVG; Dallas Fort Worth, DFW; Des Moines, DSM; Detroit, DTW; Houston, IAH; Kansas City, MCI; Las Vegas, LAS; Minneapolis St Paul, MSP; New York, LGA; Portland, PDX; Philadelphia, PHL; Sacramento, SAC; Salt Lake City, SLC, and Tuscon, TUS. We exclude Houston, Kansas City, and Sacramento, however, due to large amounts of missing data, leaving fifteen cities. Presently eight cities are traded (Atlanta, Chicago, Cincinnati, Dallas, Las Vegas, Minneapolis, New York, and Sacramento), and all but Sacramento are in our fifteen. 10  We focus on urban airports not only because they are focal points for the financial weather derivatives community, but also because they are focal points for the meteorological / climatological measurement community. The U.S. National Weather Service surface weather observation network, in particular, is located largely at airports. On balance, airports, with large open spaces of grassland without surrounding buildings, are regarded as fairly representative locations for measuring temperature. We view the sequential modeling approach employed in section 2 -fitting a trend and then characterizing seasonality in the de-trended data -as intuitive and transparent. We now consolidate and extend various aspects of that approach, to arrive at a simple yet powerful joint model. Regarding consolidation, we move from a multi-step sequential modeling approach to a single-step joint approach with a single conditional mean estimation. Regarding extension, we now include an autoregressive lag in the model, which facilitates simple assessment of the strength of serial correlation in the deviations from the trend/seasonal. The autoregressive lag also provides valuable pre-whitening for serial correlation in HAC covariance matrix estimation (Andrews and Monahan, 1992). We proceed by regressing AVG or DTR on an intercept, a linear trend, a first-order autoregressive lag, 11 monthly seasonal dummies to capture seasonal intercept variation (we drop July, so the included constant captures July and all estimated seasonal effects are relative to July), and 11 trend-seasonal interactions to capture seasonal trend slope variation (we drop the July interaction):

10 See https://www.cmegroup.com/trading/weather/temperature-based-indexes.html. where Y is AVG or DTR, T IM E t = t, Y (-1) denotes a 1-day lag, and D it = 1 if day t is in month i and 0 otherwise. The joint model (4) allows for different intercepts each month, with the different intercepts potentially trending linearly at different rates, and for serially correlated deviations from the trend/seasonal. We summarize the estimation results in Tables 1 and2, in which we show the weather station identifier (airport code) in column 1, and various aspects of the estimation results in subsequent columns.11

As shown in column 2 of Table 1, the estimated AVG trend movements over the full sample are large and positive in each city. They are also all highly statistically significant (column 3), with a median p-value of 0.00 for Wald tests of the null hypothesis of no trend. These p-values are denoted p(nt), where "nt" stands for "no trend", which corresponds to zero coefficients on TIME and all TIME interactions in regression (4) (in which case it collapses to seasonal intercepts with serial correlation). The median estimated trend movement is 3.38 • F, greater than the consensus estimate of the increase in the mean global temperature over the same period, as U.S. airports have warmed more quickly than the global average. Similarly, in column 2 of Table 2, we report the estimated full-sample trend movements for DTR. All but one are negative, and most are significant at the one percent level. The median estimated trend movement is -1.45 • F, with a median p-value, p(nt), of 0.00 for the no-trend null hypothesis (column 3). Interestingly, LAS, which has the steepest upward AVG trend, also has the steepest downward DTR trend. In column 4 of Tables 1 and2, we report p-values for Wald tests of the hypothesis of no AVG and DTR seasonality, respectively. These p-values are denoted p(ns), where "ns" stands for "no seasonality", which corresponds to zero coefficients on all included seasonal dummies and dummy interactions in regression (4) (in which case it collapses to linear trend with serial correlation). There is of course strong evidence of seasonality in AVG with all p(ns)'s equal to 0.00. Less well known is the similarly strong seasonality in DTR with all p(ns)'s again equal to 0.00. In column 5 of Tables 1 and2, we report p-values for Wald tests of the hypothesis of no evolving (i.e., trending) AVG and DTR seasonality, respectively.
These p-values are denoted p(nts), where "nts" stands for "no trending seasonality", which corresponds to zero coefficients on all seasonal dummy interactions in regression (4) (in which case it collapses to linear trend and fixed seasonal dummies with serial correlation). The results are striking. There is no evidence for evolving seasonality in AVG; the median AVG p(nts) is 0.73. In contrast, there is strong evidence of evolving seasonality in DTR; the median DTR p(nts) is 0.03. Estimated AVG and DTR serial correlation coefficients appear in column 6 of Tables 1 and2, respectively. 12 All are positive and significant at the one percent level. Their magnitudes, however, are very different. All those for AVG are around 0.75, whereas all those for DTR are around 0.35. It is interesting to note that, although the signals in both AVG and DTR are clearly driven by trend, seasonal, and cyclical components, the DTR signal is burried in much more noise. As shown in column 7 of Tables 1 and2, respectively, all AVG regression R 2 values are around 0.9, whereas all those for DTR are around 0.2. 12 One could also allow for the possibility of more sophisticated forms of serial correlation. There is, for example, some borderline evidence of long memory in the AVG equation ( 4), as discussed in Campbell and Diebold (2005). They did not, however, study DTR, and to the best of our knowledge long memory has not been explored in the DTR equation (4). Some preliminary exploration allowing for fractionally-integrated ARFIMA(p,d,q) disturbances produced mixed results, but overall there was some evidence of long memory (0<d<0.5). We leave a full exploration to future work, particularly as there are many issues and nuances regarding long memory and its relationship to structural change, as emphasized by Diebold and Inoue (2001). Climate change is one of the most consequential and pressing issues of our time. We have focused on DTR as an important summary statistic for characterizing climate change, and we have estimated new stochastic representations of DTR that can capture its changing seasonality over time and space. Throughout we have also provided parallel contrasting results for AVG. Many extensions are possible, such as full bivariate modeling of AVG and DTR, or MAX and MIN, and decomposing the shift in the DTR seasonal pattern into underlying shifts in the MAX and MIN patterns. More generally, our results may prove useful for assessing and improving structural climate models. Previous research shows that DTR is a useful metric to help assess the accuracy and degree of fit of global climate models (Braganza et al., 2010;Zhou et al., 2010;Lewis and Karoly, 2013;Rader et al., 2018). They generally found that these models persistently underestimated the trend in DTR, which was likely related to deficiencies in modeling water vapor and cloud cover processes. Our new results on the evolving seasonality of DTR may provide an additional, more refined, benchmark for such evaluations. Our results may also prove useful for assessing financial market efficiency, that is, for assessing whether the temperature forecasts embedded in financial asset prices accurately reflect temperature's underlying dynamics. It may be of interest to extend existing work based on AVG (Schlenker and Taylor, 2019) to incorporate our more complete model of AVG dynamics, or to consider multivariate modeling of AVG and DTR, extending the univariate approach undertaken in this paper. Notes to table: All results are based on daily data, 1960-2017. Column 1 reports measurement station by airport code. Column 2 reports the estimated trend movement over the entire 57-year sample in degrees Fahrenheit, using a simple regression on linear trend. The remaining columns report results from the conditional-mean regression (4). p(nt) is the robust p-value for a Wald test of no trend (all coefficients on TIME and D•T IM E interactions are 0), p(ns) is the robust p-value for a Wald test of no seasonality (all coefficients on D's and D•T IM E interactions are 0), and p(nts) is the robust p-value for Wald a test of no trend in seasonality (all coefficients on D•T IM E interactions are 0). ρ is the estimated autoregressive coefficient, and R 2 is the coefficient of determination. Asterisks denote significance at the one percent level. See text for details. Notes to table: See Table1. Notes to table: See Table1

Notes to table: See Table

An interesting and novel exception isSturm et al. (2017), based on options-theoretic thinking. Wigglesworth (2019) finds an important role of DTR in a panel study of U.S. state-level agricultural production over and above standard covariates like capital, labor, and AVG. EViews code is available at https://www.sas.upenn.edu/ ~fdiebold/papers/paper122/DTRcode.txt

The data are available at https://www.ncdc.noaa.gov/ghcn-daily-description. For details, seeMenne et al. (2012) andJaffres (2019). AVG and DTR are standard measures, used almost universally. Of course one could entertain more sophisticated measures of central tendency than AVG, for example, but again AVG and DTR are standard. There is of course no need for an intercept, which would be completely redundant and hence cause perfect multicollinearity. We emphasize that we study only urban areas. They are of maximal interest because a large fraction of the population is concentrated there. However, urban temperature patterns and their evolution (e.g., for AVG or DTR) may (and do) of course differ from those in non-urban areas. There were a (very) few missing observations, in which case we interpolated using an average of the immediately previous and subsequent days' values, rounded to the nearest integer. The missing observations are: BWI max: 1/7/04, min: 1/6/04, DSM max: 9/15/96, min: 9/15/96, and TUS max: 5/10/10, 8/18/17, 8/19/17, min: 5/11/10, 8/18/17, 8/19/17. Detailed regression results for all cities are in the online Appendix ? ? (https://www.sas.upenn.edu/ ~fdiebold/papers/paper122/OnlineAppendix.pdf, and underlying EViews code is at https://www.sas. upenn.edu/ ~fdiebold/papers/paper122/DTRcode.txt.