Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
79,710 characters · 14 sections · 63 citation commands
Nonlinear Temperature Sensitivity of Residential Electricity Demand: Evidence from a Distributional Regression Approach
JEL Classification: C51, C53, Q41
\noindentKey words and phrases: Functional regression, electricity demand, temperature sensitivity, distributional data
\setcounter{footnote}{0} \pagenumbering{arabic} \pagestyle{myheadings} \setcounter{page}{0} \thispagestyle{empty}
\baselineskip19.5pt \pagenumbering{arabic}
Accurately estimating electricity demand during extreme weather events is vital for ensuring energy security and enhancing grid reliability. For energy practitioners and policymakers, precise demand estimation is crucial for balancing supply and demand, preventing power outages, and ensuring that the energy infrastructure can withstand increased loads during such critical periods. As climate change drives more frequent and severe weather events, the importance of precise demand estimation has grown significantly, making it a cornerstone of sustainable energy management.
The relationship between temperature and electricity demand is inherently nonlinear, as demand typically rises sharply during hot or cold conditions. Traditional approaches to modeling this relationship rely on nonlinear temperature variables such as Heating Degree Days (HDD) and Cooling Degree Days (CDD). However, these methods have notable limitations, including the need to arbitrarily define threshold temperatures. To address this limitation, nonlinear econometric models or density-based approaches that incorporate the full temperature distribution as an explanatory variable have been investigated. These methodologies effectively capture the seasonal dynamics of electricity demand while allowing the threshold temperature to be estimated endogenously within the model.
Nonetheless, unexpected spikes in electricity demand during extreme temperature events present crucial challenges for grid stability, necessitating more advanced econometric modeling techniques. Residential electricity demand, in particular, exhibits heightened sensitivity to temperature extremes, increasing unpredictably under such conditions. Existing nonlinear and density-based econometric approaches often fail to adequately capture these nonlinear demand responses, leading to significant underestimation during critical periods. As climate change amplifies temperature extremes, reassessing and refining demand estimation methodologies is essential for ensuring grid stability and optimizing operational efficiency.
Beyond the inherent challenges of modeling nonlinear temperature-demand relationships, a significant source of error arises from measurement and estimation inaccuracies. In particular, constructing a representative temperature regressor from high-frequency or regional raw data introduces estimation errors, leading to biased and inconsistent demand response estimates. Additionally, the density-based approach assumes that nonparametrically estimated temperature distributions represent true distributions. This assumption further amplifies estimation bias, as errors in the temperature distribution can distort demand response projections. These errors are particularly problematic during extreme temperature events, where precise demand forecasting is crucial for ensuring energy security and grid stability.
Building on the foundational work of engle1986semiparametric, a substantial body of literature has applied the cointegration framework for analyzing electricity demand, with noticeable contributions from engle1989merging, silk1997short, beenstock1999demand, chang2014time, chang2016disentangling, and chang2021forecasting, among others. Subsequently, the literature has focused on incorporating nonlinear temperature variables, such as HDD and CDD, or employing nonlinear econometric methods to evaluate the impact of temperature fluctuations on electricity demand, including sailor1997sensitivity, henley1997non, henley1998residential, valor2001daily, pardo2002temperature, and apadula2012relationships. In particular, smooth transition, panel threshold, and switching regression models have been employed to capture nonlinear demand responses (moral2005modelling; bessec2008non). The nonlinear relationship between electricity demand and temperature has been further explored through various advanced methodologies. Notably, artificial neural networks have been applied to improve load forecasting accuracy (teixeira2017; caro2020; sharma2020), while machine learning approaches have been employed to enhance short-run demand predictions (al2019short; sultana2022).
Rather than analyzing the overall effects of weather variables, recent literature has increasingly focused on the impact of extreme temperatures on peak electricity demand. Early studies employed seasonal autoregressive moving average models with conditional heteroskedasticity to model fluctuations in peak demand (sigauke2011prediction; rallapalli2012forecasting). Later research incorporated tail-quantile estimation and extreme value theory to improve the accuracy of predicting extreme peak demand (sigauke2013extreme; chan2015extreme; sigauke2017modelling; sigauke2020modelling). More recently, density-based functional approaches have emerged to analyze these extreme effects across sectors, providing deeper insights into the nonlinearity between temperature and electricity demand (chang2016new; miller2022modeling).
Compared to the existing literature, we investigate distribution-to-scalar regression models that effectively account for the nonlinear response of residential electricity demand, even in the presence of measurement errors and extreme temperature events—an aspect that has not been comprehensively explored. By addressing the complexities introduced by measurement errors, our approach provides a more robust estimation, setting it apart from traditional techniques that may overlook these critical factors. Furthermore, our approach differs significantly from existing studies by focusing on the nonlinear response of residential electricity demand to changes in temperature distribution, rather than just temperature levels. This allows for a more comprehensive analysis, capturing the broader distributional effects of temperature fluctuations, particularly during extreme temperature events.
Accordingly, we incorporate both density and hazard rate functions (as distributional summaries of temperature), each offering distinct advantages, to construct a comprehensive analytical framework. The density-based approach provides direct interpretability within a compositional data framework, facilitating a clear assessment of temperature's impact on electricity demand across the entire temperature distribution. In contrast, the hazard rate is particularly well-suited for analyzing extreme temperature events, as it incorporates the conditional probability of such occurrences to some extent, thereby enhancing the precision of demand estimations under extreme conditions. To ensure that distributional information is properly incorporated into our econometric framework, we will consider relevant transformations of the temperature densities and hazard rates (see Section (ref)). For comparison, we additionally consider the quantile-to-demand model, which is not only a popular approach for utilizing distributional information (see, e.g., yang2020quantile) but also aligns conceptually with density- or hazard-based models.
The remainder of this paper is structured as follows. Section (ref) outlines the econometric methodology used to estimate the nonlinear response of residential electricity demand while addressing measurement errors. Section (ref) describes the residential electricity demand and local temperature data from the Republic of Korea. Section (ref) presents the empirical temperature response functions based on density-, hazard-, and quantile-based approaches under hypothetical scenarios and two historical extreme temperature events. This section also includes a discussion of the estimation results derived from these three distributional predictors. Finally, Section (ref) provides concluding remarks, with proofs provided in the Appendix.
We let $y_t$ represent electricity demand, a practical measurement of which for our empirical study is detailed in Section (ref), and let $X_t^\circ$ be a function that summarizes the distributional properties of temperature. A possible (and popular) candidate for $X_t^\circ$ may be the probability density function (PDF) of temperature, denoted hereafter $\phi_t$. However, as well documented in the literature, directly using the PDF as a predictor is not advisable in applications of statistical methods developed within a standard Hilbert space setting for functional data analysis. As proposed by petersen2016, we thus consider appropriate transformations of $\phi_t$. There are various potential choices for $X_t^\circ$, but in this paper, we mainly focus on the following: (i) the Centered-Log-Ratio (CLR) transformation and (ii) the Log-Hazard-Rate (LHR) transformation of $\phi_t$. Assuming that the PDF of temperature, $\phi_t$, is supported on $[a,b]$, the CLR transformation $X_t^\circ$ is defined as follows:
Under some mathematical conditions on the underlying PDF $\phi_t$, the transformation from $\phi_t$ to its CLR not only embeds the probability densities on $[a,b]$ into a linear subspace of the $L^2[a,b]$-Hilbert space (the Hilbert space of square-integrable functions on $[a,b]$) but also is invertible; see e.g., Egozcue2006 and Boogaart2014. It is known that the inverse CLR transformation is given by $X_t^\circ (r) \mapsto \exp(X_t^\circ (r)) \big/ \int_{a}^b \exp(X_t^\circ (r)) dr$ for $r \in [a,b]$, which recovers $\phi_t$. The CLR transformation of $\phi_t$ has been considered a standard choice in models involving distributional functional data, as discussed in several recent articles to be mentioned shortly.
As an alternative, we also consider the LHR transformation of $\phi_t$, which is given as follows: for some small positive $\delta>0$,
where $\Phi_t(s)$ is the cumulative distribution function (CDF) given by $\Phi_t(s)=\int_{0}^s \phi_t(r)dr$. Note that $h_t(r)=\phi_t(r)/(1-\Phi_t(r))$ is the so-called hazard rate of $\phi_t$ and for this to be well defined, we need to restrict the support into $[a,b-\delta]$ for some small positive $\delta$. It is known that the LHR is also invertible transformation of $\phi_t$ into a linear space under appropriate conditions (see petersen2016).
It appears to be more common to consider the CLR transformation or its similar alternatives in the literature on density-valued functional data (see, e.g., kokoszka2019forecasting, SEO2019, seong2021functional), and this is partly a reason that we also consider it in the present paper. On the other hand, cases using the LHR transformation seem scarce in the literature, to the best of the authors' knowledge, except for petersen2016, who proposed the LHR as a way to avoid potential issues when directly using the PDF. As a slight and obvious modification of the LHR, we define the following, obtained by replacing the hazard rate ($\phi_t(r)/(1-\Phi_t(r))$) with the reversed hazard rate (see block1998reversed), defined by $\phi_t(r)/\Phi_t(r)$, as follows:
We hereafter call the above as the Log-Reversed-Hazard-Rate (LRHR) transformation.
In the sequel, we will consider a functional linear model connecting the electricity demand $y_t$ with the distributional predictor $X_t^\circ$. It is worth noting that while the CLR, LHR, and LRHR transformations characterize the temperature distribution, they differ in how they summarize distributional information. In the CLR model, the distributional predictor is fully characterized by the PDF $\phi_t(s)$, which can be approximated as the average change of the distribution function of temperature on a small interval $[s,s+\Delta)$ for small $\Delta>0$, i.e, $\phi_t(s)\approx \Delta^{-1}\mathbb{P}(s \leq S_t < s+\Delta)$, where $S_t$ denote the random variable of temperature level at time $t$. On the other hand, the hazard (resp.\ reversed) rate may be approximately understood as ${\Delta}^{-1}{\mathbb{P}(s \leq S_t < s+\Delta | S_t \geq s)} = {\Delta}^{-1}{\mathbb{P}(s \leq S_t < s+\Delta)}/{\mathbb{P}(s \geq S_t)}$ (resp.\ $\Delta^{-1}{\mathbb{P}(s-\Delta < S_t \leq s | S_t\leq s)}={\Delta}^{-1}{\mathbb{P}(s \leq S_t < s+\Delta)}/{\mathbb{P}(S_t\leq s)}$). For any $s$, ${\mathbb{P}(S_t \geq s)}$ (resp.\ ${\mathbb{P}(S_t\leq s)}$) indicates the exposure probability of the temperature being higher (lower) than $s$ at time $t$, and thus the hazard (resp.\ reversed hazard) rate, which we consider, may be understood as the PDF weighted by this exposure probability. If we consider any small density shock, which is particularly concentrated on the upper (lower) tail of the temperature PDF and hence may have a significant impact on the electricity demand despite its small magnitude,\footnote{We are mainly interested in the effect of extreme temperature events on electricity demand, and such events are expected to happen with low probability; thus those will only make small changes in the PDF in the lower or upper tails.} then the LHR (resp.\ LRHR) will be much more responsive to this change than the PDF. It is thus expected that the LHR (resp.\ LRHR) will have a higher tendency to comove with electricity demand, particularly during extremely high (resp.\ low) temperature events, and hence, the model involving the LHR (resp.\ LRHR) may be more suitable for capturing and explaining the effect of extreme temperature events at high (resp.\ low) quantiles on electricity demand.
In the literature on density-valued functional data, it seems to be popular to regard the quantile function corresponding to $\phi_t$ as a distributional summary (see e.g., yang2020quantile) and use it for statistical inference. We also briefly discuss this choice of distributional predictor in Section (ref).
In practice, $\phi_t$ and $\Phi_t$ cannot directly be observed, and thus the distributional predictor $X_t^\circ$, considered in Section (ref) is not observable. For empirical analysis, it must be constructed from the raw temperature data. For example, in the CLR case, we may replace $\phi_t$ or $\log \phi_t$ with its relevant nonparametric estimate in (ref) (as in SEO2019). However, this replacement necessarily introduces estimation/smoothing errors, and in this case, it is well known that standard estimation methods for functional linear models may lead to inconsistent estimation of the model; see Benatia2017, Chen_et_al_2020, and seong2021functional for a more detailed discussion. Section (ref) details how the distributional predictors used in our empirical analysis are constructed.
Subsequently, we discuss econometric methods for studying the distributional impact of temperature on electricity demand that are robust to these estimation/smoothing errors, which naturally exist, in the distributional predictor.
For the subsequent discussion, it will be convenient to introduce some notation. We let $\mathcal H$ denote $L^2[a,b]$ (i.e., the Hilbert space of square-integrable functions defined on $[a,b]$) with inner product $\langle g,h \rangle = \int_{a}^b g(s)h(s)ds$ for and norm $\|g\|=\sqrt{\langle g,g \rangle}$, where $g,h \in \mathcal H$; $\mathcal H$ is a commonly considered Hilbert space in the literature on functional data analysis. The distributional predictor $X_t^\circ$ and its measurement $X_t$, which are introduced in (ref), are understood as random elements in $\mathcal H$. Essential concepts on $\mathcal H$, relevant to the subsequent discussion, are reviewed in Section (ref).
To study temperature sensitivity of electricity demand, we consider the following model:
where $\mathbb{E}[ \varepsilon_t] = 0$ and $f:\mathcal H\mapsto \mathbb{R}$ is a linear map, connecting $y_t$ and $X_t^{\circ}$.\footnote{The linear map $f$ can equivalently be expressed as an integral transformation associated with a kernel function $\mathrm k_f$, as follows: $ f(X_t^{\circ}) = \int_{a}^b \mathcal \mathrm k_f (r) X_t^{\circ}(r)dr$. Obviously, estimation $\mathrm{k}_f$ is identical to estimation of $f$. In this paper, we focus on $f$ but the subsequent results can be rephrased for $k_f$ with only a minor modification.} Equation (ref) corresponds to the Function-to-Scalar Regression (FSR) model, which is widely discussed in the literature. Following the terminology used in conventional linear models, $f$ is typically referred to as the (functional) coefficient or slope parameter. However, in our context, since $f$ captures how the distributional properties of $X_t^\circ$ influence $y_t$, we refer to it as the distributional coefficient hereafter, and we refer to (ref) as the Distribution-to-Scalar Regression (DSR) model.
As discussed in Section (ref), in practice, $X_t^\circ$ is not directly observed and has to be estimated or smoothed from discrete samples $\{\phi_t^{\circ} (s_i)\}_{i=1}^{N_t}$ of $\phi_t$, and practitioners often have no choice but to replace $X_t^{\circ}$ with its feasible estimate $X_t$. As earlier pointed out by seong2021functional, this is a common source of endogeneity in the regression model involving functional data, and ignoring such errors may be detrimental to the use of standard estimators of $f$, requiring the exogeneity of the predictor, developed in the literature. Even if more detailed discussion can be found in recent articles (e.g., Chen_et_al_2020,seong2021functional), it will be helpful to illustrate this issue with a concrete example for the subsequent discussion.
Suppose that $X_t = X_t^{\circ} + e_t$, where $e_t$ represents an additive deviation of $X_t$ from $X_t^{\circ}$ and satisfies $\mathbb{E}[e_t] = 0$. From (ref), we obtain the following relationship:
where $\mathbb{E}[u_t] = 0$. Given that $X_t^\circ$ is unobserved, while $X_t$ is observable, practitioners will work with (ref) rather than (ref). In (ref), however, $X_t$ and $u_t$ are generally correlated (meaning that $\mathbb{E}[\langle X_t,x \rangle u_t] \neq 0$ for some $x\in \mathcal H$) and thus $X_t$ is endogenous. As is well documented in the aforementioned articles, this endogeneity makes the standard estimators, developed under the exogeneity condition ($\mathbb{E}[\langle X_t,x\rangle u_t]=0$ for every $x\in \mathcal H$), inconsistent and invalidates the inferential methods based on those estimators (see also Benatia2017).
To estimate the distributional coefficient $f$ from (ref) and implement valid statistical inference in our DSR model, which inevitably involves measurement errors, we adapt the functional IV approach, considered in seong2021functional for Function-to-Function Regression (FFR) models, to our context. To this end, we let $Z_t$ be the functional Instrumental Variable (IV), which satisfies the following:
$C_{XZ}$ (resp.\ $C_{Zu}$) defined above is called the (cross-)covariance operator of $X_t$ and $Z_t$ (resp.\ $Z_t$ and $u_t$). The conditions given in (ref) and (ref) consist of only the minimal requirements for the IV $Z_t$ in our DSR model. More conditions for the consistency and (local) asymptotic normality of our estimator will be detailed in Section (ref) of the Appendix. Assuming that the measurement errors are not serially correlated and that $X_t$ tends to exhibit time series dependence in the considered empirical model, we may let $Z_t$ be a lagged distributional predictor in our empirical analysis, as in seong2021functional; this is also a special case of the IV considered in Chen_et_al_2020. To obtain our proposed estimator, we need to compute the sample covariance operators $\widehat{C}_{XZ}$ and its adjoint $\widehat{C}_{XZ}^\ast$ defined as follows: for any $x \in \mathcal H$,
where $\bar{X}_T=T^{-1}\sum_{t=1}^T X_t$ and $\bar{Z}_T=T^{-1}\sum_{t=1}^T Z_t$. $\widehat{C}_{XZ}$ and $\widehat{C}_{XZ}^\ast$ are, respectively, the sample counterparts of $C_{XZ}$ and its adjoint $C_{XZ}^\ast$ (defined by $C_{XZ}^\ast(x) = \mathbb{E}[\langle Z_t-\mathbb{E}[Z_t],x \rangle(X_t-\mathbb{E}[X_t])])$. As is well known in the literature, the composite map $\widehat{C}_{XZ}^\ast \widehat{C}_{XZ}$ allows the eigendecomposition with respect to its nonnegative eigenvalues (hereafter denoted $\{\widehat{\lambda}_j^2\}_{j\geq 1}$) and the corresponding eigenvectors (denoted $\{\widehat{g}_j\}_{j\geq 1}$). That is,
where $\widehat{\lambda}_1^2\geq \widehat{\lambda}_2^2 \geq \ldots \geq 0$. The eigenelements of $\widehat{C}_{XZ}^\ast \widehat{C}_{XZ}$ can be computed using the standard functional principal component analysis (FPCA) as in Bosq2000 and seong2021functional. Then our proposed estimator $\widehat{f}$, which is a map from $\mathcal H$ to $\mathbb{R}$, is defined as follows: for any $x \in \mathcal H$,
where $\bar{y}_T=T^{-1}\sum_{t=1}^T y_t$ and ${\rm K}$ is sufficiently smaller than $T$; a more detailed and technical requirement on the choice of ${\rm K}$ will be postponed to Section (ref) of the Appendix. If ${\rm K}$ is determined, it is straightforward to compute the estimator using $\widehat{C}_{XZ}$ and the sample eigenelements $\{\widehat{\lambda}_j^2,\widehat{g}_j\}$. The above estimator turns out to be viewed as a sample analogue estimator based on the population equation $C_{Zy}C_{XZ}= f C_{XZ}^\ast C_{XZ}$, which holds when $Z_t$ satisfies (ref) (see Section (ref) of the Appendix for more detailed discussion). From this equation, the estimator (ref) is obtained by (i) replacing the population operators $C_{Zy}$ and $C_{XZ}$ with their sample counterparts, $\widehat C_{Zy}$ and $\widehat C_{XZ}$, respectively, and then (ii) employing a regularized inverse of $\widehat{C}_{XZ}^\ast \widehat{C}_{XZ}$, which is necessary to compute the estimator. Through our study of the asymptotic properties of the proposed estimator, we establish the consistency of $\widehat{f}$ under appropriate assumptions as follows:
Proposition (ref) is presented without the mathematical details intentionally; a more rigorous version is provided in the Appendix (see Proposition (ref)). The consistency stated in Proposition (ref) implies that for every $x \in \mathcal{H}$, $\widehat{f}(x)$ converges in probability to $f(x)$.
It is often of interest to practitioners to make inference on $f(\zeta)$, which may be understood as the effect of an additive distributional perturbation $\zeta$ given to $X_t^\circ$ on $y_t$. {Specifically, in our empirical analysis, $\zeta$ will be set to various extreme temperature events and we are interested in the effect of these events on $y_t$.} The consistency result in Proposition (ref) implies that this effect can be estimated by $\widehat{f}(\zeta)$. Beyond the consistent estimation, the following result, which can be used for statistical inference on $f(\zeta)$, is established under some appropriate conditions (see Proposition (ref) in the Appendix):
where $\widehat{\sigma}^2_u = T^{-1} \sum_{t=1}^T \widehat{u}_t^2$ with $\widehat{u}_t = y_t-\widehat{f}(X_t)$ and
The result (ref) may be understood as the local asymptotic normality of the functional IV estimator, which was earlier obtained by seong2021functional in the FFR model.
It is well known that $f(\zeta)$ for any $\zeta$ may be understood as the integral transformation
for some uniquely identified function $\psi_f \in \mathcal H$. $\psi_f(s)$ is sometimes interpreted as the temperature effect when the distributional predictor $X^\circ$ hypothetically concentrated at a point $s$ (chang2014time). It might be of interest to practitioners to implement statistical inference on $\psi_f(s)$, but in our setup, it is not possible to implement direct statistical inference with a perturbation which is fully concentrated at a point.\footnote{In the standard $L^2[a,b]$-Hilbert space setting, such a function is essentially equivalent to the zero function.} However, it is, instead, possible to consider a weighted average of $\psi_f(s)$ near $s$. For example, if we let $\zeta_s(r) = w_h(|r-s|)/{\int_{a}^{b} w_h(|r-s|)dr}$, then statistical inference on the following quantity is feasible based on the asymptotic normality result (ref):
where $w_h(|r-s|)$ is the standard kernel function defined near $s$ with bandwidth $h$. The above $\tilde{\psi}_f(s)$ is dependent on the choice of kernel function, but it may be understood as the local weighted average of $\psi_f$ near $s$ for any reasonable choice of the kernel function. We hereafter call $\tilde{\psi}_f(s)$ as the kernel-weighted temperature response function.
To investigate the nonlinear temperature sensitivity of electricity demand, we analyze an extensive dataset that includes residential electricity demand (in GWh) and local temperature (in degrees Celsius) readings from the Republic of Korea. The dataset spans from January 1999 to December 2023, with 300 monthly observations. The residential electricity demand data is provided by the Korea Electric Power Corporation (KEPCO), the state-owned utility responsible for electricity sales across the country. Due to the inherent variability in the number of days within each month, we have normalized the raw monthly electricity demand by dividing it by the number of effective days per month. This adjustment produces a metric known as electricity demand per effective day, as discussed in chang2014time, which underscores the importance of accounting for temporal differences in energy demand analysis.
chang2016new found no empirical evidence suggesting that non-climate variables significantly influence the demand response to temperature in the residential sector. Building on this finding, and to isolate the short-run fluctuations in residential electricity demand, we calculate the deviation from the 12-month moving average of the natural logarithm of electricity demand per effective day, as illustrated in Figure (ref). This deviation serves as a proxy for the short-run component of electricity demand, allowing us to focus on the short-run effects of temperature changes while filtering out the long-run effects of income and electricity price. By employing this detrending approach, we ensure that our analysis effectively captures the nonlinear responses of residential electricity demand to temperature fluctuations, independent of other non-climate factors (chang2016new).
To model the seasonality, we retrieve the hourly temperature data from the KMA National Climate Data Center.\footnote{Downloaded from https://data.kma.go.kr/cmmn/main.do on July 1, 2024.} We then construct the monthly temperature PDF over the common support, $[-20,40]$, using hourly temperature observations from the five largest cities (Seoul, Daejeon, Daegu, Gwangju, and Busan) covering the territory of the Republic of Korea.\footnote{Instead of estimating the range of the common support, we adopt the common support of $[-20,40]$ with a grid difference, $0.1$, as provided in the literature (chang2014time; chang2016new).} To represent the temperature PDF of the Republic of Korea, we calculate the weighted average of the estimated time series for these five PDFs, assigning weights based on the annual proportion of electricity sales in each of the five regions: the Seoul Metropolitan Area (Seoul, Incheon, Gyeonggi), Daejeon, Daegu, Gwangju, and Busan.\footnote{The combined annual electricity consumption of the five regions accounts for approximately 55.3% of the national total on average during the period from 1999 to 2021. Due to data availability constraints, the proportions for the years 2022 and 2023 have been substituted with the values from 2021.}\footnote{The annual regional consumption data was obtained from the Korea Energy Statistics Information System (KESIS) and can be accessed at \url{https://www.kesis.net/sub/subChart.jsp?report_id=33150&reportType=0} (data downloaded on July 1, 2024).}
Figure (ref) illustrates the time series of temperature PDFs, estimated using Gaussian kernel and Silverman's rule of thumb bandwidth, and the dynamics of their descriptive statistics. Note that the generated monthly temperature PDFs exhibit significant heterogeneity, particularly between summer and winter. The time series of monthly temperature distributions in the Republic of Korea from 1999 to 2023 reveals stable seasonal cycles in the mean temperature, reflecting consistent annual patterns. However, the standard deviation indicates increasing variability, particularly post-2010, suggesting a growing spread in temperature values around the mean. The skewness exhibits seasonal fluctuations, with positive skewness during warmer months indicating a higher frequency of extreme high temperatures, and negative skewness during colder months reflecting an increase in extreme low temperatures. However, this skewness does not show a clear trend of increasing asymmetry over time but rather indicates a recurring seasonal pattern. The kurtosis remains relatively stable but experiences significant spikes, particularly after 2010. These spikes suggest periods with more frequent and intense extreme temperature events. Thus, while the skewness reflects consistent seasonal variability, the spikes in kurtosis indicate a rise in temperature extremes, likely influenced by climate change.
As discussed in Section (ref), distributional predictors in practice need to be estimated from the raw data. For the CLR transformations, we use the estimated densities, denoted hereafter as $\widehat{\phi}_t$, obtained in Section (ref) (see Figure (ref)), and construct $X_t$ as follows:
For the model using the LRH and LRHR predictors, $X_t$ can be constructed as
where $\widehat{\Phi}_t(r) = \int_{0}^{r} \widehat{\phi}_t(s)ds$. However, in practical computation for the LHR predictor, instead of $1-\widehat{\Phi}_t(s)$, we take the maximum of a small positive constant ($10^{-3}$ in the subsequent analysis) and $1-\widehat{\Phi}_t(s)$ in (ref), which is done to avoid the case where $1-\widehat{\Phi}_t(s)$ equals zero, and hence the hazard rate is undefined. The LRHR predictor is similarly constructed by taking the maximum of a small positive constant and $\widehat{\Phi}_t(s)$.\footnote{As an alternative, we also considered adding 0.001 to $1-\widehat{\Phi}_t(s)$ or $\widehat{\Phi}_t(s)$, but we found that the estimation results change little.} Note that this convenient method of constructing the LHR and LRHR results in the constructed predictors not perfectly matching their theoretical counterparts ((ref) and (ref)). However, for our analysis, we only require that the constructed predictors serve as reasonable proxies, and thus, this does not invalidate the empirical results that will be discussed.
In this section, we estimate the DSR model (ref) with potential endogeneity, using the distributional predictors described in (ref) and (ref). For the computation of our estimator with finite samples, we represent the distributional predictor using 80 orthonormal Fourier basis functions on the support.\footnote{Moderate changes in the number of basis functions result in only minor numerical differences in the estimation results.} Noting that measurement errors mostly arise from smoothing \( X_t \) from its discrete realizations, we assume that these errors are serially uncorrelated. Under this assumption, and given that the sequence of \( X_t^\circ \) (and also \( X_t \)) is serially correlated, the lagged variable \( X_{t-\kappa} \) for \( \kappa \geq 1 \) may serve as a candidate IV for our empirical analysis. Similar approaches were previously adopted in Chen_et_al_2020 and seong2021functional concerning dependent functional data. Since our approach relies on a high degree of dependence, where possible, of the IV on the endogenous distributional predictor \( X_t \), as in standard IV methods (see Section (ref) for more details), it is not advisable to use a large value of \( \kappa \), as this could make \( X_{t-\kappa} \) only weakly correlated with \( X_t \). Therefore, we choose to set \( \kappa \) to 1. We believe that this selection allows us to account for potential errors in the density/quantile estimation process or data quality while enhancing the robustness and reliability of our analysis. Lastly, in our estimation procedure, ${\rm K}$ (see (ref)) is set to 3.\footnote{Following the results given in Section (ref) of the Appendix and Remark (ref), we set ${\rm K}=1+\max_{j\geq 1}\{\tilde{\lambda}_j^2 \geq \alpha\}$ for $\alpha = 0.01 \times T^{-0.2}$, where $\tilde{\lambda}_j^2=\lambda_j^2/\sum_{j=1}^\infty \widehat{\lambda}_j^2$ and this is to obtain a choice of ${\rm K}$ which is independent of the scale of $X_t$ and $Z_t$; see Section S5 of seong2021functional.}
We hereafter focus on statistical inference regarding the impact of temperature events on electricity demand, given by $f(\zeta)$, where $\zeta$ denotes the distributional shocks applied to the distributional predictor. This quantity is estimated using our proposed estimator, $\widehat{f}(\zeta)$, and is accompanied by an asymptotically valid confidence interval constructed based on the results in (ref). Taking advantage of the distributional flexibility of our DSR model, we examine hypothetical scenarios and two historical events for the specification of $\zeta$, as detailed in the following section.
In this paper, we present two temperature response functions. The first is the kernel-weighted temperature response function (see (ref)), which serves as a benchmark and is similar to those estimated by chang2014time and chang2016new. It captures the point-wise temperature sensitivity across the extreme temperature domain. The second is derived from historical extreme temperature events. While the benchmark temperature response function provides intuitive insights into temperature sensitivity by isolating the effects of temperature changes from distributional shifts, the response function based on historical extreme events enhances the practical relevance of extreme temperature fluctuations.
More specifically, the key distinction between these two approaches lies in their treatment of distributional changes. The analysis of historical extreme events explicitly accounts for shifts in the distribution as normal temperature distributions transition toward extreme ones. In contrast, the benchmark temperature response function hypothetically assumes that the temperature stays at a specified level without consideration of the distributional shift that appeared in the past. Although we provide representative cases of extreme temperature events, the specific distributional characteristics of past cold or heat waves may not necessarily recur. Therefore, the difference between the benchmark response function and the date-specific response function can be interpreted as the additional distributional effect unique to that date.
To estimate the benchmark response function, $\tilde{\psi}_f(s)$ (see (ref)), we employ Gaussian kernel, which is truncated on the support; more specifically, for each temperature level \( s \), $w_h(|r-s|)$ is defined as follows: \[ w_h(|r-s|) = \frac{1}{\sqrt{2\pi}} \exp\mathopen{}\mathclose\bgroup\originalleft\{-0.5\mathopen{}\mathclose\bgroup\originalleft(\frac{|r-s|}{h}\aftergroup\egroup\originalright)^2\aftergroup\egroup\originalright\} \times {1}\{ -20 \leq s \leq 40\}, \] where the bandwidth \( h \) is set to span 15°C around \( s \) to align with the range of density increases observed in cold (17.3°C) and heat (12.2°C) wave shocks, as illustrated by the black lines of Figure (ref).
For the response function based on historical extreme events, we consider two significant extreme weather events that occurred in the Republic of Korea during the sample period.\footnote{See e.g., \url{https://www.wunderground.com/cat6/Hottest-Day-Korean-History} and \url{https://watchers.news/2021/01/11/historic-cold-wave-and-heavy-snow-hit-south-korea/}.} In August 2018, the Republic of Korea experienced one of its most intense heat waves on record, with temperatures soaring to an average of 29.2°C, the highest since record-keeping began in 1907. This extreme heat led to the residential electricity demand of 8,851.04 GWh during that month, underscoring the strain placed on the energy grid by cooling needs. Similarly, January 2021 was marked by a significant cold wave, with temperatures plunging to an average of -6.1°C, well below the seasonal norm. During this period, residential electricity demand also spiked to 7,163.28 GWh as households increased heating usage to cope with the severe cold.
Figure (ref) illustrates the year-over-year changes in monthly temperature PDFs associated with extreme weather events, specifically the heat wave of August 2018 and the cold wave of January 2021. In the left panel of Figure (ref), the blue solid line represents the temperature PDF for January 2021, compared to the blue dashed line for January 2020, with the black line showing the (vertical) difference between the two PDFs, capturing the effect of the cold wave. Similarly, in the right panel of Figure (ref), the red solid line represents the temperature PDF for August 2018, while the red dashed line shows the temperature PDF for August 2017, with the black line indicating the difference between these two PDFs.
The temperature distributions during the 2018 heat wave and the 2021 cold wave exhibit distinct shifts, indicative of extreme weather conditions. The August 2018 distribution is significantly right-skewed, suggesting a higher frequency of exceptionally high temperatures, with a pronounced right tail representing extreme heat events. In contrast, the January 2021 distribution is heavily left-skewed, highlighting the prevalence of unusually low temperatures, with a long left tail representing severe cold conditions. These graphs reveal significant changes in temperature distributions during these historical extreme events, illustrating how both the heat wave and cold wave deviated sharply from the typical conditions of the preceding months.
Our empirical analysis aims to provide practitioners with insights into the potential increase in future electricity demand when similar extreme temperature events occur. By understanding extreme temperature events as potential perturbations in the distributional predictor \( X_t \), we can estimate their impact on electricity demand. We henceforth refer to the function \( z_{\text{extr}}(s) \), corresponding to the black line in the left (resp.\ right) panel of Figure (ref), as the observed cold wave (resp.\ heat wave). This function is computed by subtracting \( \mathrm{Z}_{\text{norm}}(s) \) from \( \mathrm{Z}_{\text{extr}}(s) \), where \( \mathrm{Z}_{\text{norm}}(s) \) denotes the reference (normal) temperature PDF in January 2020 (resp.\ August 2017), and \( \mathrm{Z}_{\text{extr}}(s) \) denotes the extreme temperature PDF in January 2021 for the cold wave (resp.\ August 2018 for the heat wave). Consequently, \( \mathrm{Z}_{\text{extr}}(s) = \mathrm{Z}_{\text{norm}}(s) + z_{\text{extr}}(s) \).
Due to the flexibility of our model in the choice of potential perturbations, we may also consider fractions of these observed heat and cold waves and their effects on electricity demand; specifically, we let
for some large positive integer \( M \). Note that these fractions of observed heat or cold waves satisfy (i) \( |z_a(s)| \leq |z_{\text{extr}}(s)| \) and (ii) \( \int_{-20}^{40} z_a(s) \, ds = 0 \) for all considered values $a$. Thus, adding these fractions of heat or cold waves to the reference PDF \( \mathrm{Z}_{\text{norm}}(s) \) results in a new PDF, hereafter denoted as \( \mathrm{Z}_a(s)\) (i.e., \(\mathrm{Z}_a(s) = \mathrm{Z}_{\text{norm}}(s) + z_a(s)) \). Since the application of \( z_a(s) \) induces a globally more moderate shock compared to that of \( z_{\text{extr}}(s) \), \( \mathrm{Z}_a(s) \) represents less extreme temperature conditions than the extreme temperature PDF \( \mathrm{Z}_{\text{extr}}(s) \) observed in January 2021 for the cold wave or August 2018 for the heat wave. Thus, \( z_a(s) \) can naturally be interpreted as a possible, but less extreme, heat or cold wave. Observe also that, as $a$ grows, $\mathrm{Z}_a(s)$ constructed from $z_a(s)$ moves away from the normal temperature PDF ($\mathrm{Z}_{\text{norm}}(s)$) and closer to the extreme temperature PDF ($\mathrm{Z}_{\text{extr}}(s)$) by a constant functional increment $(1/M) \times z_{\text{extr}}(s)$ at each step; the left panel of Figure (ref) illustrates $\mathrm{Z}_a$ for a few selected values of $a$.
Similarly, Figure (ref) illustrates the corresponding year-over-year changes in monthly temperature LRHR and LHR for these cold and heat wave events, respectively. Specifically, in January 2021, the LRHR gradually increased as temperatures decreased, peaking at -19.5°C, which represents the highest conditional probability of occurrence at this temperature. This indicates that, given temperatures had already dropped below -19.5°C, the likelihood of further decreases to even lower temperatures was greatest. This peak highlights -19.5°C as a critical threshold during the cold wave events. In contrast, January 2020 exhibited a peak LRHR at -7.6°C, reflecting a focus on milder cold conditions and a rapid decline in LRHR at more extreme lows. Similarly, in August 2018, the LHR peaked at 39.5°C, indicating the dominance of intense heatwave conditions. In comparison, August 2017 exhibited a peak at 36.4°C, indicating a milder concentration of high-temperature events. It is worth noting that the exponentially increasing pattern of the monthly temperature LHR in August 2018 (LRHR in January 2021) indicates that the probability of the temperature staying around 39.5°C (-19.5°C) is higher than for any temperature below (above) this value. This would demonstrate the cumulative temperature effect on electricity demand during peak times to some extent (see Section (ref)). The bottom left (resp.\ right) panel of the figure illustrates the transition from the normal temperature LHR (resp.\ LRHR) to the extreme temperature LHR (resp.\ LRHR).
Given that the quantile function corresponding to the temperature distribution can serve as a distributional summary, it may be of interest to practitioners to examine changes in quantile functions during the considered extreme temperature events. The left panel of Figure (ref) shows the quantile functions corresponding to the normal and extreme temperature distributions and also illustrates the cold wave event, expressed as a quantile shock. As shown, the impact is most pronounced in the lower quantiles, particularly those associated with the coldest temperatures. This reflects that the January 2021 cold wave caused a significant downward shift in the quantile function compared to January 2020, with the most substantial changes at the lowest quantiles. In the middle panel, the quantile shocks associated with the August 2018 heat wave exhibit a relatively uniform increase across the quantiles, with a more pronounced effect at both extremes. This implies that the heat wave not only elevated overall temperature levels but also disproportionately intensified occurrences at both the lower and upper quantiles. The fractions of observed cold and heat waves, expressed as quantile shocks, can be constructed as in (ref). These quantile shocks are used to illustrate the transition from the normal temperature quantile function to the extreme one in each scenario, as shown in the right panel of the figure.
In this section, we analyze the impact of extreme temperature events on residential electricity demand using the PDF-based DSR model, hereafter referred to as the density-to-demand model. We first apply this model to estimate point-wise temperature sensitivity through the benchmark response function. Subsequently, we analyze the effects of observed heat and cold wave fractions on electricity demand using the CLR predictor. Note that we generate the relevant functional shocks $\zeta_a$, corresponding to \( z_a(s) \), in the CLR domain. This can be implemented by applying the CLR transformation to both the reference PDF \( \mathrm{Z}_{\text{norm}}(s) \) and \( \mathrm{Z}_a(s) \), and then computing the subtraction of the former from the latter in each of the heat and cold wave scenarios.\footnote{If $\mathrm{X}_{\text{extr}}$ (resp.\ $\mathrm{X}_{a})$ is the CLR corresponding to $\mathrm{Z}_{\text{extr}}$ (resp.\ $\mathrm{Z}_{a})$, then $\zeta_a=\mathrm{X}_{\text{extr}} - \mathrm{X}_{a}$.} The effect of $\zeta_a$ on electricity demand is then given by $f(\zeta_a)$.
In the left panel of Figure (ref), we present our benchmark response function on the range $s \in [-15,35]$ and its 95% confidence interval. Consistent with the findings of chang2016new, the estimated temperature response function exhibits a nonlinear U-shaped pattern across the temperature domain. However, the estimates with $\kappa=1$ reveal a more pronounced degree of nonlinearity. Note also that, unlike existing approaches, our approach explicitly accounts for the interdependence of temperature changes, ensuring that adjustments at one level are accompanied by compensatory changes in other temperatures to maintain the unit integral requirement of a PDF.\footnote{The function $f(\zeta_s)$, where $\zeta_s$ is defined near $s$, represents the effect of an additional perturbation applied to the distributional predictor $X_t$, as indicated by $f(X_t+\zeta_s) - f(X_t)$. If $X_t$ is the PDF itself, $X_t+\zeta_s$ may not be a valid PDF, as it could violate the unit integral requirement. However, this issue does not arise if $X_t$ is the CLR transformation of a PDF and $\zeta_s$ is a shock expressed in the CLR.}
The right two panels of Figure (ref) present the estimated response of electricity demand on the $y$-axis with respect to $a$, denoting the fraction parameter defined in (ref), of $\mathrm{Z}_a(s)$, depicted on the $x$-axis. As earlier discussed, $\mathrm{Z}_a(s)$ approaches $\mathrm{Z}_{\text{extr}}(s)$ with a constant functional increment as $a$ increases. The solid lines show the estimates, $\widehat{f}(\zeta_a)$, which can be interpreted as the estimated demand responses to progressively but constantly intensifying heat or cold wave shocks, as described in Section (ref). The dotted lines indicate the 95% confidence intervals for $f(\zeta_a)$ based on the asymptotic normality result given in (ref).
The empirical results from the estimation with $\kappa = 1$, which is designed to account for potential measurement errors, indicate an increase in electricity demand response as the considered shock gets closer to the observed cold or heat wave. For comparison, we also report estimation results without accounting for potential endogeneity, with the $y$-axis scaled to half of that in the case of $\kappa=1$ to enhance visualization. Similar to the benchmark response function, the empirical results for observed heat and cold wave fractions with $\kappa=0$ exhibit a similar upward trend. However, the magnitude is substantially lower compared to the case with $\kappa=1$. Since the estimation results for $\kappa=0$ do not account for endogeneity, we focus on interpreting the estimation results with $\kappa=1$ as our main empirical findings.
Note that the right panel of Figure (ref) provides statistical insights into how the considered cold and heat wave events change temperature distributions. As the considered cold wave shock approaches the observed severe cold wave, the mean temperature experiences a linear decline from 2.31°C to -1.53°C, with a substantial increase in variance ($17.30 \rightarrow 51.37$), indicating greater temperature dispersion. Concurrently, skewness decreases ($0.14 \rightarrow -0.04$), and kurtosis declines ($2.62 \rightarrow 2.40$), signifying a more symmetrical and less peaked distribution. In contrast, as the considered heat wave shock approaches the observed severe heat wave, the mean temperature rises linearly from 26.02°C to 28.54°C, with a modest increase in variance ($15.0 \rightarrow 17.11$), rising skewness ($-0.06 \rightarrow 0.20$), and declining kurtosis ($2.90 \rightarrow 2.50$), reflecting a more right-skewed but less peaked distribution. As illustrated in the right panel of Figure (ref), these distributional changes, characterized by greater variance shifts per unit change in mean temperature, reveal that cold wave shocks have a significantly greater impact on electricity demand, with a 28.06% increase, compared to a 9.89% increase observed during heat wave shocks.
The observation that cold wave shocks lead to a greater increase (7.31%) in electricity demand per unit change in mean temperature, following the transition from $\mathrm{Z}_{\text{norm}}$ to $\mathrm{Z}_{\text{extr}}$, compared to heat wave shocks (3.93%) highlights that electricity demand is more sensitive to the extreme and intense temperature conditions associated with cold waves. This implies that demand response is influenced not only by temperature extremes but also by the specific nature of these extremes. This finding departs from existing literature, which predominantly focused on the demand response to changes in temperature levels alone, and also highlights the importance of understanding how shifts in temperature distributions impact electricity demand. Such insights emphasize the need for energy planning that accounts for the varying effects of temperature distribution changes, particularly in the context of increasingly severe cold wave events.
Notably, as the cold waves intensify, electricity demand rises at a diminishing rate as incremental heating needs taper off. More specifically, the demand response values for the cold wave event at fraction parameters $a$ of 0.33, 0.66, and 1.0 are estimated at 0.1677, 0.2328, and 0.2806, respectively. Given that electricity demand is expressed in natural logarithms, these values indicate that a functional shock causing normal temperatures to approach extreme levels by 33% corresponds to a 16.77% year-over-year increase in monthly electricity demand. An additional 33% progression toward extreme temperatures results in a further 6.51% increase, followed by an additional 4.77% increase with a complete transition to extreme conditions. Collectively, the full transition from normal temperatures in January 2020 to extreme temperatures in January 2021 corresponds to a total 28.06% year-over-year increase in monthly electricity demand.
Conversely, electricity demand is estimated to increase at an accelerating rate as the heat waves intensify, driven by heightened cooling requirements. More specifically, the demand response values for the heat wave events at fraction parameters $a$ of 0.33, 0.66, and 1.0 are 0.0222, 0.0498, and 0.0989, respectively. Using the same approach, these results indicate that a functional shock causing normal temperatures to approach extreme levels by 33% corresponds to a 2.22% year-over-year increase in monthly electricity demand. An additional 33% progression toward extreme temperatures results in a 2.76% year-over-year increase, followed by an additional 4.91% increase with the complete transition to extreme conditions. Collectively, the full transition from normal temperatures in August 2017 to extreme temperatures in August 2018 corresponds to a total 9.89% year-over-year rise in monthly electricity demand. This contrasting behavior can be attributed to the extensive reliance on air conditioning during heat waves and the substitution effect between natural gas and electricity for heating purposes during cold waves.
In this section, we analyze the impact of extreme temperature events on residential electricity demand using the DSR model with the LHR (resp.\ LRHR) predictor for both hypothetical and historical heat (resp.\ cold) wave scenarios, collectively referred to as the hazard-to-demand model. Figure (ref) illustrates the benchmark response function $\widetilde{\psi}(s)$, focusing on the extreme temperature ranges of $[-20, 5]$ and $[25, 40]$, and the corresponding estimated electricity demand response as a function of parameter $a$ within the hazard-to-demand model. While the estimation results from the hazard-to-demand model are broadly consistent with those from the density-to-demand model, they provide distinct interpretative insights into how variations in the conditional probability of extreme temperature events distinctly influence electricity demand under cold and heat wave scenarios.
While both models capture increased temperature sensitivity, the hazard-to-demand model produces higher response function values in the extreme heat range but lower values in the extreme cold range compared to the density-to-demand model. Notably, it provides a more reliable representation of the response function at the lower temperature boundary, offering a structured depiction of cold extremes. Specifically, electricity demand rises as temperatures decline within the cold temperature range. However, at extremely low temperatures, demand exhibits diminishing sensitivity to further declines, reflecting saturation effects driven by persistent heating requirements and substitution effects with natural gas. In contrast, in the hot temperature range, rising temperatures correspond to an increased conditional probability of extreme heat events, leading to a steady rise in electricity demand, particularly at higher temperatures. This increased sensitivity is primarily driven by intensive cooling needs, highlighting the asymmetric response of electricity demand to temperature extremes. Further discussion on these distinctions is provided in the following section.
Although with a smaller overall magnitude, the hazard-to-demand model aligns with the density-to-demand model for historical extreme temperature events. The cold wave shocks result in a greater increase in electricity demand per 1°C change in mean temperature (3.01%) following the transition from $\mathrm{Z}_{\text{norm}}$ to $\mathrm{Z}_{\text{extr}}$, compared to the heat wave shocks (1.80%). Moreover, as the cold waves intensify, the rate of increase in electricity demand diminishes, reflecting tapering incremental heating needs. Specifically, the LRHR-based demand response values for the cold wave events at fraction parameters $a$ of 0.33, 0.66, and 1.0 are estimated at 0.0987, 0.1113, and 0.1155, respectively. These estimates indicate that a functional shock causing normal temperatures to shift toward extreme cold levels by 33% corresponds to a 9.87% year-over-year increase in monthly electricity demand. An additional 33% progression toward extreme cold conditions results in a further 1.26% increase, followed by an additional 0.42% increase with a complete transition to extreme cold conditions.
For the heat wave events, the LHR-based demand response values at fraction parameters $a$ of 0.33, 0.66, and 1.0 are estimated at 0.0111, 0.0228, and 0.0453, respectively. These results indicate that a functional shock causing normal temperatures to shift toward extreme heat levels by 33% corresponds to a 1.11% year-over-year increase in monthly electricity demand. An additional 33% progression toward extreme heat conditions results in a further 1.17% increase, followed by an additional 2.25% increase with a complete transition to extreme heat conditions. Both the density-to-demand and hazard-to-demand models, despite differences in their estimated magnitudes, consistently reveal distinct demand response patterns for the heat wave and cold wave events. The cold wave shocks result in more significant demand increases overall while the heat wave shocks exhibit a stronger incremental effect as temperatures approach extreme levels.
Lastly, similar to the density-to-demand model, the empirical results from the hazard-to-demand model for $\kappa=0$ exhibit an upward trend but with a substantially lower magnitude compared to the case with $\kappa=1$. As the estimation results for $\kappa=0$ do not account for endogeneity, they are not considered for interpretation.
Our empirical findings for extreme temperature events provide more practical insights than those derived from existing functional approaches. A key limitation of conventional methods is their tendency to yield inaccurate estimates at the boundaries, particularly in data-sparse regions near critical temperature thresholds such as -20°C and 40°C. This boundary issue can lead to misinterpretations of demand responses at extreme temperatures, potentially resulting in significant misjudgments in energy policy and management.
By leveraging both representative heat and cold wave scenarios, as well as demand responses at specific temperature levels, our approach provides broader distributional insights into electricity demand dynamics under extreme temperature events. Specifically, we incorporate higher-order distributional changes associated with extreme temperatures, illustrating how such conditions result in significantly greater increases in electricity demand compared to normal scenarios--an effect that has been largely overlooked in previous studies. Additionally, we integrate the LHR and LRHR predictors, which are expected to more effectively capture electricity demand dynamics under extreme temperature conditions, into the hazard-to-demand model, a framework not previously explored in the literature. Furthermore, our estimation results remain robust to endogeneity concerns arising from measurement errors in the distributional predictors.
The estimation results from both models, as presented in Figures (ref) and (ref), indicate that the benchmark response function from the hazard-to-demand model more effectively captures electricity demand dynamics under extreme temperature conditions than that of the density-to-demand model. In the extreme cold range, the benchmark response function from the hazard-to-demand model (top left panel of Figure (ref)) remains non-decreasing, except below $-16.2$°C, a threshold not observed in historical cold wave events. In contrast, the benchmark response function from the density-to-demand model (left panel of Figure (ref)) begins to decline at $-13.3$°C, suggesting an implausible reduction in electricity demand under extreme cold temperatures. Similarly, in the extreme heat range, the benchmark response function from the hazard-to-demand model suggests a steeper increase in demand, surpassing the growth rate implied by that of the density-to-demand model.
While benchmark response functions from both models consistently show that cold waves result in greater demand increases than heat waves, that of the density-to-demand model exhibits excessive sensitivity, predicting disproportionately large demand fluctuations in response to imperceptibly small temperature variations (e.g., 0.1°C) under extreme cold conditions. In contrast, the benchmark response function from the hazard-to-demand model under extreme cold scenarios demonstrates more moderate sensitivity but still exhibits slightly larger responses compared to extreme heat scenarios. In this context, the hazard-to-demand model aligns more closely with empirical observations from historical heat and cold wave events.
In the reported estimation results for the functional shocks generated from the observed cold and heat waves (the right panels of Figures (ref) and (ref)), it is important to note that both the density-to-demand and hazard-to-demand models capture the incremental effects of extreme heat on electricity demand, whereas the benchmark response function \( \widetilde{\psi}(s) \) from either of the models fails to capture such effects as temperature approaches extreme heat. This suggests that distributional characteristics associated with heat wave events, such as a more right-skewed distributional change with significantly reduced kurtosis, further amplify electricity demand in a nonlinear manner. In contrast, incorporating distributional characteristics of historical cold wave events, such as significantly increased variance, into the benchmark response function does not qualitatively affect its overall shape.
Based on the derived implications above, the hazard-to-demand model appears to provide an empirically more sensible framework for analyzing the impact of extreme temperature events on electricity demand. This may be because, unlike the PDF, which considers the entire temperature distribution and may consequently dilute the influence of rare but critical extreme events that occur as tail events, the LHR predictor assigns substantial importance to tail events, as is evident from its construction. This allows for a more targeted analysis, focusing on the likelihood of temperatures persisting within an extreme range given that they have already exceeded a critical threshold (e.g., -15°C or 35°C). By emphasizing tail events, the LHR approach enhances the reliability of demand estimations under extreme temperature conditions.
An essential factor in extreme temperature analysis would be the cumulative temperature effect, wherein prolonged exposure to high or low temperatures significantly increases electricity demand for cooling or heating, even if the actual temperature remains unchanged. This phenomenon results in an apparent paradox where electricity demand continues to rise despite stable temperature levels, as extended exposure amplifies the need for indoor climate control. For instance, during a heatwave, electricity demand increases not only due to immediate high temperatures but also because prolonged exposure intensifies cooling requirements. In this context, the PDF, which captures only the relative frequency of temperatures without considering exposure duration, may exhibit weaker comovement with electricity demand, particularly under extreme conditions. Conversely, the LHR framework, which inherently accounts for the persistence of extreme temperature states, is expected to provide a more robust methodological foundation for assessing the dynamic relationship between extreme temperature events and electricity demand.
It might be of interest to consider quantile function corresponding to $\phi_t$ as the distributional predictor in our DSR model. The model with the quantile function predictor shares a conceptual similarity with the density-to-demand and hazard-to-demand models in that all three utilize the distributional properties of temperature as predictive variables. However, these models differ fundamentally in their interpretation of the demand response function \( f \). In the density-to-demand and hazard-to-demand models, the kernel-weighted benchmark response function \( \widetilde{\psi}(s) \) measures the impact of more frequently occurring temperatures around a given temperature level \( s \in [-20,40] \), either unconditionally or conditionally, while allowing for potential variations in associated temperature quantiles. In contrast, a similar quantity, which can similarly be defined in the quantile-to-demand model, represents the effect of changes around the \( s \)-quantile, where \( s \in [0,1] \), with the corresponding temperature level varying due to monthly (seasonal) fluctuations.
To illustrate the difference between the quantile-to-demand model and the previously considered models, we additionally examine the impact of temperature quantiles on electricity demand during extreme weather events by estimating the quantile-to-demand model, a more detailed description of which is provided in Section (ref) of the Appendix. Notably, the benchmark response function from the quantile-to-demand model, reported in Figure (ref), remains largely insignificant, except within the 0.63 to 0.87 quantile range for $\kappa=1$. This contrasts with the temperature density and LHR sensitivity presented in the left panels of Figures (ref) and (ref). The temperature quantile time series exhibits a monthly pattern, with lower values in winter and higher values in summer, whereas residential electricity demand exhibits a dual-peak structure, rising sharply in both summer and winter. This divergence introduces nonlinearities in the relationship between temperature and electricity demand, as an increase in temperature quantile time series simultaneously reduces heating demand while increasing cooling demand. Consequently, unlike an increase in the temperature density or LHR time series, which is unambiguously associated with higher cooling (heating) demand around 40°C (-20°C), applying fixed-quantile shocks in the quantile-to-demand framework would result in ambiguous slope estimates due to these offsetting effects. Similarly, incorporating the distributional characteristics into the quantile shocks in Figure (ref) would not change the result of producing insignificant or ambiguous temperature sensitivity.
In this paper, we study the nonlinear temperature sensitivity of residential electricity demand using distributional predictors of temperature, which are likely to contain measurement errors. Given that residential electricity demand is highly sensitive to temperature fluctuations, accurately capturing these nonlinearities is critical, especially during extreme temperature events where underestimation could jeopardize grid stability. Additionally, measurement errors may stem from the assumption that nonparametrically estimated temperature distributions represent true distributions, compounded by data quality issues and limitations within the estimation process.
We address these challenges by applying our DSR framework with distributional predictors. This approach allows for the effective identification of nonlinearities in the relationship between temperature distribution and electricity demand, even in the presence of measurement errors. Through the analysis of historical heat wave and cold wave shocks of varying intensities within a distributional context, we find that the hazard-to-demand model more reliably captures nonlinear temperature sensitivity than the density-to-demand model, particularly in extreme conditions. Our results indicate that residential electricity demand is significantly more sensitive to intense cold waves, with demand increasing sharply as temperatures decline. However, as cold wave severity intensifies, the rate of electricity demand growth diminishes, whereas for heat waves, demand growth accelerates with increasing intensity.
Given our findings, energy practitioners should prioritize incorporating detailed temperature distribution analysis, particularly when planning for extreme weather events. Our empirical findings suggest that residential electricity demand exhibits greater sensitivity to cold waves compared to heat waves, albeit with a slower rate of increase in demand, highlighting the need for advanced demand estimating models that account for these nonlinear sensitivities. Moreover, it is crucial to address potential measurement errors by adopting robust econometric techniques, in order to ensure that demand response strategies are based on accurate and reliable estimates. Our econometric approach would enhance grid stability and support effective energy management during extreme temperature events.